Numerics with Infinities


Aubrey Jaffer

This SRFI is currently in ``draft'' status. To see an explanation of each status that a SRFI can hold, see here. It will remain in draft status until 2005/07/13, or as amended. To provide input on this SRFI, please mailto:srfi-70@srfi.schemers.org. See instructions here to subscribe to the list. You can access previous messages via the archive of the mailing list.


When Scheme was an experimental and research language, the Scheme reports were written sparsely to provide a framework which allowed unusual implementations. With many experiments having been tried and hopefully having gained knowledge from them, and with a growing user base wanting portability for their programs, this SRFI attempts to tighten the specifications for Scheme's numeric subsystems.

Checks for division by zero (and their absence) remain one of the most common programming errors. In numerical code the conscientious programmer must deeply nest conditionals to assure that checks for zero preceed each division.

With infinities added to the number system we find that division by zero "works". It lets initialization of variables precede bounds checks and gives flexibility in placement of those checks. Limit and corner cases need not crash programs.

This SRFI reworks section 6.2 "Numbers" of R5RS to:




This SRFI strengthens the definition of exact numbers by constraining each exact number to represent a single mathematical number. With that proviso, exact number systems have limits in how precise a number can be represented in a finite amount of storage. Operations which would return too precise exact numbers can either return an inexact number (assuming the range of inexacts is sufficient) or report a violation of an implementation restriction.

Note that the iterative series often used in floating-point calculations cannot be summarily converted to work on exact numbers without the danger of exceeding the implementation limits just mentioned.


The definition of inexact numbers is strengthened by delineating their approximate nature:
Every mathematical number within the (convex) range of inexacts supported by an implementation will round to an inexact number on input or as a result of computation. The neighborhood of mathematical numbers rounding to a particular inexact number must be simply connected.
With this clear distinction between exact and inexact numbers, the last sentence of section 6.2.1 Numerical types, which seems to imply that each internal number representation is instantiated both in exact and inexact forms, is struck:
... In order to catch uses of inexact numbers where exact numbers are required, Scheme explicitly distinguishes exact from inexact numbers. This distinction is orthogonal to the dimension of type.
Section 6.2.2 Exactness states:
A number is inexact if it was written as an inexact constant, if it was derived using inexact ingredients, or if it was derived using inexact operations. Thus inexactness is a contagious property of a number.
But this section also states:
An operation may, however, return an exact result if it can prove that the value of the result is unaffected by the inexactness of its arguments. For example, multiplication of any number by an exact zero may produce an exact zero result, even if the other argument is inexact.
By those sentences, inexactness is a contagious property of all numbers except 0, the only number whose exactness is contagious! To resolve this conflict, the latter two sentences are struck.

Integer versus Exact

R5RS describes quotient, remainder, and modulo as integer operations. They are required to work for integers, whether exact or inexact. The use of N for names of their arguments similarly declares gcd, lcm, odd?, and even? to be integer operations.

But these operations become poorly behaved for inexacts outside the range where each integer has a unique representation. For instance, all large magnitude IEEE-754 numbers are integers and even (a multiple of 2) because the radix, and hence the scaling factors, are even.

Distinguishing these functions on the basis of integer arguments overlooks their generalizations to exact rational numbers:

    modulo(x/y, w/z) = modulo(lcm(y, z)*x/y, lcm(y, z)*w/z)/lcm(y, z)
    remainder(x/y, w/z) = remainder(lcm(y, z)*x/y, lcm(y, z)*w/z)/lcm(y, z)
    gcd(x/y, w/z) = gcd(x, w)/lcm(y, z)
    lcm(x/y, w/z) = lcm(x, w)/gcd(y, z)
Common-Lisp mod and rem are general to rational and real arguments:
(mod 2/3 1/5)               ==>  1/15
(mod .666 1/5)              ==>  0.065999985
This SRFI extends quotient, remainder, and modulo to work for exact rationals and inexact reals.

In inexact calculations, one is not guaranteed that results will yield integers. Thus odd? and even? will sometimes signal errors or give false negatives when inexact integer calculations are slightly wide of their integer marks. If these predicates instead accept only exact integers, then their behavior is predictable; and conversion from inexacts to exacts must be explicit.

Convenience Functions

Grepping through SLIB and several Scheme applications' source for occurrences of `round', `ceiling', `floor', and `truncate' finds that in essentially all 30+ cases `inexact->exact' was called with the results of these functions.

One possibility is having those functions also convert to exact, like the Common-Lisp functions of the same names do. But implementation of the extension of `quotient', `modulo', and `remainder' to finite real numbers requires inexact `floor' and `truncate'.

To capture these common patterns of usage, `round->exact', `ceiling->exact', `floor->exact', and `truncate->exact' are added.


This SRFI introduces inexact real infinities for implementations with inexact numbers. These infinities arise as the result of: This SRFI attempts to strike a balance between tolerance and signaling of range errors. Having +/0. and -/0. allows one level of calculation to be executed before bounds checking. Subsequent numerical calculations on infinities can signal an error or return a non-real infinity: 0/0.. 0/0. can be considered an error waiting to happen.

This SRFI mandates the real infinities for implementations with inexacts, but leaves to the implementations the decision of whether to implement 0/0. or not.

Are there predicates for infinities?
The `infinite?' predicate will return `#t' for infinities, `#f' for other numbers;
The `finite?' predicate will return `#f' for infinities, `#t' for other numbers;
`real?' will return `#f' for `0/0.'; and
the two real infinities have opposite signs.

What about exact infinities?
+/0. + 1 = +/0.; the only numbers which can be their own increments are inexact numbers.

6.2.2 Exactness gives semantics for inexact numbers, which correspond to mathematical neighborhoods. Its standard analysis with (mathematical) limits.

In contrast, each exact number corresponds to exactly one mathematical number. An exact infinity would have no additive inverse and a non-unique multiplicative inverse (0). If one gives the negative infinity an infinitesimal reciprocal, then the field is not closed under addition or multiplication.

How do infinities behave in comparisons?

(= +/0. +/0.)               ==>  #t
(= -/0. +/0.)               ==>  #f
(= -/0. -/0.)               ==>  #t
For any finite real number x:
(= -/0. x))                 ==>  #f
(= +/0. x))                 ==>  #f
(< -/0. x +/0.))            ==>  #t
(> +/0. x -/0.))            ==>  #t

Are real infinities rational?
No; (mathematical) rational numbers have non-zero denominators.

Are real infinities integers?
Because they are not rational they cannot be integers.

What about complex (non-real) infinities?
Infinities in rectangular notation leads to strange beasts like -/0.+5i. I know of no mathematical system having such numbers. Polar notation allows an infinitude of numbers having infinite magnitude and finite angle. Multiplication and division of complex infinities would work, but at the price of forcing their representations to be polar.

Chongkai Zhu points out that -/0. cannot be encoded by polar notation unless pi is a number in the implementation.

What about +0.0 and -0.0?
Dealing with zero is certainly thorny; but having more varieties of zero does not make dealing with unsigned zero any easier. The limit procedure provides a mathematical alternative to non-standard analysis half-measures.

0/0. = 0/0.

The bad boy of IEEE Standard 754 is not-a-number, which does not equal itself. This was a cute trick to provide a means for software to detect it without requiring literal constants or new comparison primitives.

But reflexive identity is one of the most basic mathematical principals. In implementations where inexact numbers are boxed

    (let ((narn (/ 0.0 0.0))) (eq? narn narn))
is required to be true, which forces eqv? to return true in the same situation:
Eq?'s behavior on numbers and characters is implementation-dependent, but it will always return either true or false, and will return true only when eqv? would also return true.
The `=' procedure is similarly forced to find narn equal to itself.


