Saturday, August 21, 2010

SICP 1.45: Computing nth roots

From SICP section 1.3.4 Procedures as Returned Values

Exercise 1.45 asks us to do some experiments to find out how many average damps are needed to compute nth roots as a fixed-point search based on repeated average damping of y → x / yn-1. Once we know that, we need to write a general procedure for finding nth roots using fixed-point, average-damp, and the repeated procedure from exercise 1.43.

We already saw in section 1.3.3 that computing a square root by finding a fixed point of y → y/x doesn't converge, and that this can be fixed by average damping (see exercise 1.36). The same method (a single average damp) works for finding cube roots as fixed points of y → x / y2.

Here's a first attempt at an nth-root procedure, which works when n is 2 or 3.
(define (nth-root x n)
(fixed-point
(average-damp
(lambda (y) (/ x (expt y (- n 1)))))
1.0))

> (nth-root 100 2)
10.0
> (nth-root 1000 3)
10.000002544054729

The same procedure does not work for fourth roots. Try running nth-root with the following parameters and you'll see that the procedure never finishes.
> (nth-root 10000 4)

This is because the fixed-point search for y → x / y3 doesn't converge when only one average damp is used. We can fix this by damping the average twice.
(define (nth-root x n)
(fixed-point
((repeated average-damp 2)
(lambda (y) (/ x (expt y (- n 1)))))
1.0))

We'll need to run the new procedure several times to find out where average damping twice fails to converge. The code above works for finding nth roots when n is less than 8.
> (nth-root 10000 4)
10.0
> (nth-root 100000 5)
9.99999869212542
> (nth-root 1000000 6)
9.999996858149522
> (nth-root 10000000 7)
9.9999964240619

Let's increase the repetitions of average-damp one more time and see if we can spot a pattern.
(define (nth-root x n)
(fixed-point
((repeated average-damp 3)
(lambda (y) (/ x (expt y (- n 1)))))
1.0))

> (nth-root 256 8)
2.0000000000039666
> (nth-root 512 9)
1.9999997106840102
> (nth-root 1024 10)
2.000001183010332
> (nth-root 2048 11)
1.999997600654736
> (nth-root 4096 12)
1.9999976914703093
> (nth-root 8192 13)
2.0000029085658984
> (nth-root 16384 14)
1.9999963265447058
> (nth-root 32768 15)
2.0000040951543023

When we average damp 3 times, the nth-root procedure converges for n up to 15. That worked longer than I expected, but it gives us enough information to see the pattern. When n (roughly) doubles, we need to increase the number of average damps by one.

maximum n: 3, 7, 15
average damps: 1, 2, 3

Where a is the number of average damps, we can say

nmax = 2a+1 - 1

We can test this by checking to see that by increasing the number of average damps to four, our nth-root procedure works for values of n up to 31.
(define (nth-root x n)
(fixed-point
((repeated average-damp 4)
(lambda (y) (/ x (expt y (- n 1)))))
1.0))

> (nth-root 2147483648 31)
1.9999951809750391

Just as expected, the procedure fails to converge if we run it with n = 32. I'm convinced that this is the right pattern.

In order to calculate the number of average damps from n, we just need to take the log2 of n then floor the result. Since Scheme only has a log procedure that finds logarithms in base e, this must be what the authors meant when they said in the exercise "Assume that any arithmetic operations you need are available as primitives."

It's easy enough to write a procedure to compute log2 if we just remember that

logn(x) = log(x) / log(n)

(define (log2 x)
(/ (log x) (log 2)))

> (log2 16)
4.0
> (log2 64)
6.0
> (log2 65536)
16.0

With that we can fix the nth-root procedure so it works for any value of n.
(define (nth-root x n)
(fixed-point
((repeated average-damp (floor (log2 n)))
(lambda (y) (/ x (expt y (- n 1)))))
1.0))

> (nth-root 4294967296 32)
2.000000000000006
> (nth-root 18446744073709551616 64)
2.0000000000000853
> (nth-root 340282366920938463463374607431768211456 128)
2.0000000000082006


Related:

For links to all of the SICP lecture notes and exercises that I've done so far, see The SICP Challenge.

2 comments:

Anonymous said...

Bill,

Did you try using 3 average damps beyond 15? 3 average damps works all the way up to 31-roots. For 32-roots and beyond, 4 average damps are required.

In fact, 4 average damps seems to work fine for 64-roots, 128-roots, 256-roots, and so on.

Actually, 3 average damps works for 512-roots too. Perhaps I skipped over the n-root where 4 average damps first fails.

Anonymous said...

It seems that if the index any integer
2^n. it works with less damps. Nice thing with the 'floor' procedure, I didn't know about it, and wrote a little procedure to
give me the integer lesser than the log2.