;;; -*- Mode:LISP; Package:USER; Base:10; Readtable:CL; Fonts:(CPTFONT CPTFONTB) -*-
1;;;FIXED-POINT.LISP
;;;2/23/88 - Keith Corbett
;;;The purpose of this example is to illustrate a few problems
;;;that LISP* 1programmers frequently run into with floating-point
;;;calculations.* 1Actually, these are concerns in any programming
;;;language, but the* 1variety of datatypes and numeric functions
;;;in LISP makes these* 1problems even more critical.
;;;
;;;The most important feature of this example is the trick of
;;;using an* 1EPSILON-value in calculations involving
;;;floating-point numbers. Any* 1repeated calculation that involve
;;;convergence upon a value -- even* 1simple comparisons, such as
;;;(ZEROP ) -- may yield poor results* 1unless an appropriate
;;;EPSILON-value is used.
;;;
;;;It is important to note that this number, which embodies the
;;;concept* 1"as close as we're going to get to the desired
;;;result", is the* 1same as "the smallest number we can
;;;express with a floating-point* 1number. For example,
;;;Common-LISP defines a constant value,* 1SINGLE-FLOAT-EPSILON.
;;;This value is the absolute value of the* 1smallest difference
;;;that can be detected between two numbers. In* 1other words, if
;;;the absolute value of the difference between two* 1floats is
;;;more than SINGLE-FLOAT-EPSILON, then the two numbers are not
;;;(floating-point) equivalent.
;;;
;;;We can demonstrate that the two kinds of EPSILON-values are
;;;not the* 1same by looking at an example of a numerical algorithm
;;;that uses a* 1convergence technique. The example code below,
;;;with the comments* 1preceding it, is from BYTE magazine,
;;;February 1988, from the* 1article "LISP: A Language for
;;;Stratified Design" by Abelson and* 1Sussman. When I implemented
;;;the example, I had to make two* 1changes: 1) I translated from
;;;the Scheme dialect to Common-LISP, and* 12) I had to use a
;;;greater EPSILON-value than they did. (They used* 11.0e-10, I
;;;used 1.0e-9). Attempting to use an EPSILON of lower* 1magnitude
;;;on the Lambda caused my code to loop, swapping* 1alternately
;;;between two* 1values. This is a typical result* 1that can* 1arise
;;;from a* 1failure to determine, and use, EPSILON-values* 1properly.
;;;
;;;Note* 1that I* 1had to* 1determine a* 1workable* 1EPSILON-value by
;;;experimentation; standard mathematics texts* 1on this problem
;;;would probably show a better technique for* 1calculating the
;;;appropriate value in a given circumstance.
;;;
;;;Abelson & Sussman illustrate a second problem that* 1can cause a
;;;numerical algorithm to loop; see AVERAGE-DAMP, below.
;;;
;;;The third important feature of this example is to illustrate
;;;the pitfalls of integral and rational arithmetic where
;;;floating-point* 1calculations are desired. In particular, using
;;;the* 1Common-LISP function `/' with integral arguments can yield
;;;a* 1strange-looking* 1ratio,* 1rather than the nice floating-point
;;;result you probably* 1want.* 1See my comments at the end.
;;;
;;; -KmC*
;;1Abelson & Sussman: The "Boring Meeting Algorithm":*
;;
;;A clear way to formulate {certain numeric algorithms, such as SQRT,}
;;is as a process of computing a FIXED POINT. {For example,} the square
;;root of a radicand X is the number Y such that Y = X / Y, or, in other
;;words, Y is a fixed point of the procedure (LAMBDA(Y) (/ X Y)).
;;
;;How do you find a fixed point of a function? In favorable cases, you
;;can iterate the function until the result is close to the input. For
;;example, during boring meetings many of us have noticed that we can
;;find the fixed point of the COSINE function by entering 1 on a pocket
;;calculator and repeatedly pressing the COSINE button. After a while,
;;the calculated value converges to approximately .739085. You can
;;capture that general idea as the FIXED-POINT procedure shown in listing
;;2 {below}. Procedure FIXED-POINT, given a one-argument procedure F and
;;an INITIAL-VALUE, keeps applying F until successive values are
;;sufficiently close to each other.
;;
;;Listing 2:
(defun fixed-point(f initial-value
&optional (epsilon 1.0e-9))
(labels((close-enough? (v1 v2)
(< (abs (- v1 v2)) epsilon))
(loopnext (value)
(let((next-value (funcall f value)))
(if (close-enough? value next-value)
next-value
(loopnext next-value)))))
(loopnext initial-value)))
1;;;Note: On the Lambda, (fixed-point #'cos 1) --> .7390851313*,
1;;;as described in the article.*
;;1Abelson & Sussman, continued:*
;;
;;You can attempt to find square roots by using the following definition
;;of SQRT: (DEFINE (SQRT X) (FIXED-POINT (LAMBDA(Y) (/ X Y)) 1))
;;
;;Unfortunately, this doesn't work. Unlike the COSINE function, applying
;;the indicated procedure over and over does not converge to a fixed
;;point, but rather alternates between the same two values, which are on
;;opposite sides of the square root.
;;
;;In situations like this, you can often force convergence by averaging.
;;The AVERAGE-DAMP procedure shown in listing 3 {below} takes as its
;;argument a procedure that computes a function F and returns as its
;;result a procedure that computes a function with the same fixed points
;;as F, but whose oscillations are damped out by averaging successive
;;values. The new definition of SQRT shows how to use AVERAGE-DAMP to
;;express the square root method as a process of finding the fixed point
;;of an AVERAGE-DAMPed function.
;;
;;Listing 3:
(defun average-damp(f)
(lambda(x)
(quotient (+ x (funcall f x)) 2)))
(defun sqrt2(x)
(fixed-point (average-damp (lambda(y) (/ x y))) 1.))
1;;;KmC: Note that for SQRT2, above, I used the Common-LISP #'/
;;;function to* 1divide X by Y. What happens if this function is
;;;passed two* 1integer-valued arguments?
;;;
;;; (SQRT2 9) -->
;;;340282366920938463463374607431768211457/113427455640312821154458202477256070485
;;;
;;;In other words, the #`/ function returns a RATIO which is
;;;numerically* 1equivalent to the `correct answer', 3.0. Confirm
;;;this for yourself by* 1evaluating (FLOAT(SQRT2 9)).
;;;
;;;You can get a better looking result from SQRT2 by passing it* 1a
;;;floating-point argument, i.e. (SQRT2 9.0).
;;;
;;;This problem is not a bug, but it is something you have to
;;;consider* 1in your programs. There are two things to do to
;;;correct this example: 1]* 1set X to the FLOAT of itself at the
;;;beginning of SQRT2, and 2] use the* 1#'QUOTIENT function instead
;;;of #'/ [which may be overkill, but it* 1is more in the* 1spirit of
;;;what we want].
;;;
;;;For example:*
(defun sqrt3(x)
(setq x (float x))
(fixed-point (average-damp (lambda(y) (quotient x y))) 1.))
1;;;Now, (SQRT3 9) --> 3.0, which is what we wanted.