The description of sqrt has been wrong since R3RS. Prepending "For real z" to the second sentence corrects it.


In R3RS and R4RS the description of expt is very simple:
procedure: expt z1 z2

Returns z1 raised to the power z2:

z_1z_2 = ez_2 log z_1

00 is defined to be equal to 1.

The simplicity is nice, but in the absence of infinities, the straightforward scheme implementation
    (define (expt b x) (exp (* x (log b))))
is incapable of operating when b is zero.

According to G. Sussman, the change in R5RS to

0z is 1 if z = 0 and 0 otherwise.
was "... attempting to follow Common Lisp and APL when it came to boundary conditions and branch cuts."

About expt, Common-Lisp says:

expt is defined as bx = ex log b. This defines the principal values precisely. The range of expt is the entire complex plane. Regarded as a function of x, with b fixed, there is no branch cut. Regarded as a function of b, with x fixed, there is in general a branch cut along the negative real axis, continuous with quadrant II. The domain excludes the origin. By definition, 00=1. If b=0 and the real part of x is strictly positive, then bx=0. For all other values of x, 0x is an error.
With log(0) being -/0., and because exp(-/0.) is 0.0, Common-Lisp behavior would follow from:
    (define (expt b x) (exp (* x (log b))))
except when b is zero.

With (log 0) being -/0., the product of x and -/0. determines the trajectory of exp. Exp's domain repeats every 2*pi radians in the imaginary direction. This multiplication rotates the negative real infinity. If it rotates into a direction other than pure real, then the limit of (exp (* b (log x))) does not exist (because of the 2*pi periodicity).

Common-Lisp's criterion of testing the sign of the real part of x can be seen as snapping the direction of the product to either a positive or negative real infinity.

(exp (* 1 -/0.))            ==>  0.0
(exp (* -1 -/0.))           ==>  +/0.
Because Common-Lisp doesn't have explicit infinities, it signals an error for negative real parts of x. The Common-Lisp float behavior except, for (expt 0.0 0.0), can be produced by
(define (expt z1 z2) (exp (* (if (zero? z1) (real-part z2) z2) (log z1))))

Having real infinities enables these definitions to work when z1 is zero. Although the simpler definition is mathematically satisfying, the "real-part" variant allows (expt 0.0 z) to return something other than 0/0. when the exponent has small imaginary parts due to inexact inaccuracies.

(define (expt z1 z2)
  (cond ((and (exact? z2) (not (and (zero? z1) (negative? z2))))
	 (integer-expt z1 z2))
	((zero? z2) (+ 1 (* z1 z2)))
	 (exp (* (if (zero? z1) (real-part z2) z2) (log z1))))))
In the (above) definition from the reference implementation, (expt z 0.0) returns 1 or 1.0 when z is finite; and 0/0. otherwise.

(expt 0.0 z), where z has a negative real-part, returns +/0., where Common-Lisp signals an error. Otherwise the behavior is the same as Common-Lisp's expt.

Inexactness Considered Harmful?

Sebastian Egner, Richard Kelsey, Michael Sperber:
Cleaning up the Tower: Numbers in Scheme,
The 2004 Scheme Workshop, Snowbird, Utah, October 2004.


The R5RS specification of numerical operations leads to unportable and intransparent behavior of programs. Specifically, the notion of "exact/inexact numbers" and the misleading distinction between "real" and "rational" numbers are two primary sources of confusion. Consequently, the way R5RS organizes numbers is significantly less useful than it could be. Based on this diagnosis, we propose to abandon the concept of exact/inexact numbers from Scheme altogether. In this paper, we examine designs in which exact and inexact rounding operations are explicitly separated, while there is no distinction between exact and inexact numbers. Through examining alternatives and practical ramifications, we arrive at an alternative proposal for the design of the numerical operations in Scheme.
In my view the essential distinction of number-theoretic operations is that they are exact, not that they are integer as R5RS holds. This SRFI rectifies the problem by extending most of the affected procedures. Restricting the integer operations to exact integer arguments would also work.

Cleaning up the Tower demotes the exactness property relative to types while this SRFI promotes it. Yet both find that the exactness of procedures should play a more important role in determining the exactness of numerical results.

In their proposed system:

Any program that does not contain calls to floating-point operations always computes exactly and reproducibly, independent of the Scheme implementation it runs on.
This will require either As an author of a Scheme implementation, I find the first option unappealing; coding a complicated feature to solve a problem that the implementation doesn't have. If some feature must be jettisoned, then removing exact rational non-integers might be a simpler and less wrenching alternative than removing exactness as an organizing principle. Experience with SCM shows that the R5RS number system can work well without exact rational non-integers.

Full reproducibility will require extensive changes to R5RS:

The reproducibility condition is too strong; it needs restrictions so that it doesn't expose properties of the host file system.

Clean Infinities

The treatment of infinities by "Cleaning up the Tower" also differs from this SRFI.
... an exact division by zero is virtually always a symptom of a genuine programming error or of illegal input data, and the introduction of infinity will only mask this error.
This is not the case in the R5RS model, where substitution of exact for inexact does not change the computation.

The advantage of returning NaN instead of raising an error is that the computation still continues, postponing the interpretation of the results to a more convenient point in the program. In this way, NaN is quite useful in numerical computations.
While finding +/0. and -/0. to be very useful in computation, I cannot say the same for 0/0.. It usually arises as the result of operating on a real infinity. Tolerating one error is useful; tolerating more than one error in a computation masks programming errors. That is why this SRFI leaves to the implementation whether to signal an error or return 0/0..
The problem with NaN is that the program control structure will mostly not recognize the NaN case explicitly. Assume we define comparisons with NaN always to result in #f, as IEEE 754 does, then
(do ((x NaN (+ x 1))) ((> x 10)))
will hang but
(do ((x NaN (+ x 1))) ((not (<= x 10))))
will stop, which is counter-intuitive and may be surprising.
Because it has no sensible place in the total-order of real numbers, 0/0. is not a real number. Thus it is an illegal argument to the comparison procedures `<', `<=', `>', and `>='.
While +inf, -inf, and NaN are quite useful for inexact computations, there is a high price to pay when they are carried over into the exact world: The rational numbers must be extended by the special objects, and the usual algebraic laws will not hold for the extension anymore.
This need not be the case. This "high price" is a consequence of the proposed (section 5) "type permeability" kinds #2 and #3, which prohibit exact arguments to inexact procedures.

As shown earlier, infinities are always inexact. Infinities being inexact makes detecting them in exact calculations easier. Not only will the `finite?', `rational?', and `exact?' procedures return `#f' for infinities; but passing infinities to exact-only or integer-only procedures is an error.

