Square roots by Newton’s method

Question

The following Scheme program implements Newton’s method for computing the square root of a number:

(import (scheme small))

(define (sqrt x)
  (define (sqrt-iter guess)
    (if (good-enough? guess)
      guess
      (sqrt-iter (improve guess))))
  (define (good-enough? guess)
    (define tolerance 0.001)
    (< (abs (- (square guess) x)) tolerance))
  (define (improve guess)
    (if (= guess 0) guess (average guess (/ x guess))))
  (define (average x y)
    (/ (+ x y) 2))
  (define initial-guess 1.0)
  (sqrt-iter initial-guess))

(display (sqrt 0)) (newline)
(display (sqrt 1e-12)) (newline)
(display (sqrt 1e-10)) (newline)
(display (sqrt 1e-8)) (newline)
(display (sqrt 1e-6)) (newline)
(display (sqrt 1e-4)) (newline)
(display (sqrt 1e-2)) (newline)
(display (sqrt 1e0)) (newline)
(display (sqrt 1e2)) (newline)
(display (sqrt 1e4)) (newline)
(display (sqrt 1e6)) (newline)
(display (sqrt 1e8)) (newline)
(display (sqrt 1e10)) (newline)
(display (sqrt 1e12)) (newline)
(display (sqrt 1e13))

Output:

0.03125
0.031250000010656254
0.03125000106562499
0.03125010656242753
0.031260655525445276
0.03230844833048122
0.10032578510960605
1.0
10.000000000139897
100.00000025490743
1000.0000000000118
10000.0
100000.0
1000000.0
[Hanging forever…]

As we can see, this naive program does not perform well:

for small numbers (below x = 1e-2), the tolerance 0.001 is too large;
for large numbers (from x = 1e13), the program hangs forever.

Both problems can be solved by redefining the good-enough? procedure like this:

(define (good-enough? guess)
  (= (improve guess) guess))

But this solution is not the purpose of my post. Instead, I would like to understand why the naive program fails the way it fails for large numbers.

I have not read the IEEE Standard for Floating-Point Arithmetic (IEEE 754), but to my understanding, floating-point numbers cannot represent all reals and have an absolute precision that is very high for small numbers and very low for large numbers (this Wikipedia figure seems to confirm this). In other words, small floating-point numbers are dense and large floating-point numbers are sparse. The consequence of this, is that the naive program will hang forever, trapped in an infinite recursion, if guess has not reached the tolerance range yet and if (improve guess) cannot improve the guess anymore because the distance between the new guess and the old guess is below the absolute precision of the old guess (so the new guess is the same floating-point number as the old guess, meaning that a fixed point of (improve guess) has been reached).

To guarantee that, for a given x, the naive program returns, it seems to me that this predicate must hold:

tolerance > absolute_precision(sqrt(x)).

If we choose a tolerance of 0.001 = 1e-3, that means that the absolute precision of sqrt(x) should be less than 1e-3. Consequently, according to the Wikipedia figure above for binary64 floating-point numbers, sqrt(x) should be less than 1e13 and therefore x should be less than 1e26.

Here are my questions:

Is my predicate correct?
Why does the program hang forever from x = 1e13 instead of the expected x = 1e26?
Why does the program hang forever with x in {1e13, 1e15, 1e17, …, 1e49} and any x greater than 1e49 but still return with x in {1e14, 1e16, 1e18, …, 1e48}?

Note. — I am using the Chibi Scheme interpreter on a 64-bit MacBook Pro so the IEEE 754 binary64 format.

Sylwester Sylwester · Accepted Answer · 2020-01-25T20:54:41

When you square 3162277.6601683795 you get 10000000000000.002 (in binary this is infinite so looses precision compared to the argument). If you convert 10000000000000002 to hex you get 2386F26FC10002 which is 4x13 nibbles and 2 uses 2 bits. In float the msb is always 1 and is omitted so I'm guessing 53 bits is used and lots of decimals are skipped since that is how floating point works. You cannot get closer to 10000000000000 by tolerance 0.001 from either side

I'm thinking you could compare the guess with the previous guess and see if that is below a tolerance. eg.

(define (sqrt x)
  (define (sqrt-iter guess prev-guess)
    (if (good-enough? guess prev-guess)
      guess
      (sqrt-iter (improve guess) guess)))
  (define (good-enough? guess prev-guess)
    (define tolerance 1e-20)
    (< (abs (- guess prev-guess)) tolerance))
  (define (improve guess)
    (if (= guess 0) guess (average guess (/ x guess))))
  (define (average x y)
    (/ (+ x y) 2))
  (define initial-guess 1.0)
  (sqrt-iter initial-guess 0))

This results in 3162277.6601683795. Unlike your version when you substract the guesses missing bits will make the difference smaller so even with a tolerance of 1e-20 it will produce an answer for 1e270.

Now if you look at the guesses it makes it will jump quite large making the new guess always much better so only when the guess is close to the actual result it will be close to the previous guess and I'm guessing you would get as good estimate if your tolerance were just squared from the previous.

I agree that the fixed tolerance is no good and we should have had something that looks at the guesses divided with each other instead. When they are close the result will always be close to 1 and while it is off by much it will be close to 2 or 1/2 depending. Just comparing without tolerance might put you in a an infinite loop when the improve will jump over and under the mark and never match so the tolerance was put in there for a reason.

Square roots by Newton’s method

1 Answers