SICP Exercise 1.19: Fibonacci on Steroids

There is a clever algorithm for computing the Fibonacci numbers in a logarithmic number of steps. Recall the transformation of the state variables a and b in the fib-iter process of section 1.2.2: a ← a + b and b ← a. Call this transformation T, and observe that applying T over and over again n times, starting with 1 and 0, produces the pair Fib(n + 1) and Fib(n). In other words, the Fibonacci numbers are produced by applying T_n, the n^th power of the transformation T, starting with the pair (1,0). Now consider T to be the special case of p = 0 and q = 1 in a family of transformations T_pq, where T_pq transforms the pair (a,b) according to a ← bq + aq + ap and b ← bp + aq. Show that if we apply such a transformation T_pq twice, the effect is the same as using a single transformation T_p'q' of the same form, and compute p' and q' in terms of p and q. This gives us an explicit way to square these transformations, and thus we can compute T_n using successive squaring, as in the fast-expt procedure. Put this all together to complete the following procedure, which runs in a logarithmic number of steps:

(define (fib n)
  (fib-iter 1 0 0 1 n))
(define (fib-iter a b p q count)
  (cond ((= count 0) b)
        ((even? count)
         (fib-iter a
                   b
                         ; compute p'
                         ; compute q'
                   (/ count 2)))
        (else (fib-iter (+ (* b q) (* a q) (* a p))
                        (+ (* b p) (* a q))
                        p
                        q
                        (- count 1)))))

We're given the transformation:

Tpq(a, b) = ((bq + aq + ap), (bp + aq))

If we apply the transformation T_pq to itself then we get the following:

  Tpq(Tpq(a, b))
= Tpq((bq + aq + ap), (bp + aq))
= (((bp + aq)q + (bq + aq + ap)q + (bq + aq + ap)p), ((bp + aq)p + (bq + aq + ap)q))
= ((bpq + aq2  + bq2 + aq2 + apq + bpq + apq + ap2), (bp2 + apq + bq2 + aq2 + apq))
= ((b(q2 + 2pq) + a(q2 + 2pq) + a(p2 + q2)), (b(p2 + q2) + a(q2 + 2pq)))

Now if we define p' and q' as:

p' = p2 + q2
q' = q2 + 2pq

...then this gives:

Tpq(Tpq(a, b)) = ((bq' + aq' + ap'), (bp' + aq'))

...which is of the same form as the transform T_pq, and so we can define T_p'q'(a, b) as:

  Tp'q'(a, b)
= ((bq' + aq' + ap'), (bp' + aq'))
= ((b(q2 + 2pq) + a(q2 + 2pq) + a(p2 + q2)), (b(p2 + q2) + a(q2 + 2pq)))

We can now simply plug the mappings for p' and q' into the template procedure above to give:

(define (fib n) (fib-iter 1 0 0 1 n))
(define (fib-iter a b p q count)
  (cond ((= count 0) b)
        ((even? count)
         (fib-iter a
                   b
                   (+ (square p) (square q))
                   (+ (* 2 p q) (square q))
                   (/ count 2)))
        (else (fib-iter (+ (* b q) (* a q) (* a p))
                        (+ (* b p) (* a q))
                        p
                        q
                        (- count 1)))))

Finally, we can check it works as expected:

> (fib 0)
0
> (fib 1)
1
> (fib 2)
1
> (fib 3)
2
> (fib 4)
3
> (fib 5)
5
> (fib 6)
8
> (fib 7)
13
> (fib 8)
21