The Limit

Many Scheme implementations treat an exact 0 as stronger than an inexact 0.0; returning 0 for (/ 0 0.0) while signaling an error for (/ 0.0 0). This leads to the use of exactness coercions to select behavior of `/' and `expt'; a practice which is very opaque.

This SRFI specifies that division of zeros be 0/0. or signal an error so long as one of the zeros is inexact. But why be restricted to that single choice? This SRFI introduces the `limit' procedure to produce the mathematically expected behavior of finite and infinite one-sided limits of continuous functions. Its mathematical basis and algorithm are detailed in the article Computing the Limit.


Here is the text proposed to replace section "6.2 Numbers" of R5RS. Deleted text is marked with a line through it. Additions and changes are marked in red.
Note: The type restriction for the naming convention `r' is "exact rational number":
r, r1, ... rj, ...
exact rational number

6.2 Numbers

Numerical computation has traditionally been neglected by the Lisp community. Until Common Lisp there was no carefully thought out strategy for organizing numerical computation, and with the exception of the MacLisp system [Pitman83] little effort was made to execute numerical code efficiently. This report recognizes the excellent work of the Common Lisp committee and accepts many of their recommendations. In some ways this report simplifies and generalizes their proposals in a manner consistent with the purposes of Scheme.

It is important to distinguish between the mathematical numbers, the Scheme numbers that attempt to model them, the machine representations used to implement the Scheme numbers, and notations used to write numbers. This report uses the types number, complex, real, rational, and integer to refer to both mathematical numbers and Scheme numbers. Machine representations such as fixed point and floating point are referred to by names such as fixnum and flonum.

6.2.1 Numerical types

Mathematically, numbers may be arranged into a tower of subtypes in which each level is a subset of the level above it:


For example, 3 is an integer. Therefore 3 is also a rational, a real, and a complex. The same is true of the Scheme numbers that model 3. For Scheme numbers, these types are defined by the predicates number?, complex?, real?, rational?, and integer?.

There is no simple relationship between a number's type and its representation inside a computer. Although most implementations of Scheme will offer at least two different representations of 3, these different representations denote the same integer.

Scheme's numerical operations treat numbers as abstract data, as independent of their representation as possible. Although an implementation of Scheme may use fixnum, flonum, and perhaps other representations for numbers, this should not be apparent to a casual programmer writing simple programs.

It is necessary, however, to distinguish between numbers that are represented exactly and those that may not be. For example, indexes into data structures must be known exactly, as must some polynomial coefficients in a symbolic algebra system. On the other hand, the results of measurements are inherently inexact, and irrational numbers may be approximated by rational and therefore inexact approximations. In order to catch uses of inexact numbers where exact numbers are required, Scheme explicitly distinguishes exact from inexact numbers. This distinction is orthogonal to the dimension of type.

6.2.2 Exactness

Scheme numbers are either exact or inexact. A number is exact if it was written as an exact constant or was derived from exact numbers using only exact operations. A number is inexact if it is infinite, if it was written as an inexact constant, if it was derived using inexact ingredients, or if it was derived using inexact operations. Thus inexactness is a contagious property of a number.

Each exact number corresponds to a single mathematical number. It is the programmer's responsibility to avoid using exact numbers with magnitude or precision too large to be represented in the implementation.

Inexact numbers are approximate. Every mathematical number within the (convex) range of inexacts supported by an implementation will round to an inexact number on input or as a result of computation. The neighborhood of mathematical numbers rounding to a particular inexact number must be simply connected.

Because real infinities are mandated (in implementations supporting inexacts), all mathematical real numbers map to inexact numbers. It is the programmer's responsibility to avoid using non-real complex numbers with magnitude too large to be represented in the implementation.

If two implementations produce exact results for a computation that did not involve inexact intermediate results, the two ultimate results will be mathematically equivalent. This is generally not true of computations involving inexact numbers since because approximate methods such as floating point arithmetic may be used, but it is the duty of each implementation to make the result as close as practical to the mathematically ideal result.

Rational operations such as `+' should always produce exact results when given exact arguments. If the operation is unable to produce an exact result, then it may either report the violation of an implementation restriction or it may silently coerce its result to an inexact value. See section 6.2.3 Implementation restrictions.

With the exception of inexact->exact, `round->exact', `ceiling->exact', `floor->exact', and `truncate->exact', the operations described in this section must generally return inexact results when given any inexact arguments. An operation may, however, return an exact result if it can prove that the value of the result is unaffected by the inexactness of its arguments. For example, multiplication of any number by an exact zero may produce an exact zero result, even if the other argument is inexact.

6.2.2x Infinities

The interpretation of real infinities is that +/0. represents real numbers greater than can be encoded by finite inexacts in the implementation (> 179.76931348623157e306 for IEEE-754 64-bit flonums) and that -/0. represents numbers less than can be encoded by finite inexacts in the implementation (< -179.76931348623157e306 for IEEE-754 64-bit flonums). This preserves the total ordering of the (mathematical) real numbers and extends Scheme's representation to cover the entire real line. +/0. and -/0. represent the half-lines beyond either end of the implementation's inexact rational range. For any finite real number x:
(= -/0. x))                 ==>  #f
(= +/0. x))                 ==>  #f
(< -/0. x +/0.))            ==>  #t
(> +/0. x -/0.))            ==>  #t

Implementations of Scheme which provide inexact real numbers shall implement positive infinity and negative infinity as unique inexact real numbers.

An optional third infinity, which is not real, may be returned by a numerical function when no inexact number (including infinities) is the correct answer. An implementation may report a violation of an implementation restriction in any calculation for which the result would be an unreal infinity.

6.2.3 Implementation restrictions

Implementations of Scheme are not required to implement the whole tower of subtypes given in section 6.2.1 Numerical types, but they must implement a coherent subset consistent with both the purposes of the implementation and the spirit of the Scheme language. For example, an implementation in which all numbers are real may still be quite useful.

