## Friday, January 22, 2010

### 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 Tn, starting with the pair (1, 0).

We are then asked to consider the family of transformations Tpq 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 Tpq, 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 Tpq twice, the effect is the same as using a single transformation Tp'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 Tn using successive squaring, just like the fast-expt procedure from exercise 1.16.

We can apply Tpq twice by defining new variables and using substitution. Let's define a1 and b1 as the results of applying transformation Tpq once

a1 = bq + aq + ap
b1 = bp + aq

The next step is to define a2 and b2 and apply the tranformation a second time, this time using a1 and b1 in place of a and b.

a2 = b1q + a1q + a1p
b2 = b1p + a1q

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

a2 = (bp + aq)q + (bq + aq + ap)q + (bq + aq + ap)p
b2 = (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

b2 = bp' + aq'

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

b2 = (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' = p2 + q2
q' = 2pq + q2

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

a2 = 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'.

a2 = (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.

Eric Hsu said...

I suppose it's missing the point to use the closed form of the Fibonacci numbers...?

Bill the Lizard said...

Eric,
Yes, definitely. The first code sample I showed in the article is what's given in the book, and all we have to do is fill in the blanks.

I am planning a separate article on the closed form (Binet's formula) though, since it doesn't seem to be mentioned explicitly in the book. (The authors did lead us very close to it in exercise 1.13, though.)

Alfred Hitchcook said...

You might love this http://www.catonmat.net/blog/using-fibonacci-numbers-to-convert-from-miles-to-kilometers/

Martin Schönert said...

Your post is of course quite old, but anyhow ;-).

My copy of SICP is apparently too old and does not include this excercise. So I don't know whether the authors mention the following fact.

The logarithmic algorithm (which is a special case of a more general algorithm to compute Lucas sequences) is optimal in the sense that the number of bit operations required is proportional to the size of the result.

The same can be achieved for Binet's formula, but it is a bit tricky working out the required precision for each step (so that you have enough precision to get the correct result but not too much and ruin the bit complexity).

Bartosz Prokop said...

That is way too much unnecessary math. All you have to see is that applying your transformation T^n that is defined by p and q to those parameters will effectively square it. This is hinted by the authors (successive squaring) and I think it's overall better solution then your brute-force one.