; 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