# How does the Fibonacci exponentiation by squaring algorithm work?

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;
}
``````

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.

• 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). – rcgldr Feb 20 at 8:44

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
``````
• 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. – High Performance Mark Feb 20 at 16:42
• @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. – rcgldr Feb 20 at 17:39