SICP exercise 1.19, transformations

From Structure and Implementation of Computer Programs, second edition:

Exercise 1.19: 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ⁿ, the nth 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:(5)
(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)))))
(5) This exercise was suggested to us by Joe Stoy, based on an example in Kaldewaij 1990.

I don’t know where to start with how confused I am. T is a “special case” of a transformation that is totally different and involves two extra variables? (Well, it says it’s a “special case” of a “family of transformations,” and I’m guessing those are both mathematical terms of art of which I am unaware and for which I cannot find any definitions.)

Also, apparently if you use p and q to transform a and b once, according to a ← bq + aq + ap and b ← bp + aq, a becomes the next number in the Fibonacci sequence (while b holds a’s previous value), just like in the original fib-iter process, but if you apply the transformation twice, you can halve the number of steps it takes to get the desired value in the Fibonacci sequence?!

Given all this, I have absolutely no idea what p' and q' are supposed to be; I mean, judging from the text and the code, I have to assume that p' and q' are the values of p and q that, when used to transform a and b, would end up getting them halfway down the road to b becoming the desired value in the Fibonacci sequence (and a becoming the value after that), but not only do I have no idea how to start figuring out the method for computing p' and q', I also can’t understand how the transformation a ← bq + aq + ap and b ← bp + aq is supposed to be equivalent to a ← a + b and b ← a when you apply it once but something totally different when you apply it twice.

(I suppose “apply the transformation twice” must really mean something to do with “squaring the transformation,” of which I can find something online saying “T²(x)=T(T(x))”, which doesn’t really help here, as T_(pq) only transforms the pair (a,b) and I’m supposed to be figuring out how to calculate p' and q' from p and q… let me know if I’m barking up the wrong tree, here.)

Anyway, I suppose I must be missing some tools in my mental toolbox for understanding any of this, so please let me know what they are. Please DO NOT give me the answer to the exercise, or directly describe how to calculate p' and q'; I’m trying to get through this book legit. I just need to know what it is I’m supposed to understand and don’t. Thank you.

Solution

First addressing the elephant in the room of fib-iter being recursive, not iterative, as defined in your book (this is a programming forum after all).

The reason why it is a valid way of naming the function is because this is a tail recursion. All reasonable compilers/interpreters (all the ones I know at least) will transform the recursion into an iteration by tail call elimination.

Secondly, while the book claims the number of steps is logarithmic, the author does not mean the time complexity is logarithmic. The reason is that each of these steps is not O(1).
To be a bit more specific, fib(n) having O(n) digits, it is not possible to calculate them in time faster than O(n). In practice, most of that time is spent doing multiplications and the time complexity is greater than O(n).

The exercise presents an implementation of the problem: How many distinct ways are there to climb stairs made of n steps with the constraint that at each step, you can only climb 1 or 2 steps up.
I am simplifying a bit but considering the function stairs, returning the count:

Proving it corresponds to the Fibonacci sequence is easy: if you have n steps to climb, the first step you make is either:
- 1 step up, then stairs of length n-1 remain.
- Or 2 steps up, then stairs of length n-2 remain.
At the end of the day, stairs(n) = stairs(n-1) + stairs(n-2) (= fibonacci(n+1)).
From there, you can find a better way of making the calculation, e.g. by counting the ways to climb the first half, then second half of the stairs. As paths chosen on both halves are independent, the result is a product.
Again, I am simplifying quite a bit; the formula depends on whether n is even or odd and must account for paths that land the middle step or skip it.
I leave it to you to figure out both formulas and you should find the same as what the answer for this exercise.

The implementation is clever in the sense that it is implemented with an iteration, like I said above, rather than pure recursion.

Back to the problem itself:

If T_pq(a,b) does: a ← bq + aq + ap, b ← bp + aq,
it should be easy to see how T₀₁(a,b) does: a ← 1b + 1a + 0a = a + b, b ← 0b + 1a = a, which is exactly what T does.

Now, you are right that T²(x) = T(T(x)) here, meaning we are to calculate T_pq(T_pq(a,b)):