Implementations may also support only a limited range of numbers of any type, subject to the requirements of this section. The supported range for exact numbers of any type may be different from the supported range for inexact numbers of that type. For example, an implementation that uses flonums to represent all its inexact real numbers may support a practically unbounded range of exact integers and rationals while limiting the range of inexact reals (and therefore the range of inexact integers and rationals) to the dynamic range of the flonum format. Furthermore the gaps between the representable inexact integers and rationals are likely to be very large in such an implementation as the limits of this range are approached.

An implementation of Scheme must support exact integers throughout the range of numbers that may be used for indexes of lists, vectors, and strings or that may result from computing the length of a list, vector, or string. The length, vector-length, and string-length procedures must return an exact integer, and it is an error to use anything but an exact integer as an index. Furthermore any integer constant within the index range, if expressed by an exact integer syntax, will indeed be read as an exact integer, regardless of any implementation restrictions that may apply outside this range. Finally, the procedures listed below will always return an exact integer result provided all their arguments are exact integers and the mathematically expected result is representable as an exact integer within the implementation:

+            -             *
quotient     remainder     modulo
max          min           abs
numerator    denominator   gcd
lcm          floor         ceiling
truncate     round         rationalize

Implementations are encouraged, but not required, to support exact integers and exact rationals of practically unlimited size and precision, and to implement the above procedures and the `/' procedure in such a way that they always return exact results when given exact arguments. If one of these procedures is unable to deliver an exact result when given exact arguments, then it may either report a violation of an implementation restriction or it may silently coerce its result to an inexact number. Such a coercion may cause an error later.

An implementation may use floating point and other approximate representation strategies for inexact numbers.

This report recommends, but does not require, that the IEEE 32-bit and 64-bit floating point standards be followed by implementations that use flonum representations, and that implementations using other representations should match or exceed the precision achievable using these floating point standards [IEEE].

In particular, implementations that use flonum representations must follow these rules: A flonum result must be represented with at least as much precision as is used to express any of the inexact arguments to that operation. It is desirable (but not required) for potentially inexact operations such as `sqrt', when applied to exact arguments, to produce exact answers whenever possible (for example the square root of an exact 4 ought to be an exact 2). If, however, an exact number is operated upon so as to produce an inexact result (as by `sqrt'), and if the result is represented as a flonum, then the most precise flonum format available must be used; but if the result is represented in some other way then the representation must have at least as much precision as the most precise flonum format available.

Although Scheme allows a variety of written notations for numbers, any particular implementation may support only some of them. For example, an implementation in which all numbers are real need not support the rectangular and polar notations for complex numbers. If an implementation encounters an exact numerical constant that it cannot represent as an exact number, then it may either report a violation of an implementation restriction or it may silently represent the constant by an inexact number.

6.2.4 Syntax of numerical constants

The syntax of the written representations for numbers is described formally in section 7.1.1 Lexical structure. Note that case is not significant in numerical constants.

A number may be written in binary, octal, decimal, or hexadecimal by the use of a radix prefix. The radix prefixes are `#b' (binary), `#o' (octal), `#d' (decimal), and `#x' (hexadecimal). With no radix prefix, a number is assumed to be expressed in decimal.

A numerical constant may be specified to be either exact or inexact by a prefix. The prefixes are `#e' for exact, and `#i' for inexact. An exactness prefix may appear before or after any radix prefix that is used. If the written representation of a number has no exactness prefix, the constant may be either inexact or exact. It is inexact if it contains a decimal point, an exponent, a "#" character in the place of a digit; otherwise it is exact.

Negative infinity is written `-/0.'. Positive infinity is written `+/0.'. The (optional) non-real infinity is written `0/0.'.

In systems with inexact numbers of varying precisions it may be useful to specify the precision of a constant. For this purpose, numerical constants may be written with an exponent marker that indicates the desired precision of the inexact representation. The letters `s', `f', `d', and `l' specify the use of short, single, double, and long precision, respectively. (When fewer than four internal inexact representations exist, the four size specifications are mapped onto those available. For example, an implementation with two internal representations may map short and single together and long and double together.) In addition, the exponent marker `e' specifies the default precision for the implementation. The default precision has at least as much precision as double, but implementations may wish to allow this default to be set by the user.

       Round to single --- 3.141593
       Extend to long --- .600000000000000

6.2.5 Numerical operations

The reader is referred to section 1.3.3 Entry format for a summary of the naming conventions used to specify restrictions on the types of arguments to numerical routines.

The examples used in this section assume that any numerical constant written using an exact notation is indeed represented as an exact number. Some examples also assume that certain numerical constants written using an inexact notation can be represented without loss of accuracy; the inexact constants were chosen so that this is likely to be true in implementations that use flonums to represent inexact numbers.

procedure: number? obj
procedure: complex? obj
procedure: real? obj
procedure: rational? obj
procedure: integer? obj

These numerical type predicates can be applied to any kind of argument, including non-numbers. They return #t if the object is of the named type, and otherwise they return #f. In general, if a type predicate is true of a number then all higher type predicates are also true of that number. Consequently, if a type predicate is false of a number, then all lower type predicates are also false of that number.

If z is an inexact complex number, then `(real? z)' is true if and only if `(zero? (imag-part z))' is true. If x is an inexact real number, then `(integer? x)' is true if and only if

    (and (finite? x) (= x (round x)))
(complex? 3+4i)                        ==>  #t
(complex? 3)                           ==>  #t
(real? 3)                              ==>  #t
(real? -2.5+0.0i)                      ==>  #t
(real? #e1e10)                         ==>  #t
(rational? 6/10)                       ==>  #t
(rational? 6/3)                        ==>  #t
(integer? 3+0i)                        ==>  #t
(integer? 3.0)                         ==>  #t
(integer? 8/4)                         ==>  #t
(number? 0/0.)                         ==>  #t
(complex? 0/0.)                        ==>  #f
(complex? +/0.)                        ==>  #t
(real? 0/0.)                           ==>  #f
(real? -/0.)                           ==>  #t
(rational? +/0.)                       ==>  #f
(rational? 0/0.)                       ==>  #f
(integer? -/0.)                        ==>  #f

Note: The behavior of these type predicates on inexact numbers is unreliable, since because any inaccuracy may affect the result.

Note: In many implementations the rational? procedure will be the same as real?, and the complex? procedure will be the same as number?, but unusual implementations may be able to represent some irrational numbers exactly or may extend the number system to support some kind of non-complex numbers.

procedure: exact? z
procedure: inexact? z

These numerical predicates provide tests for the exactness of a quantity. For any Scheme number, precisely one of these predicates is true.

(exact? 5)                   ==>  #t
(inexact? +/0.)              ==>  #t
(inexact? 0/0.)              ==>  #t

procedure: = z1 z2 z3 ...
procedure: < x1 x2 x3 ...
procedure: > x1 x2 x3 ...
procedure: <= x1 x2 x3 ...
procedure: >= x1 x2 x3 ...

These procedures return #t if their arguments are (respectively): equal, monotonically increasing, monotonically decreasing, monotonically nondecreasing, or monotonically nonincreasing.

(= +/0. +/0.)               ==>  #t
(= -/0. +/0.)               ==>  #f
(= -/0. -/0.)               ==>  #t
(= 0/0. 0/0.)               ==>  #t
For any finite real number x:
(< -/0. x +/0.))            ==>  #t
(> +/0. x -/0.))            ==>  #t

