Wednesday, August 27, 2014

A use of Newton's method

I've seen more than one book claim that computing with real numbers inevitably involves round-off errors because real numbers can have an infinite number of digits after the decimal point and no finite representation can hold them. This is false. Instead of representing a real number as a nearby rational number with an error introduced by rounding, we'll represent a real number as computer program that generates the digits. The number of digits generated is potentially infinite, but the program that generates them is definitely finite.

Here is Gosper's algorithm for computing the square root of a rational number.
(define (gosper-sqrt a b c d)
  ;;   Solve for
  ;; ax + b 
  ;; ------  = x
  ;; cx + d
  (define (newtons-method f f-prime guess)
    (let ((dy (f guess)))
      (if (< (abs dy) 1)
          guess
          (let ((dy/dx (f-prime guess)))
            (newtons-method f f-prime (- guess (/ dy dy/dx)))))))

  (define (f x)
    (+ (* c x x)
       (* (- d a) x)
       (- b)))

  (define (f-prime x)  
    (+ (* 2 c x)
       (- d a)))

  (let ((value (floor (newtons-method f f-prime b))))
    (cons-stream value
                 (gosper-sqrt (+ (* c value) d)
                              c
                              (+ (* (- a (* value c)) value) 
                                 (- b (* value d)))
                              (- a (* value c))))))

1 ]=> (cf:render (gosper-sqrt 0 17 10 0))
1.303840481040529742916594311485836883305618755782013091790079369...

;; base 10, 100 digits
1 ]=> (cf:render (gosper-sqrt 0 17 10 0) 10 100)
1.303840481040529742916594311485836883305618755782013091790079369896765385576397896545183528886788497...

No comments: