| From: Sebastian Egner <xxxxxx@philips.com>
| Date: Wed, 25 May 2005 09:35:25 +0200
|
| Example: An important function from information theory is
|
| f(x) = -x log(x).
|
| This function is in principle well behaved (smooth, analytic, etc.)
| on (0,1], but its derivative does not exist at x = 0. Moreover,
| f(0) cannot directly be computed numerically because the underflow
| from log(x) is not cancelled by the multiplication with zero.
| Practical numerical code: IF x < xmin THEN 0 ELSE x log(x), where
| xmin is chosen minimal such that log(xmin) is an ordinary number
| and not -infinity.
|
| Using LIMIT in this case is not a good idea for two reasons:
| a) It is expensive and unnecessary, except for very small x.
| b) At least the reference implementation of LIMIT doesn't get it right:
|
| (limit (lambda (x) (* -1 x (log x))) 0 1e-9) => -inf.0
|
| This may be a bug in the reference implementation, but it is certainly
| a violation of the new specification as f(x) is monotonic on [0,1/e].
That bug has been corrected:
> (limit (lambda (x) (* -1 x (log x))) 0 1e-9)
102.70578127633066e-12
> (limit (lambda (x) (* -1 x (log x))) 0 1e-222)
0.0
http://savannah.gnu.org/cgi-bin/viewcvs/scm/scm/Transcen.scm?rev=HEAD&content-type=text/vnd.viewcvs-markup
has the new limit code.
http://swiss.csail.mit.edu/~jaffer/III/Limit.html
explains the mathematics and algorithm.
The documentation for the new limit procedure is appended.
| When you try to fix the reference implementation, you will find that
| it cannot be fixed because it comes from the "black box" procedure:
It took some hard work, but the new limit procedure does treat its
procedure argument as a black box, never calling it at the limit
point. It essentially detects whether the function is convex or
concave near the limit point and uses polynomial extrapolation, either
from the function or its inverse, to estimate the limit value. The
inverse function extrapolation is currently limited to quadratic;
higher order fitting can certainly be done, but the equations get big.
| At a certain moment f(x) becomes -inf.0, so that must be the limit.
The new limit procedure calls f at points which are controlled by
arguments to limit: X1+X2/K, X1+2*X2/K, ..., X1+X2. The stipulation
that f be monotonic and continuously differentiable over the range
(X1, X1+X2] makes polynomial extrapolation a reasonable approximation.
| Bottom line: In the end, LIMIT might do more harm than it is worth.
| You might want to reconsider if it is a feature that is essential for
| the Scheme programming language itself.
The new limit procedure works better than I had expected; even
handling a nasty example from an old Macsyma manual (see appended
examples). Arguments that limits are too esoteric or new
counterexamples may yet scuttle them from srfi-70.
-=-=-=-=-
Function: limit proc x1 x2 k
Function: 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) ==> 1/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 + 1/0 1.0e9) ==> 1/0
(limit + -1/0 -1.0e9) ==> -1/0
(limit / 0 1.0e-9) ==> 1/0
(limit / 0 -1.0e-9) ==> -1/0
(limit tan (atan 1/0) -1.0e-15) ==> 1/0
(limit tan (atan -1/0) 1.0e-15) ==> -1/0
(limit (lambda (x) (/ (log x) x)) 0 1.0e-9) ==> -1/0
(limit (lambda (x) (/ (magnitude (log x)) x)) 0 -1.0e-9)
==> -1/0
;; limit doesn't exist
(limit sin 1/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)) 1/0 1.0e222) ==> #f
;; asymptotes
(limit / -1/0 -1.0e222) ==> 0.0
(limit / 1/0 1.0e222) ==> 0.0
(limit (lambda (x) (expt x x)) 0 1.0e-18) ==> 1.0
(limit (lambda (x) (sin (/ x))) 1/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))) 1/0 40)
==> 7.3890560989306504
