SICP Exercise 1.19: Computing Fibonacci numbers

From SICP section 1.2.4 Exponentiation

Exercise 1.19 asks us to complete a procedure for computing Fibonacci numbers in a logarithmic number of steps. The following code is given:

(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 reminded of the transformation of the state variables a and b in the original fib-iter procedure from section 1.2.2: a ← a + b and b ← a. If these state changes are labeled transformation T, then applying T repeatedly for n iterations starting with a = 1 and b = 0 produces the pair a = Fib(n + 1) and b = Fib(n). So the Fibonacci numbers are produced by the nth power of the transformation T, or Tⁿ, starting with the pair (1, 0).

We are then asked to consider the family of transformations T_pq which transforms the pair (a, b) according to the following rules:

a ← bq + aq + ap
b ← bp + aq

We can verify by quick substitution that the original transformation T is just a special case of T_pq, where p = 0 and q = 1.

a ← b(1) + a(1) + a(0)
a ← b + a

b ← b(0) + a(1)
b ← a

We are asked to show that if we apply 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 will give us an explicit way to square these transformations, which we can use to compute Tⁿ using successive squaring, just like the fast-expt procedure from exercise 1.16.

We can apply T_pq twice by defining new variables and using substitution. Let's define a₁ and b₁ as the results of applying transformation T_pq once

a₁ = bq + aq + ap
b₁ = bp + aq

The next step is to define a₂ and b₂ and apply the tranformation a second time, this time using a₁ and b₁ in place of a and b.

a₂ = b₁q + a₁q + a₁p
b₂ = b₁p + a₁q

Now that we have a system of equations defined, we can use substitution on our way to simplifying.

a₂ = (bp + aq)q + (bq + aq + ap)q + (bq + aq + ap)p
b₂ = (bp + aq)p + (bq + aq + ap)q

The second equation is shorter, so it should be easier to manipulate. Remember, we're trying to find p' and q', so we need to rewrite the equation to fit the form

b₂ = bp' + aq'

where p' and q' can be computed in terms of q and p.

b₂ = (bp + aq)p + (bq + aq + ap)q
= (bpp + apq) + (bqq + aqq + apq)
= bpp + apq + bqq + aqq + apq
= (bpp + bqq) + (2apq + aqq)
= b(pp + qq) + a(2qp + qq)

From this we can see that p' and q' can be computed using the following equations:

p' = p² + q²
q' = 2pq + q²

Manipulating the equation for a₂ in the same way, we can verify those results. This time we're trying to fit the form

a₂ = bq' + aq' + ap'

The groupings required for this manipulation are made even easier by the fact that we now already know the formulas for p' and q'.

a₂ = (bp + aq)q + (bq + aq + ap)q + (bq + aq + ap)p
= (bpq + aqq) + (bqq + aqq + apq) + (bpq + apq + app)
= bpq + aqq + bqq + aqq + apq + bpq + apq + app
= (bpq + bpq + bqq) + (apq + apq + aqq) + (app + aqq)
= b(pq + pq + qq) + a(pq + pq + qq) + a(pp + qq)
= b(2pq + qq) + a(2pq + qq) + a(pp + qq)

Now that we've verified the formulas for p' and q' in terms of p and q, we can use them to complete the procedure we were given.

(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
                    (+ (* p p) (* q q))     ; compute p'
                    (+ (* 2 p q) (* q q))   ; compute q'
                    (/ count 2)))
         (else (fib-iter (+ (* b q) (* a q) (* a p))
                         (+ (* b p) (* a q))
                         p
                         q
                         (- count 1)))))

You can test several known values in an interpreter to verify that this actually works.

> (fib 0)
0
> (fib 1)
1
> (fib 2)
1
> (fib 3)
2
> (fib 5)
5
> (fib 10)
55
> (fib 20)
6765

Related:

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