These predicates are required to be transitive.

Note: The traditional implementations of these predicates in Lisp-like languages are not transitive.

Note: While it is not an error to compare inexact numbers using these predicates, the results may be unreliable because a small inaccuracy may affect the result; this is especially true of = and zero?. When in doubt, consult a numerical analyst.

library procedure: finite? z
library procedure: infinite? z
library procedure: zero? z
library procedure: positive? x
library procedure: negative? x
library procedure: odd? n
library procedure: even? n

These numerical predicates test a number for a particular property, returning #t or #f. See note above.

(positive? +/0.)            ==>  #t
(negative? -/0.)            ==>  #t
(finite? -/0.)              ==>  #f
(infinite? +/0.)            ==>  #t
(finite? 0/0.)              ==>  #f
(infinite? 0/0.)            ==>  #t

library procedure: max x1 x2 ...
library procedure: min x1 x2 ...

These procedures return the maximum or minimum of their arguments.

(max 3 4)                              ==>  4    ; exact
(max 3.9 4)                            ==>  4.0  ; inexact
For any real number x:
(max +/0. x)                           ==>  +/0.
(min -/0. x)                           ==>  -/0.

Note: If any argument is inexact, then the result will also be inexact (unless the procedure can prove that the inaccuracy is not large enough to affect the result, which is possible only in unusual implementations). If `min' or `max' is used to compare numbers of mixed exactness, and the numerical value of the result cannot be represented as an inexact number without loss of accuracy, then the procedure may report a violation of an implementation restriction.

procedure: + z1 ...
procedure: * z1 ...

These procedures return the sum or product of their arguments.

(+ 3 4)                                ==>  7
(+ 3)                                  ==>  3
(+)                                    ==>  0
(+ +/0. +/0.)                          ==>  +/0.
(+ +/0. -/0.)                          ==>  0/0.

(* 4)                                  ==>  4
(*)                                    ==>  1
(* 5 +/0.)                             ==>  +/0.
(* -5 +/0.)                            ==>  -/0.
(* +/0. +/0.)                          ==>  +/0.
(* +/0. -/0.)                          ==>  -/0.
(* 0 +/0.)                             ==>  0/0.
For any finite number z:
(+ +/0. z)                             ==>  +/0.
(+ -/0. z)                             ==>  -/0.
For any number z:
(+ 0/0. z)                             ==>  0/0.
(* 0/0. z)                             ==>  0/0.

procedure: - z1 z2
procedure: - z
optional procedure: - z1 z2 ...
procedure: / z1 z2
procedure: / z
optional procedure: / z1 z2 ...

With two or more arguments, these procedures return the difference or quotient of their arguments, associating to the left. With one argument, however, they return the additive or multiplicative inverse of their argument.

(- 3 4)                                ==>  -1
(- 3 4 5)                              ==>  -6
(- 3)                                  ==>  -3
(- +/0. +/0.)                          ==>  0/0.

(/ 3 4 5)                              ==>  3/20
(/ 3)                                  ==>  1/3
(/ 0.0)                                ==>  +/0.
(/ 1.0 0)                              ==>  +/0.
(/ -1 0.0)                             ==>  -/0.
(/ +/0.)                               ==>  0.0
(/ 0 0.0)                              ==>  0/0.
(/ 0.0 0)                              ==>  0/0.
(/ 0.0 0.0)                            ==>  0/0.

library procedure: abs x

`Abs' returns the absolute value of its argument.

(abs -7)                               ==>  7
(abs -/0.)                             ==>  +/0.

procedure: quotient n1 n2 x1 x2
procedure: remainder n1 n2 x1 x2
procedure: modulo n1 n2 x1 x2

These procedures implement number-theoretic (integer) division. x2 should be non-zero. All three procedures return integers. `quotient' returns an integer. If x1/x2 is an integer:

    (quotient x1 x2)                   ==> x1/x2
    (remainder x1 x2)                  ==> 0
    (modulo x1 x2)                     ==> 0

If x1/x2 is not an integer:

    (quotient x1 x2)                   ==> n_q
    (remainder x1 x2)                  ==> x_r
    (modulo x1 x2)                     ==> x_m

where n_q is x1/x2 rounded towards zero, 0 < |x_r| < |x2|, 0 < |x_m| < |x2|, x_r and x_m differ from x1 by a multiple of x2, x_r has the same sign as x1, and x_m has the same sign as x2.

From this we can conclude that for integers n1 and n2 with x2 not equal to 0,

     (= x1 (+ (* x2 (quotient x1 x2))
           (remainder x1 x2)))
                                       ==>  #t

provided all numbers involved in that computation are exact.

(modulo 13 4)                          ==>  1
(remainder 13 4)                       ==>  1

(modulo -13 4)                         ==>  3
(remainder -13 4)                      ==>  -1

(modulo 13 -4)                         ==>  -3
(remainder 13 -4)                      ==>  1

(modulo -13 -4)                        ==>  -1
(remainder -13 -4)                     ==>  -1

(remainder -13 -4.0)                   ==>  -1.0  ; inexact

(quotient 2/3 1/5)                     ==>  3
(modulo 2/3 1/5)                       ==>  1/15

(quotient .666 1/5)                    ==>  3
(modulo .666 1/5)                      ==>  65.99999999999995e-3

library procedure: gcd n1 r1 ...
library procedure: lcm n1 r1 ...

These procedures return the greatest common divisor or least common multiple of their arguments. The result is always non-negative.

For exact integer arguments, these procedures are the familiar number theoretic operators:

(gcd 32 -36)                           ==>  4
(gcd)                                  ==>  0
(lcm 32 -36)                           ==>  288
(lcm)                                  ==>  1
For exact rational arguments, gcd returns the largest rational that divides into each of its arguments a whole number of times, while lcm returns the smallest rational that is an integer multiple of its arguments.
(gcd 1/6 1/4)                          ==>  1/12
(lcm 1/6 1/4)                          ==>  1/2
(gcd 1/6 5/4)                          ==>  1/12
(lcm 1/6 5/4)                          ==>  5/2

procedure: numerator q
procedure: denominator q

These procedures return the numerator or denominator of their argument; the result is computed as if the argument was represented as a fraction in lowest terms. The denominator is always positive. The denominator of 0 is defined to be 1.

(numerator (/ 6 4))                    ==>  3
(denominator (/ 6 4))                  ==>  2
  (exact->inexact (/ 6 4)))            ==> 2.0

procedure: floor x
procedure: ceiling x
procedure: truncate x
procedure: round x

These procedures accept finite real numbers and return integers. `Floor' returns the largest integer not larger than x. `Ceiling' returns the smallest integer not smaller than x. `Truncate' returns the integer closest to x whose absolute value is not larger than the absolute value of x. `Round' returns the closest integer to x, rounding to even when x is halfway between two integers.

Rationale: `Round' rounds to even for consistency with the default rounding mode specified by the IEEE floating point standard.

Note: If the argument to one of these procedures is inexact, then the result will also be inexact. If an exact value is needed, the result should be passed to the `inexact->exact' procedure.

(floor -4.3)                           ==>  -5.0
(ceiling -4.3)                         ==>  -4.0
(truncate -4.3)                        ==>  -4.0
(round -4.3)                           ==>  -4.0

(floor 3.5)                            ==>  3.0
(ceiling 3.5)                          ==>  4.0
(truncate 3.5)                         ==>  3.0
(round 3.5)                            ==>  4.0  ; inexact

(round 7/2)                            ==>  4    ; exact
(round 7)                              ==>  7

library procedure: floor->exact x
library procedure: ceiling->exact x
library procedure: truncate->exact x
library procedure: round->exact x

These procedures are the compositions of `inexact->exact' with `floor', `ceiling', `truncate', and `round'.

library procedure: rationalize x y

`Rationalize' returns the simplest rational number differing from x by no more than y. A rational number r_1 is simpler than another rational number r_2 if r_1 = p_1/q_1 and r_2 = p_2/q_2 (in lowest terms) and |p_1|<= |p_2| and |q_1| <= |q_2|. Thus 3/5 is simpler than 4/7. Although not all rationals are comparable in this ordering (consider 2/7 and 3/5) any interval contains a rational number that is simpler than every other rational number in that interval (the simpler 2/5 lies between 2/7 and 3/5). Note that 0 = 0/1 is the simplest rational of all.

  (inexact->exact .3) 1/10)            ==> 1/3    ; exact
(rationalize .3 1/10)                  ==> #i1/3  ; inexact

(rationalize 3 +/0.)                   ==>  0

library procedure: limit proc x1 x2 k
library procedure: limit proc x1 x2

Proc must be a procedure taking a single inexact real argument. K is the number of points on which proc will be called; it defaults to 8.

If x1 is finite, then Proc must be continuous on the half-open interval:

( x1 .. x1+x2 ]

And x2 should be chosen small enough so that proc is expected to be monotonic or constant on arguments between x1 and x1 + x2.

Limit computes the limit of proc as its argument approaches x1 from x1 + x2. Limit returns a real number or real infinity or `#f'.

If x1 is not finite, then x2 must be a finite nonzero real with the same sign as x1; in which case limit returns:

(limit (lambda (x) (proc (/ x))) 0.0 (/ x2) k)

Limit examines the magnitudes of the differences between successive values returned by proc called with a succession of numbers from x1+x2/k to x1.

If the magnitudes of differences are monotonically decreasing, then then the limit is extrapolated from the degree n polynomial passing through the samples returned by proc.

If the magnitudes of differences are increasing as fast or faster than a hyperbola matching at x1+x2, then a real infinity with sign the same as the differences is returned.

If the magnitudes of differences are increasing more slowly than the hyperbola matching at x1+x2, then the limit is extrapolated from the quadratic passing through the three samples closest to x1.

If the magnitudes of differences are not monotonic or are not completely within one of the above categories, then #f is returned.

;; constant
(limit (lambda (x) (/ x x)) 0 1.0e-9)           ==> 1.0
(limit (lambda (x) (expt 0 x)) 0 1.0e-9)        ==> 0.0
(limit (lambda (x) (expt 0 x)) 0 -1.0e-9)       ==> +/0.
;; linear
(limit + 0 976.5625e-6)                         ==> 0.0
(limit - 0 976.5625e-6)                         ==> 0.0
;; vertical point of inflection
(limit sqrt 0 1.0e-18)                          ==> 0.0
(limit (lambda (x) (* x (log x))) 0 1.0e-9)     ==> -102.70578127633066e-12
(limit (lambda (x) (/ x (log x))) 0 1.0e-9)     ==> 96.12123142321669e-15
;; limits tending to infinity
(limit + +/0. 1.0e9)                            ==> +/0.
(limit + -/0. -1.0e9)                           ==> -/0.
(limit / 0 1.0e-9)                              ==> +/0.
(limit / 0 -1.0e-9)                             ==> -/0.
(limit (lambda (x) (/ (log x) x)) 0 1.0e-9)     ==> -/0.
(limit (lambda (x) (/ (magnitude (log x)) x)) 0 -1.0e-9)
                                                ==> -/0.
;; limit doesn't exist
(limit sin +/0. 1.0e9)                          ==> #f
(limit (lambda (x) (sin (/ x))) 0 1.0e-9)       ==> #f
(limit (lambda (x) (sin (/ x))) 0 -1.0e-9)      ==> #f
(limit (lambda (x) (/ (log x) x)) 0 -1.0e-9)    ==> #f
;; conditionally convergent - return #f
(limit (lambda (x) (/ (sin x) x)) +/0. 1.0e222)
                                                ==> #f
;; asymptotes
(limit / -/0. -1.0e222)                         ==> 0.0
(limit / +/0. 1.0e222)                          ==> 0.0
(limit (lambda (x) (expt x x)) 0 1.0e-18)       ==> 1.0
(limit (lambda (x) (sin (/ x))) +/0. 1.0e222)   ==> 0.0
(limit (lambda (x) (/ (+ (exp (/ x)) 1))) 0 1.0e-9)
                                                ==> 0.0
(limit (lambda (x) (/ (+ (exp (/ x)) 1))) 0 -1.0e-9)
                                                ==> 1.0
(limit (lambda (x) (real-part (expt (tan x) (cos x)))) (/ pi 2) 1.0e-9)
                                                ==> 1.0
;; This example from the 1979 Macsyma manual grows so rapidly
;;  that x2 must be less than 41.  It correctly returns e^2.
(limit (lambda (x) (expt (+ x (exp x) (exp (* 2 x))) (/ x))) +/0. 40)
                                                ==> 7.3890560989306504
;; LIMIT can calculate the proper answer when evaluation
;; of the function at the limit point does not:
(tan (atan +/0.))                               ==> 16.331778728383844e15
(limit tan (atan +/0.)-1.0e-15)                 ==> +/0.
(tan (atan +/0.))                               ==> 16.331778728383844e15
(limit tan (atan +/0.)1.0e-15)                  ==> -/0.
((lambda (x) (expt (exp (/ -1 x)) x)) 0)        ==> 1.0
(limit (lambda (x) (expt (exp (/ -1 x)) x)) 0 1.0e-9)
                                                ==> 0.0

