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:
Post a Comment