; Newton's method square root termination conditions for large and small numbers
(define (square x) (* x x))
; Note: average overflows for large integers
(define (average x y) (/ (+ x y) 2))
(define (make-sqrt good-enough?)
(define (sqrt x)
(define (improve guess) (average guess (/ x guess)))
(define (sqrt-iter guess)
(if (good-enough? guess x) guess
(sqrt-iter (improve guess))))
(sqrt-iter 1.0))
sqrt)
(define (good-enough-absolute? guess x)
(define epsilon 0.001)
(< (abs (- (square guess) x)) epsilon))
(define sqrt-absolute (make-sqrt good-enough-absolute?))
(define (pow a n)
(cond ((> n 0) (* (pow a (- n 1)) a))
((< n 0) (/ (pow a (+ n 1)) a))
(else 1.0)))
; Generate list with the sequence of integers a, a + 1, a + 2, ..., b - 1
(define (range a b)
(if (< a b) (cons a (range (+ a 1) b)) '()))
(define (print . items) (map display items))
(define (println . items)
(print items)
(newline))
(define (evaluate-sqrt target-root sqrt)
(define target-square (square target-root))
(println target-square " = " target-root "^2")
(define root (sqrt target-square))
(println "sqrt(" target-square ") = " root)
(println "absolute difference: " (abs (- target-root root)))
(println "percent difference: " (/ root target-root)))
; For small numbers the percent error can be huge even if the absolute error is
; tiny. This means |guess^2 - x| < epsilon is an ineffective selection criteria
; for small numbers.
(println "Absolute error threshold for small numbers")
(evaluate-sqrt (pow 10.0 -20) sqrt-absolute)
; (1.0000000000000001e-40 = 1.0000000000000001e-20 ^2)
; (sqrt( 1.0000000000000001e-40 ) = .03125)
; (absolute difference: .03125)
; (percent difference: 3.1249999999999995e18)
; For large numbers the absolute error can be large even when the percent error
; is small. This means the computation will require increasingly large numbers
; of iterations to converge to a square root for large numbers.
; Returns a function which counts how many iterations are required to compute
; the square root of its input using good-enough?
(define (make-iteration-count-sqrt good-enough?)
(define (sqrt-iteration-count x)
(define (improve guess) (average guess (/ x guess)))
(define (sqrt-iter guess num-iterations)
(if (good-enough? guess x) num-iterations
(sqrt-iter (improve guess) (+ num-iterations 1))))
(sqrt-iter 1.0 1))
sqrt-iteration-count)
(define sqrt-absolute-iteration-count
(make-iteration-count-sqrt good-enough-absolute?))
(define (evaluate-count x sqrt-iteration-count)
(println "sqrt(" x ") takes " (sqrt-iteration-count x) " iterations"))
(println "Iteration counts for |guess^2 - x| < epsilon termination condition")
(map (lambda (x) (evaluate-count x sqrt-absolute-iteration-count))
(cons (pow 10 100) '(0.00000001 0.001 10 100 1000 1000000 1000000000)))
; Percent error |1 - guess^2/x| < epsilon
; Note: this breaks for 0
(define (good-enough-relative? guess x)
(define epsilon 0.001)
(< (abs (- (/ (square guess) x) 1.0)) epsilon))
(define sqrt-relative (make-sqrt good-enough-relative?))
(println "Percent error threshold for small and large numbers")
(evaluate-sqrt (pow 10.0 -20) sqrt-relative)
(evaluate-sqrt (pow 10.0 20) sqrt-relative)
(define sqrt-relative-iteration-count
(make-iteration-count-sqrt good-enough-relative?))
(println "Iteration counts for percent error termination condition")
(map (lambda (x) (evaluate-count x sqrt-relative-iteration-count))
(cons (pow 10 100) '(0.00000001 0.001 10 100 1000 1000000 1000000000)))
; Using percent error solves the problem of inaccuracies for small numbers but
; does not reduce the number of iterations for large numbers by very much.
; Another problem with with using absolute difference is for large numbers the
; improve step will use the previous guess as the next guess causing the
; computation to infinite loop.
; For large floating point numbers, the distances between numbers increases.
; When gn ~ sqrt(x)
; then (gn + x/gn)/2 = (gn + ~gn)/2 can equal gn because there is no floating
; point number between gn and x/gn.
; I could not produce this behavior with mit-scheme but I could with biwascheme
; (sqrt-absolute 1000000000000000)
;
; The computation start to infinite loop at guess = 31622776.601683795
>>> guess, x = 31622776.601683795, 1000000000000000
>>> guess, x/guess, guess - x/guess
(31622776.601683795, 31622776.60168379, 3.725290298461914e-09)
>>> next_guess = (guess + x/guess)/2
>>> next_guess - guess, next_guess == guess
(0.0, True)
>>> x - guess * guess
-0.125
; TODO(adrs): implement the method the book suggests