Finding Nth Fibonacci number in logarithmic time

The obvious way to compute $F_n$ is to add your way up from the bottom, which costs $O(n)$. There is a better way, and it comes from noticing that Fibonacci is secretly a matrix power.

Where the matrix comes from

The Fibonacci numbers appear as the successive convergents of the continued fraction for the golden ratio $\varphi$:

$\varphi=1+\cfrac1{1+\cfrac1{1+\cfrac1{1+\cfrac1{1+\ddots}}}}$

The matrix formed from successive convergents of any continued fraction has determinant $\pm1$, and for this one that matrix is $\begin{pmatrix}1&1\\\ 1&0\end{pmatrix}$. Raising it to the $n^{th}$ power lays out three consecutive Fibonacci numbers:

${\begin{pmatrix}1&1\\\ 1&0\end{pmatrix}}^n=\begin{pmatrix}F_{n+1}&F_n\\\ F_n&F_{n-1}\end{pmatrix}$

So computing $F_n$ reduces to computing a matrix power, and the answer falls out of position $(0,0)$.

That alone buys nothing. Multiplying the matrix by itself one step at a time is $n$ multiplications — the same $O(n)$ as just adding, with more bookkeeping.

Fast exponentiation

The saving comes from how you take the power. To compute $M^n$, compute $M^{n/2}$ once and square it, throwing in one extra multiplication when $n$ is odd:

$M^n=\begin{cases}(M^{n/2})^2 & n \text{ even}\\\ (M^{n/2})^2\times M & n \text{ odd}\end{cases}$

Each step halves $n$, so there are $O(\log n)$ multiplications instead of $n$. Since each matrix multiply is a fixed 8 multiplications and 4 additions, the whole thing is $O(\log n)$.

Implementation

#include <stdio.h>
 
void multiply(int F[2][2], int M[2][2]);
void power(int F[2][2], int n);
 
int fib(int n) {
	int F[2][2] = {{1, 1}, {1, 0}};
	if (n == 0)
		return 0;
	power(F, n - 1);
	return F[0][0];
}
 
void power(int F[2][2], int n) {
	if (n == 0 || n == 1)
		return;
	int M[2][2] = {{1, 1}, {1, 0}};
	power(F, n / 2);
	multiply(F, F);
	if (n % 2 != 0)
		multiply(F, M);
}
 
void multiply(int F[2][2], int M[2][2]) {
	int x = F[0][0] * M[0][0] + F[0][1] * M[1][0];
	int y = F[0][0] * M[0][1] + F[0][1] * M[1][1];
	int z = F[1][0] * M[0][0] + F[1][1] * M[1][0];
	int w = F[1][0] * M[0][1] + F[1][1] * M[1][1];
	F[0][0] = x;
	F[0][1] = y;
	F[1][0] = z;
	F[1][1] = w;
}
 
int main() {
	printf("%d\n", fib(15));
	return 0;
}

610

The four highlighted lines are the entire optimisation. power(F, n / 2) recurses on half the exponent, multiply(F, F) squares the result, and the odd case folds in one more copy of the base matrix.

Two details are easy to get wrong. fib calls power(F, n - 1) rather than power(F, n), because F already starts as $M^1$ rather than the identity. And multiply writes into four temporaries before assigning, since overwriting F[0][0] early would corrupt the values the later expressions still need to read.