Update: I changed the argument 1 in nth-root to 1.0, and then it works.
Why are they different?
Here's my original question.
I have
(define tolerance 0.00001)
(define (fixed-point f first-guess)
(define (close-enough? v1 v2)
(< (abs (- v1 v2)) tolerance))
(define (try guess)
(let ((next (f guess)))
(if (close-enough? guess next)
next
(try next))))
(try first-guess))
(define (average v1 v2)
(/ (+ v1 v2) 2))
(define (average-damp f)
(lambda (x)
(average x (f x))))
(define (compose f g)
(lambda (x)
(f (g x))))
(define (repeated f n)
(cond ((= n 1)
f)
(else
(compose
f
(repeated f (- n 1))))))
(define (log2 x)
(/ (log x) (log 2)))
(define (nth-root n x)
(fixed-point
((repeated average-damp (floor (log2 n)))
(lambda (y) (/ x (expt y (- n 1)))))
1))
All the procedures except nth-root work well in previous Exercises
When I enter (nth-root 2 9), it outputs 3.00...
But when I enter (nth-root 3 8), it outputs nothing. Why the recursion can not be stopped?
I can't find where is wrong
If you look at what happens when you try to compute (nth-root 3 9) you'll discover something interesting: the program is not in fact looping endlessly and failing to converge: it is looping only a few times with each iteration taking a dramatically increasing time.
Here is the log of the first few iteration times (specifically the time between eah call of try in fixed-point) on my machine, together with about how many times slower each iteration is than the last.
> (nth-root 3 8)
0.000008
0.000141 (factor of ~17.6)
0.000055 (factor of ~0.4)
0.000057 (factor of ~1)
0.000072 (factor of ~1.3)
0.000225 (factor of ~3.1)
0.001119 (factor of ~5)
0.008606 (factor of ~7.7)
0.059756 (factor of ~6.9)
0.312313 (factor of ~5.2)
2.578403 (factor of ~8.3)
22.959139 (factor of ~8.9)
204.239457 (factor of ~8.9)
1811.492695 (factor of ~8.9)
The last iteration before I killed it took more than half an hour.
So it is not failing to converge at all: rather it is (probably) converging but it is taking a very, very long time to do so: each iteration is taking something like 9 times as the previous one meaning the process is exponential in the number of iterations. So I believe that this will, eventually, yield an answer. The problem is that 'eventually' is 'after a very, very long time'.
The reason for this is that computing (nth-root 3 8) means you are doing exact arithmetic, and in particular you're doing arithmetic on exact rationals, and these rationals rapidly grow to have catastrophically large numerators and denominators, so the rational arithmetic functions become extremely slow.
In particular the assumption that arithmetic operations on numbers represented exactly is constant time is false.
I tested this in Racket (which is already not a very quick implementation) by writing wrappers for some arithmetic operators which print the lengths of the numerators & denominators of their arguments in bits (it is much better to do this than to try and print the numbers themselves!).
Here is the definition of the wrappers (this is Racket-specific code):
(define (ndl x)
(list (integer-length (numerator x))
(integer-length (denominator x))))
(define-syntax-rule (define-wrapped name op)
(define (name . args)
(cond
[(andmap exact? args)
(printf "~A ~S~%" 'op (map ndl args))
(define r (apply op args))
(printf " -> ~S~%" (ndl r))
r]
[else
(apply op args)])))
(define-wrapped plus +)
(define-wrapped minus -)
(define-wrapped divide /)
(define-wrapped exponent expt)
So, then, if you replace the appropriate functions in the program by these things and then run it, you'll get a bunch of output which starts to look like this:
expt ((65735 65734) (2 1))
-> (131470 131468)
/ ((4 1) (131470 131468))
-> (131471 131470)
+ ((65735 65734) (131471 131470))
-> (197206 197204)
/ ((197206 197204) (2 1))
-> (197206 197205)
- ((65735 65734) (197206 197205))
-> (197195 197205)
- ((2 1) (1 1))
-> (2 1)
expt ((197206 197205) (2 1))
-> (394411 394409)
/ ((4 1) (394411 394409))
-> (394412 394411)
+ ((197206 197205) (394412 394411))
-> (591618 591616)
/ ((591618 591616) (2 1))
-> (591618 591617)
- ((197206 197205) (591618 591617))
-> (591606 591617)
- ((2 1) (1 1))
-> (2 1)
expt ((591618 591617) (2 1))
-> (1183235 1183233)
/ ((4 1) (1183235 1183233))
-> (1183236 1183235)
+ ((591618 591617) (1183236 1183235))
-> (1774853 1774851)
/ ((1774853 1774851) (2 1))
The last number has about 1.7 million bits in both its numerator and denominator, or over half a million decimal digits. And you can also see that the number of bits in these numbers is growing approximately exponentially for each iteration (and in fact probably actually exponentially). In other words the log of the magnitude of these things is growing exponentially. The time taken for an operation on any of these numbers is at best linear in the number of bits, so this means the whole process is, at best, exponential (assuming the growth in numbers of bits is indeed exponential). It will also be exponential in space of course, but that's probably not hurting by the time the exponential time kills you.
Eventually you'll get an answer, I think.
In response to Will Ness's comment, here is a Racket version of the above code which uses its bigfloat library to allow you to compute things to high, but controlled, precision. To use this in a meaningful way you need to think harder than I was bothered to about the interaction berween (bf-precision), which controls how many bits of precision the bigfloats have, tolerance as defined by the program and the magnitude of the result.
(require math/bigfloat)
(define tolerance (bf "0.00001"))
(define (fixed-point f first-guess)
(define (close-enough? v1 v2)
(bf< (bfabs (bf- v1 v2)) tolerance))
(define (try guess)
(let ((next (f guess)))
(if (close-enough? guess next)
next
(try next))))
(try first-guess))
(define (average v1 v2)
(bf/ (bf+ v1 v2) 2.bf))
(define (average-damp f)
(lambda (x)
(average x (f x))))
(define (compose f g)
(lambda (x)
(f (g x))))
(define (repeated f n)
(cond ((= n 1)
f)
(else
(compose
f
(repeated f (- n 1))))))
(define (log2 x)
(/ (log x) (log 2)))
(define (nth-root n x)
(fixed-point
((repeated average-damp (floor (log2 n)))
(lambda (y) (bf/ (bf x) (bfexpt y (bf (- n 1))))))
1.bf))