Project Nayuki


Fast Fibonacci algorithms

Definition: The Fibonacci sequence is defined as \(F(0) = 0\), \(F(1) = 1\), and \(F(n) = F(n-1) + F(n-2)\) for \(n ≥ 2\). So the sequence (starting with \(F(0)\)) is 0, 1, 1, 2, 3, 5, 8, 13, 21, ….

If we want to compute a single term in the sequence (e.g. \(F(n)\)), there are a couple of algorithms to do so. Some algorithms are much faster than others.

Algorithms

Textbook recursive (extremely slow)

Naively, we can directly execute the recurrence as given in the mathematical definition of the Fibonacci sequence. Unfortunately, it’s hopelessly slow: It uses \(Θ(n)\) stack space and \(Θ(φ^n)\) arithmetic operations, where \(φ = \frac{\sqrt{5} + 1}{2}\) (the golden ratio). In other words, the number of operations to compute \(F(n)\) is proportional to the final numerical answer, which grows exponentially.

Dynamic programming (slow)

It should be clear that if we already computed \(F(k-2)\) and \(F(k-1)\), then we can add them to get \(F(k)\). Next, we add \(F(k-1)\) and \(F(k)\) to get \(F(k+1)\). We repeat until we reach \(k = n\). Most people notice this algorithm automatically, especially when computing Fibonacci by hand. This algorithm takes \(Θ(1)\) space and \(Θ(n)\) operations.

Matrix exponentiation (fast)

The algorithm is based on this innocent-looking identity (which can be proven by mathematical induction):

\( \left[ \begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix} \right]^n = \left[ \begin{matrix} F(n+1) & F(n) \\ F(n) & F(n-1) \end{matrix} \right] \).

It is important to use exponentiation by squaring with this algorithm, because otherwise it degenerates into the dynamic programming algorithm. This algorithm takes \(Θ(1)\) space and \(Θ(\log n)\) operations. (Note: We are counting the number of bigint arithmetic operations, not fixed-width word operations.)

Fast doubling (faster)

Given \(F(k)\) and \(F(k+1)\), we can calculate these:

\(\begin{align} F(2k) &= F(k) \left[ 2F(k+1) - F(k) \right]. \\ F(2k+1) &= F(k+1)^2 + F(k)^2. \end{align}\)

These identities can be extracted from the matrix exponentiation algorithm. In a sense, this algorithm is the matrix exponentiation algorithm with the redundant calculations removed. It should be a constant factor faster than matrix exponentiation, but the asymptotic time complexity is still the same.

Summary: The two fast Fibonacci algorithms are matrix exponentiation and fast doubling, each having an asymptotic complexity of \(Θ(\log n)\) bigint arithmetic operations. Both algorithms use multiplication, so they become even faster when Karatsuba multiplication is used. The other two algorithms are slow; they only use addition and no multiplication.

Source code

Implementations are available in multiple languages:

  • Java: FastFibonacci.java (all 3 algorithms, timing benchmark, runnable main program)

  • Python: fastfibonacci.py (fast doubling function only)

  • JavaScript: fastfibonacci.js (fast doubling function only)

  • Haskell: fastfibonacci.hs (fast doubling function only)

  • C#: FastFibonacci.cs (fast doubling only, runnable main program)
    (requires .NET Framework 4.0 or above; compile with csc /r:System.Numerics.dll fastfibonacci.cs)

Benchmarks

Graphs

Benchmark graph

All algorithms, naive multiplication

Benchmark graph

All algorithms, Karatsuba multiplication

Benchmark graph

Fast algorithms, both multiplication algorithms

(Note: The graphs have logarithmic scales on the x and y axes.)

Table