a ← qT_pq(a,b)_b + qT_pq(a,b)_a + pT_pq(a,b)_a
b ← pT_pq(a,b)_b + qT_pq(a,b)_a

or, when applying T_pq:

a ← q(bp + aq) + q(bq + aq + ap) + p(bq + aq + ap)
b ← p(bp + aq) + q(bq + aq + ap)

In 2 steps, we can easily expand both lines, then reorganize them to find an expression of the same form, with p' and q':

a ← bpq + aqq + bqq + aqq + apq + bpq + apq + app, meaning:
a ← b(2pq + qq) + a(2pq + qq) + a(pp+qq)
b ← bpp + apq + bqq + aqq + apq, meaning
b ← b(pp+qq) + a(2pq + qq)

and so:

p' = (pp+qq) and q' = (2pq + qq)

Knowing that lets you do fast exponantiation i.e.

Calculate p'' and q'' so that T_p''q'' = T_p'q'² = T_pq⁴
then:
Calculate p''' and q''' so that T_p'''q''' = T_pq⁸

... and so on

You can indeed use these 2 expressions to complete the function and it gives the right results.

The function:
- keeps b,a and p,q as 2 pairs of consecutive fibonacci numbers (in that order)
And, depending on the parity of count:
- Ensures q is some fib(2^k) (q is first fib(1), then fib(2), fib(4), fib(8) and so on thanks to the fast exponantiation described earlier), meaning p is some fib(2^k-1)
- Advances a and b by 2^k positions in the fibonacci sequence. (e.g. if b = fib(35) and q = fib(128) then b becomes fib(163)).

If you are interested in testing the performance of that function, below is a C++ implementation that uses a big integer library for its calculation (A difference only shows for big values of n/fib(n)), to be compared with fib_linear, an implementation of applying T(a,b) n times (O(n) steps, but time complexity is O(n²)).

I have included a line (to be uncommented) to let you track the values of b, a, p and q.

#include <chrono>
#include <iostream>
#include "bigint.h"

using bigint::bigint_t;
using std::chrono::duration_cast;
using std::chrono::milliseconds;
using std::chrono::steady_clock;
using std::chrono::time_point;

bigint_t fib_iter(bigint_t a, bigint_t b, bigint_t p, bigint_t q, unsigned int count)
{
    //std::cerr << "b = " << b << '\n' <<  "a = " << a << '\n' << "p = " << p << '\n' << "q = " << q << '\n' << '\n';
    if (count == 0)
        return b;

    if (count % 2 == 0) {
        auto p2 = p * p, q2 = q * q, pq2 = (p * q) << 1;
        return fib_iter(a, b, p2 + q2, q2 + pq2, count / 2);
    }
    else {
        auto ap = a * p, aq = a * q, bp = b * p, bq = b * q;
        return fib_iter(
            bq + aq + ap,
            bp + aq,
            p,
            q,
            count - 1
        );
    }
}
bigint_t fib(unsigned int n)
{
    return fib_iter(1, 0, 0, 1, n);
}

bigint_t fib_linear(unsigned int n)
{
    if (n <= 1)
        return bigint_t(static_cast<int>(n));
    bigint_t a = 1, b = 0;
    for (unsigned int i = 1; i < n; ++i) {
        auto c = a;
        a = a + b;
        b = c;
    }
    return a;
}

int main(int argc, char *argv[])
{
    unsigned int n = 100000;

    steady_clock clock;
    auto start = clock.now();   
    auto f = fib_linear(n);
    auto mid = clock.now();
    bigint_t f2 = fib(n);
    auto end = clock.now();

    if (f == f2) {
        auto durationOld = duration_cast<milliseconds>(mid - start).count(), durationNew = duration_cast<milliseconds>(end - mid).count();
        std::cerr << durationNew << "ms" << " vs " << durationOld << "ms" << (durationNew < durationOld ? " (faster than original)" : " (faster than original)") << '\n';
    }
    else
        std::cerr << "Incorrect result" << '\n';
    
    std::cerr << f.toString() << '\n';
    std::cerr << f2.toString() << '\n';

    return 0;
}