[sbi]

Thank you for the reference. The algorithm in the parent comment is not actually matrix exponentiation but the "fast doubling algorithm" in Minase's terminology, where you exponentiate the matrix algebraically using the Cayley-Hamilton theorem rather than using matrix multiplication.

[ColinWright]

You have me confused. The usual "Russian Peasant Multiplication" can be adapted to produce fast exponentiation, or "exponentiation by squaring". Your code appears to be doing exactly that.

Further, that is then exactly what is quoted in reference [11], just as leephillips says. There it is then taken a step further and specialized to the calculation of the Fibonacci numbers.

So the referenced paper seems to cover what you have said, and go even further. To that end, I'm confused as to what you are saying. Are you agreeing that the paper does what your code does? Or are you disagreeing, in which case, could you be more explicit?

Thanks.

[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.