n Fast doubling, Karatsuba multiplication Fast matrix, Karatsuba multiplication Fast doubling, naive multiplication Fast matrix, naive multiplication Slow dynamic programming Slow recursive
15 4141 0424 197887104
25 6382 0924 4421 8225322
35 7082 7404 5092 3429256
45 9453 0274 7332 660133114
55 9893 6774 7873 147172219
65 9723 9564 7653 371211400
86 1913 9724 9693 4282891 161
106 2834 9525 0224 1543703 113
136 3075 6105 0464 66748813 480
166 4794 9555 1774 21060557 300
206 5425 9235 2344 985763394 000
256 6326 5655 2635 4799644 373 000
326 7945 8875 3884 9081 235127 500 000
406 8186 8805 4335 7151 5525 980 000 000
506 8067 7425 4866 4462 023
636 93110 1805 5898 3392 598
797 16211 0905 7539 1873 396
1007 2799 2255 9047 7174 472
1267 42712 4106 05910 2205 866
1587 60013 0906 14110 9007 888
2008 00611 7006 5569 96910 640
2518 14615 6606 67213 06014 280
3168 59718 8107 08916 53019 610
3989 50120 5508 07818 12027 650
5019 96424 0508 49221 34038 970
63111 07038 7909 51035 72055 540
79413 02041 81011 52039 38080 280
1 00014 66050 87013 13048 230118 000
1 25918 64099 02016 99095 640175 300
1 58525 300113 50023 660110 800263 000
1 99532 360148 10030 770144 700397 500
2 51245 540314 80043 980311 400608 800
3 16267 800372 20066 250369 000937 200
3 98198 560491 50096 780488 1001 457 000
5 012143 5001 050 000145 9001 132 0002 269 000
6 310214 1001 284 000227 7001 357 0003 546 000
7 943320 3001 662 000351 3001 821 0005 547 000
10 000466 4003 519 000538 4004 382 0008 700 000
12 589691 1004 303 000851 7005 254 00013 640 000
15 8491 007 0005 481 0001 310 0007 079 00021 440 000
19 9531 493 00011 800 0002 081 00017 260 00033 620 000
25 1192 185 00013 620 0003 296 00020 710 00053 030 000
31 6233 205 00017 570 0005 159 00027 860 00083 310 000
39 8114 637 00036 800 0008 109 00068 540 000131 500 000
50 1196 750 00042 430 00012 910 00082 230 000207 700 000
63 0969 913 00054 770 00020 410 000110 600 000326 900 000
79 43314 450 000113 300 00032 300 000275 100 000517 100 000
100 00020 800 000130 600 00051 640 000330 700 000819 700 000
125 89330 380 000168 900 00081 150 000445 200 0001 296 000 000
158 48944 090 000346 800 000129 200 0001 103 000 0002 058 000 000
199 52663 260 000405 400 000205 100 0001 325 000 0003 249 000 000
251 18992 330 000517 300 000325 100 0001 766 000 0005 153 000 000
316 228133 700 0001 055 000 000515 700 0004 413 000 0008 161 000 000
398 107191 900 0001 228 000 000815 500 0005 311 000 00012 930 000 000
501 187280 200 0001 572 000 0001 297 000 0007 059 000 00020 520 000 000
630 957404 900 0003 181 000 0002 061 000 00017 570 000 00032 570 000 000
794 328580 700 0003 691 000 0003 265 000 00021 090 000 00051 650 000 000
1 000 000846 100 0004 724 000 0005 182 000 00028 310 000 00082 000 000 000
1 258 9251 221 000 0009 570 000 0008 168 000 00070 280 000 000130 300 000 000
1 584 8931 750 000 00011 050 000 00012 970 000 00084 120 000 000207 300 000 000
1 995 2622 549 000 00014 230 000 00020 610 000 000112 700 000 000329 700 000 000
2 511 8863 676 000 00028 800 000 00032 610 000 000279 900 000 000525 100 000 000
3 162 2785 247 000 00032 980 000 00051 600 000 000335 600 000 000
3 981 0727 654 000 000

All times are given in nanoseconds (ns) with 4 significant figures. All the tests above were performed on Intel Core 2 Quad Q6600 (2.40 GHz) using a single thread, Windows XP SP 3, Java 1.6.0_22.

Proofs

Matrix exponentiation

We will use weak induction to prove this identity.

Base case

For \(n = 1\), clearly \( \left[ \begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix} \right]^1 = \left[ \begin{matrix} F(2) & F(1) \\ F(1) & F(0) \end{matrix} \right] \).

Induction step

Assume for \(n ≥ 1\) that \( \left[ \begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix} \right]^n = \left[ \begin{matrix} F(n+1) & F(n) \\ F(n) & F(n-1) \end{matrix} \right] \). Then:

\( \left[ \begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix} \right]^{n+1} \\ = \left[ \begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix} \right]^n \left[ \begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix} \right] \\ = \left[ \begin{matrix} F(n+1) & F(n) \\ F(n) & F(n-1) \end{matrix} \right] \left[ \begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix} \right] \\ = \left[ \begin{matrix} F(n+1)+F(n) & F(n+1)+0 \\ F(n)+F(n-1) & F(n)+0 \end{matrix} \right] \\ = \left[ \begin{matrix} F(n+2) & F(n+1) \\ F(n+1) & F(n) \end{matrix} \right]. \)

Fast doubling

We will assume the fact that the matrix exponentiation method is correct for all \(n ≥ 1\).

\( \left[ \begin{matrix} F(2n+1) & F(2n) \\ F(2n) & F(2n-1) \end{matrix} \right] \\ = \left[ \begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix} \right]^{2n} \\ = \left( \left[ \begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix} \right]^n \right)^2 \\ = \left[ \begin{matrix} F(n+1) & F(n) \\ F(n) & F(n-1) \end{matrix} \right]^2 \\ = \left[ \begin{matrix} F(n+1)^2+F(n)^2 & F(n+1)F(n)+F(n)F(n-1) \\ F(n)F(n+1)+F(n-1)F(n) & F(n)^2+F(n-1)^2 \end{matrix} \right]. \)

Therefore, by equating the cells in the matrix:

\(\begin{align} F(2n+1) &= F(n+1)^2 + F(n)^2. \\ F(2n) &= F(n) \left[ F(n+1) + F(n-1) \right] \\ &= F(n) \left[ F(n+1) + (F(n+1) - F(n)) \right] \\ &= F(n) \left[ 2F(n+1) - F(n) \right]. \\ F(2n-1) &= F(n)^2 + F(n-1)^2. \end{align}\)

More info