;;; -*- 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 <x>) -- 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 <NOT> 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 ide 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.