[sbi]

You are correct that this code uses "Russian Peasant Multiplication." The question is: what is being exponentiated? One choice is the matrix A = [0 1; 1 1], whose powers are A^n = [F_{n-1} F_n; F_n F_{n+1] if we compute those powers treating A as a matrix. This involves some redundant computations since F_n is being computed twice.

An alternative method is to use the following identity:

  (aA + b)(cA + d) = acA^2 + (ad + bc)A + bd
                   = (ad + bc + ac)A + (bd + ac)
                   = ((a + b)(c + d) - bd)A + (bd + ac)

(In fact, A^n = F_n A + F_{n-1} I_2). This step requires 3 multiplications and 4 additions and is more-or-less equivalent to the "fast doubling" algorithm in reference [11], which is more optimized. But this derivation is completely generic---it just relies on the characteristic polynomial of A, namely A^2 - A - 1 = 0.

In the case where A is not the Fibonacci matrix but something much larger, working with polynomials modulo the characteristic polynomial of A is much faster than working with A as a matrix.

[ColinWright]

Right - that helps. Thanks for the clarification - very helpful.