procedure: exp z
procedure: log z
procedure: sin z
procedure: cos z
procedure: tan z
procedure: asin z
procedure: acos z
procedure: atan z
procedure: atan y x

These procedures are part of every implementation that supports general real numbers; they compute the usual transcendental functions. `Log' computes the natural logarithm of z (not the base ten logarithm). `Asin', `acos', and `atan' compute arcsine (sin-1), arccosine (cos-1), and arctangent (tan-1), respectively. The two-argument variant of `atan' computes (angle (make-rectangular x y)) (see below), even in implementations that don't support general complex numbers.

In general, the mathematical functions log, arcsine, arccosine, and arctangent are multiply defined. The value of log z is defined to be the one whose imaginary part lies in the range from -pi (exclusive) to pi (inclusive). log 0 is undefined. With log defined this way, The values of sin-1 z, cos-1 z, and tan-1 z are according to the following formulae:

sin-1 z = -i log (i z + sqrt(1 - z2))

cos-1 z = pi / 2 - sin-1 z

tan-1 z = (log (1 + i z) - log (1 - i z)) / (2 i)

The above specification follows [CLtL], which in turn cites [Penfield81]; refer to these sources for more detailed discussion of branch cuts, boundary conditions, and implementation of these functions. When it is possible these procedures produce a real result from a real argument.

If the function has a real-valued limit as its argument tends toward positive infinity, then that is the value returned by the function applied to +/0.. If the function has a real-valued limit as its argument tends toward negative infinity, then that is the value returned by the function applied to -/0..

(exp +/0.)                  ==> +/0.
(exp -/0.)                  ==> 0.0
(log +/0.)                  ==> +/0.
(log 0.0)                   ==> -/0.
(log -/0.)                  ==> 0/0.
(atan -/0.)                 ==> -1.5707963267948965
(atan +/0.)                 ==> 1.5707963267948965
The functions sin, cos, tan, asin, and acos either return 0/0. or report a violation of an implementation restriction when given +/0., -/0., or 0/0. as an argument.

procedure: sqrt z

Returns the principal square root of z. For real z The result will have either positive real part, or zero real part and non-negative imaginary part.

(sqrt -5)                   ==>  0.0+2.23606797749979i
(sqrt +/0.)                 ==>  +/0.
(sqrt -/0.)                 ==>  0/0.

procedure: expt z1 z2

Returns z1 raised to the power z2. For z_1 ~= 0

z_1^z_2 = e^z_2 log z_1

0^z is 1 if z = 0 and 0 otherwise. 0^0 is 1.

For inexact arguments not both zero

(define (expt z1 z2) (exp (* (if (zero? z1) (real-part z2) z2) (log z1))))
(expt 0.0 z)
returns 1.0 for z equal to 0;
returns 0.0 for z having positive real part (including +/0.);
returns +/0. for z having negative real part (including -/0.); and
returns 0/0. or reports a violation of an implementation restriction otherwise.
(expt 5 3)                  ==>  125
(expt 5 -3)                 ==>  1/125
(expt 5 0)                  ==>  1
(expt 0 5)                  ==>  0
(expt 0 0)                  ==>  1
(expt 0 5+.0000312i)        ==>  0.0
(expt 0 -5)                 ==>  +/0.
(expt 0 -5+.0000312i)       ==>  +/0.
(expt 0 0.0)                ==>  1.0

procedure: make-rectangular x1 x2
procedure: make-polar x3 x4
procedure: real-part z
procedure: imag-part z
procedure: magnitude z
procedure: angle z

These procedures are part of every implementation that supports general complex numbers. Suppose x1, x2, x3, and x4 are real numbers and z is a complex number such that

z = x1 + i x2 = x3 ei x4


(make-rectangular x1 x2)               ==> z
(make-polar x3 x4)                     ==> z
(real-part z)                          ==> x1
(imag-part z)                          ==> x2
(magnitude z)                          ==> |x3|
(angle z)                              ==> x_angle

where -pi < x_angle <= pi with x_angle = x4 + 2pi n for some integer n.

(angle +/0.)                ==> 0.0
(angle -/0.)                ==> 3.141592653589793

Rationale: `Magnitude' is the same as abs for a real argument, but `abs' must be present in all implementations, whereas `magnitude' need only be present in implementations that support general complex numbers.

procedure: exact->inexact z
procedure: inexact->exact z

`Exact->inexact' returns an inexact representation of z. The value returned is the inexact number that is numerically closest to the argument. If an exact argument has no reasonably close inexact equivalent, then a violation of an implementation restriction may be reported.

`Inexact->exact' returns an exact representation of z. The value returned is the exact number that is numerically closest to the argument. If an inexact argument has no reasonably close exact equivalent, then a violation of an implementation restriction may be reported.

These procedures implement the natural one-to-one correspondence between exact and inexact integers throughout an implementation-dependent range. See section 6.2.3 Implementation restrictions.

`Exact->inexact' and `inexact->exact' are idempotent.

6.2.6 Numerical input and output

procedure: number->string z
procedure: number->string z radix

Radix must be an exact integer, either 2, 8, 10, or 16. If omitted, radix defaults to 10. The procedure `number->string' takes a number and a radix and returns as a string an external representation of the given number in the given radix such that

(let ((number number)
      (radix radix))
  (eqv? number
        (string->number (number->string number

is true. It is an error if no possible result makes this expression true.

If z is inexact, the radix is 10, and the above expression can be satisfied by a result that contains a decimal point, then the result contains a decimal point and is expressed using the minimum number of digits (exclusive of exponent and trailing zeroes) needed to make the above expression true [howtoprint], [howtoread]; otherwise the format of the result is unspecified.

The result returned by `number->string' never contains an explicit radix prefix.

Note: The error case can occur only when z is not a complex number or is a complex number with a non-rational real or imaginary part.

Rationale: If z is an inexact number represented using flonums, and the radix is 10, then the above expression is normally satisfied by a result containing a decimal point. The unspecified case allows for infinities, NaNs, and non-flonum representations.

procedure: string->number string
procedure: string->number string radix

Returns a number of the maximally precise representation expressed by the given string. Radix must be an exact integer, either 2, 8, 10, or 16. If supplied, radix is a default radix that may be overridden by an explicit radix prefix in string (e.g. "#o177"). If radix is not supplied, then the default radix is 10. If string is not a syntactically valid notation for a number, then `string->number' returns #f.

(string->number "100")                 ==>  100
(string->number "100" 16)              ==>  256
(string->number "1e2")                 ==>  100.0
(string->number "15##")                ==>  1500.0
(string->number "+/0.")                ==>  +/0.
(string->number "-/0.")                ==>  -/0.

Note: The domain of `string->number' may be restricted by implementations in the following ways. `String->number' is permitted to return #f whenever string contains an explicit radix prefix. If all numbers supported by an implementation are real, then `string->number' is permitted to return #f whenever string uses the polar or rectangular notations for complex numbers. If all numbers are integers, then `string->number' may return #f whenever the fractional notation is used. If all numbers are exact, then `string->number' may return #f whenever an exponent marker or explicit exactness prefix is used, or if a # appears in place of a digit. If all inexact numbers are integers, then `string->number' may return #f whenever a decimal point is used. An implementation may return #f for "0/0.".


