How does the Fibonacci exponentiation by squaring algorithm work?

Question

This is one of the best algorithms to calculate the nth Fibonacci sequence. It only needs O(log(n)) time to do its job, so it's very efficient. I found it somewhere but don't know how it works!
Can anyone tell me how this algorithm works?

int fib3 (int n) {
 int i = 1, j = 0, k = 0, h = 1, t;
 while (n > 0) {
 if (n % 2) {
 t = j * h;
 j = i * h + j * k + t;
 i = i * k + t;
 }
 t = h * h;
 h = 2 * k * h + t;
 k = k * k + t;
 n /= 2;
 }
 return j;
}

• Robert Harvey: See en.wikipedia.org/wiki/Exponentiation_by_squaring, and nayuki.eigenstate.org/page/fast-fibonacci-algorithms under the Matrix exponentiation (medium) heading.
• Rasool: Thanks @RobertHarvey, I saw that link but I didn't understand. Could you please tell me what exactly it does? step-by-step. Thanks a lot :)
• Robert Harvey: Well, the code probably corresponds to the math. You don't understand the math?
• Rasool: No I understand math. I just can't relate the formula to this code. and I don't know if I've got the formula matching exactly this code! @RobertHarvey

Answer 1

It's called the "matrix form" - take a look at Wikipedia

You can compute next Fibonacci number (k+2) by multiplying matrix on a vector of two previous elements (k + 1 and k). Hence, k + 3 can be computed by multiplying matrix on vector of (k + 2 and k + 1). This equals squared matrix multiplied on (k + 1 and k). So on.

Your code simply squares the matrix, taking into account odd powers.

• rcgldr: A true matrix form would require 12 variables. This is a variation of a Lucas sequence, which takes 5 variables (as seen in the question).

Answer 2

Late answer, but an explanation for how this works:

The algorithm is based on Lucas sequence relations for Fibonacci numbers.

https://en.wikipedia.org/wiki/Lucas_sequence#Other_relations

Specifically these equations:

F(m) = F(m-1) + F(m-2)
F(m+n) = F(m+1) F(n) + F(m) F(n-1)
F(2n) = F(n) L(n) = F(n) (F(n+1) + F(n-1))
 = F(n)((F(n) + F(n-1)) + F(n-1))
 = F(n) F(n) + 2 F(n) F(n-1)

Initial state:

i = F(-1) = 1
j = F( 0) = 0
k = F( 0) = 0
h = F( 1) = 1

n is treated as the sum of powers of 2: 2^a + 2^b + ... for each iteration e, let p = 2^i, then

h = F(p)
k = F(p-1)

To advance to the next iteration, h and k are advanced to F(next power of 2):

h' = F(2p) = F(p) F(p+1) + F(p) F(p-1)
 = F(p)(F(p) + F(p-1)) + F(p) F(p-1)
 = F(p) F(p) + F(p) F(p-1) + F(p) F(p-1)
 = F(p) F(p) + 2 F(p) F(p-1)
 = h h + 2 h k
k' = F(2p-1) = F(p + (p-1)) = F(p+1) F(p-1) + F(p) F(p-2)
 = (F(p) + F(p-1)) F(p-1) + F(p) (F(p) - F(p-1))
 = F(p) F(p-1) + F(p-1) F(p-1) + F(p) F(p) - F(p) F(p-1)
 = F(p) F(p) + F(p-1) F(p-1)
 = h h + k k

During the calculation of i and j, let c = current cumulative sum of bits of n:

i = F(c-1)
j = F(c)

To update i and j for 1 bits in n corresponding to p = current power of 2:

i' = F((c-1)+p) = F(c) F(p) + F(c-1) F(p-1)
 = j h + i k
j' = F(c+p) = F(c+1) F(p) + F(c) F(p-1)
 = (F(c)+F(c-1)) F(p) + F(c) F(p-1)
 = F(c) F(p) + F(c-1) F(p) + F(c) F(p-1)
 = j h + i h + j k

• High Performance Mark: This is a very good answer, but it's not a very good Software Engineering answer. Probably just shows that the question was asked in the wrong place.
• rcgldr: @HighPerformanceMark - it's a question about an algorithm, in this case. how it works (similarly to asking how it was derived, which is what my answer provides), and I see "related" questions about algorithms also in Software Engineering. I spend more time at SO than at SE, so I'm not that familiar with where questions belong at SE.