Most of the glibc transcendental functions do the right things when passed IEEE-754 infinities. Code which integrates glibc functions into Scheme procedures handling complex numbers can be found in scm/Transcen.scm.

Here is code for the procedures extended from R5RS:

(define (infinite? z) (and (= z (* 2 z)) (not (zero? z))))
(define (finite? z) (not (infinite? z)))

(define (ipow-by-squaring x n acc proc)
  (cond ((zero? n) acc)
	((eqv? 1 n) (proc acc x))
	(else (ipow-by-squaring (proc x x)
				(quotient n 2)
				(if (even? n) acc (proc acc x))

(define (integer-expt x n) (ipow-by-squaring x n (if (exact? x) 1 1.) *))

(define (expt z1 z2)
  (cond ((and (exact? z2) (not (and (zero? z1) (negative? z2))))
	 (integer-expt z1 z2))
	((zero? z2) (+ 1 (* z1 z2)))
	(else (exp (* (if (zero? z1) (real-part z2) z2) (log z1))))))

(define integer-quotient quotient)
(define integer-remainder remainder)
(define integer-modulo modulo)

(define (quotient x1 x2)
  (if (and (integer? x1) (integer? x2))
      (integer-quotient x1 x2)
      (truncate (/ x1 x2))))

(define (remainder x1 x2)
  (if (and (integer? x1) (integer? x2))
      (integer-remainder x1 x2)
      (- x1 (* x2 (quotient x1 x2)))))

(define (modulo x1 x2)
  (if (and (integer? x1) (integer? x2))
      (integer-modulo x1 x2)
      (- x1 (* x2 (floor (/ x1 x2))))))

(define integer-lcm lcm)
(define integer-gcd gcd)

(define (lcm . args)
  (/ (apply integer-lcm (map numerator args))
     (apply integer-gcd (map denominator args))))

(define (gcd . args)
  (/ (apply integer-gcd (map numerator args))
     (apply integer-lcm (map denominator args))))

(define (round->exact x) (inexact->exact (round x)))
(define (floor->exact x) (inexact->exact (floor x)))
(define (ceiling->exact x) (inexact->exact (ceiling x)))
(define (truncate->exact x) (inexact->exact (truncate x)))

(define (invintp f1 f2 f3)
  (define f1^2 (* f1 f1))
  (define f2^2 (* f2 f2))
  (define f3^2 (expt f3 2))
  (let ((c (+ (* -3 f1^2 f2)
	      (* 3 f1 f2^2)
	      (* (- (* 2 f1^2) f2^2) f3)
	      (* (- f2 (* 2 f1)) f3^2)))
	(b (+ (- f1^2 (* 2 f2^2)) f3^2))
	(a (- (* 2 f2) f1 f3)))
    (define disc (- (* b b) (* 4 a c)))
    (if (negative? (real-part disc))
	(/ b -2 a)
	(let ((sqrt-disc (sqrt disc)))
	  (define root+ (/ (- sqrt-disc b) 2 a))
	  (define root- (/ (+ sqrt-disc b) -2 a))
	  (if (< (magnitude (- root+ f1)) (magnitude (- root- f1)))

(define (extrapolate-0 fs)
  (define n (length fs))
  (define (choose n k)
    (do ((kdx 1 (+ 1 kdx))
	 (prd 1 (/ (* (- n kdx -1) prd) kdx)))
	((> kdx k) prd)))
  (do ((k 1 (+ 1 k))
       (lst fs (cdr lst))
       (L 0 (+ (* -1 (expt -1 k) (choose n k) (car lst)) L)))
      ((null? lst) L)))

(define (sequence->limit proc sequence)
  (define lval (proc (car sequence)))
  (if (finite? lval)
      (let ((val (proc (cadr sequence))))
	(define h_n*nsamps (* (length sequence) (magnitude (- val lval))))
	(if (finite? val)
	    (let loop ((sequence (cddr sequence))
		       (fxs (list val lval))
		       (trend #f)
		       (ldelta (- val lval))
		       (jdx (+ -1 (length sequence))))
	      (cond ((null? sequence)
		     (case trend
		       ((diverging) (and (real? val) (/ ldelta 0.0)))
		       ((bounded) (invintp val lval (caddr fxs)))
		       (else (cond ((zero? ldelta) val)
				   ((not (real? val)) #f)
				   (else (extrapolate-0 fxs))))))
		     (set! lval val)
		     (set! val (proc (car sequence)))
		     (if (finite? val)
			 (let ((delta (- val lval)))
			   (define h_j (/ h_n*nsamps jdx))
			   (cond ((case trend
				    ((converging) (<= (magnitude delta) h_j))
				    ((bounded)    (<= (magnitude ldelta) (magnitude delta)))
				    ((diverging)  (>= (magnitude delta) h_j))
				    (else #f))
				  (loop (cdr sequence) (cons val fxs) trend delta (+ -1 jdx)))
				 (trend #f)
				  (loop (cdr sequence) (cons val fxs)
					(cond ((> (magnitude delta) h_j) 'diverging)
					      ((< (magnitude ldelta) (magnitude delta)) 'bounded)
					      (else 'converging))
					delta (+ -1 jdx)))))
			 (and (eq? trend 'diverging) val)))))
	    (and (real? val) val)))
      (and (real? lval) lval)))

(define (limit proc x1 x2 . k)
  (set! k (if (null? k) 8 (car k)))
  (cond ((not (finite? x2)) (slib:error 'limit 'infinite 'x2 x2))
	((not (finite? x1))
	 (or (positive? (* x1 x2)) (slib:error 'limit 'start 'mismatch x1 x2))
	 (limit (lambda (x) (proc (/ x))) 0.0 (/ x2) k))
	((= x1 (+ x1 x2)) (slib:error 'limit 'null 'range x1 (+ x1 x2)))
	(else (let ((dec (/ x2 k)))
		(do ((x (+ x1 x2 0.0) (- x dec))
		     (cnt (+ -1 k) (+ -1 cnt))
		     (lst '() (cons x lst)))
		    ((negative? cnt)
		     (sequence->limit proc (reverse lst))))))))


Copyright (C) Aubrey Jaffer 2005. All Rights Reserved.

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.


Editor: Mike Sperber
Last modified: Mon Aug 8 07:42:55 CEST 2005