Fibonacci Computation, in depth: fast doubling, big-integer cost and modular answers

Definition and growth

The sequence is defined by F(0) = 0, F(1) = 1 and F(n) = F(n-1) + F(n-2). The first values are 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55. The ratio of consecutive terms approaches the golden ratio phi = (1 + sqrt 5) / 2, about 1.618, so F(n) grows like phiⁿ/sqrt 5. The number of decimal digits of F(n) is therefore about n times log₁₀(phi), roughly 0.209n. F(1,000,000) has 208,988 digits. That growth rate is the key to every cost argument below.

Naive recursion and the loop

The direct translation of the definition calls itself twice per call. If C(n) is the number of calls made to compute F(n), then C(n) = C(n-1) + C(n-2) + 1 with C(0) = C(1) = 1, which solves to C(n) = 2F(n+1) - 1. Computing F(40) makes 331,160,281 calls. The work grows as fast as the answer itself, because the recursion recomputes the same subproblems over and over.

Memoisation, or the bottom-up loop that dynamic programming would give you, removes the repetition and needs only two variables:

def fib_naive(n):            # exponential: never use beyond toy n
    return n if n < 2 else fib_naive(n - 1) + fib_naive(n - 2)

def fib_loop(n):             # n additions, O(1) extra space
    a, b = 0, 1
    for _ in range(n):
        a, b = b, a + b
    return a

The loop performs n additions. With fixed-width integers that is O(n) time and the right answer for small n. With arbitrary-precision integers it is not O(n): the k-th addition works on numbers of about 0.7k bits, so the total is about the sum of k over all k, which is O(n²) bit operations. For n = 1,000,000 that is noticeable; for n = 100,000,000 the loop is hopeless.

From the matrix to the doubling identities

The matrix form packs the recurrence into a 2-by-2 power. With Q = [[1, 1], [1, 0]], Qⁿ = [[F(n+1), F(n)], [F(n), F(n-1)]]. Raising Q to the n-th power by repeated squaring needs O(log n) matrix products, which is the general technique in matrix exponentiation. Fibonacci has extra structure that makes it cheaper. Writing Q^2k = Q^k Q^k and reading off the entries gives the doubling identities:

F(2k)   = F(k) * (2*F(k+1) - F(k))
F(2k+1) = F(k)^2 + F(k+1)^2

From the pair (F(k), F(k+1)) these give (F(2k), F(2k+1)) using three multiplications. One more addition steps to (F(2k+1), F(2k+2)). A full 2-by-2 matrix squaring, by contrast, needs eight multiplications naively or several fewer with symmetry tricks, so fast doubling does the same logarithmic number of steps with a smaller constant and far less bookkeeping.

Implementing fast doubling

The recursive form follows the identities directly. The iterative form scans the bits of n from the most significant end, which avoids recursion entirely and is easy to port to languages without big stacks.

def fib_pair(n):
    """Return (F(n), F(n+1)) by recursive fast doubling."""
    if n == 0:
        return 0, 1
    a, b = fib_pair(n >> 1)          # a = F(k), b = F(k+1), k = n // 2
    c = a * (2 * b - a)              # F(2k)
    d = a * a + b * b                # F(2k+1)
    return (d, c + d) if n & 1 else (c, d)

def fib(n):
    """Iterative fast doubling, most significant bit first."""
    a, b = 0, 1                      # (F(0), F(1))
    for bit in bin(n)[2:]:
        c = a * (2 * b - a)
        d = a * a + b * b
        a, b = (d, c + d) if bit == "1" else (c, d)
    return a

assert [fib(i) for i in range(11)] == [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
assert all(fib(i) == fib_pair(i)[0] for i in range(500))

Worked example: F(13) by bits

Follow the iterative version for n = 13, binary 1101, as in the diagram. Start at k = 0 with (0, 1). The first bit is 1: doubling k = 0 gives (0, 1) again, and stepping gives k = 1, (1, 1). The second bit is 1: doubling gives F(2) = 1 times (2 - 1) = 1 and F(3) = 1 + 1 = 2; stepping gives k = 3, (2, 3). The third bit is 0: doubling gives F(6) = 2 times (6 - 2) = 8 and F(7) = 4 + 9 = 13, so k = 6 with (8, 13). The last bit is 1: doubling gives F(12) = 8 times (26 - 8) = 144 and F(13) = 64 + 169 = 233, and stepping gives k = 13, (233, 377). F(13) = 233, after four rounds instead of thirteen additions.

The real cost: big-integer multiplication

With big integers, the number of steps is the wrong thing to count. Fast doubling does O(log n) rounds, but the last round multiplies numbers of about n/2 times 0.69 bits, and each round's operands are twice the size of the previous round's. The total cost is dominated by the final few multiplications, so the whole computation costs about a constant times M(n), the cost of one multiplication of n-bit numbers. That makes the integer library the real algorithm choice:

• CPython's built-in int uses schoolbook multiplication for small numbers and Karatsuba, about O(n^1.585), for large ones.
• GMP, available in Python through gmpy2, switches to Toom-Cook and FFT-based methods as sizes grow, which is much faster for millions of digits.
• Java's BigInteger also switches algorithms by size; benchmark on your JDK rather than assuming.

A practical consequence: for exact huge values, converting the result to a decimal string can cost as much as computing it. Time the computation and the printing separately.

import time
for n in (10**5, 10**6, 10**7):
    t = time.perf_counter(); x = fib(n); t1 = time.perf_counter() - t
    print(f"n={n:>9,}  bits={x.bit_length():>10,}  compute={t1:.3f}s")

Fibonacci modulo m and Pisano periods

Many problems, from competitive programming to hashing, ask for F(n) modulo m with n as large as 10¹⁸. Exact values are out of the question, but fast doubling works unchanged if every operation is reduced modulo m, and now each multiplication is constant time, so the whole thing is O(log n). In languages with fixed-width integers, keep intermediate products below overflow: with m near 10⁹, a product fits in 64 bits; with m near 10¹⁸, use 128-bit multiplication. In C, % of a negative number is negative, so compute 2*b - a as (2*b - a + m) % m. The same reduction pattern appears in modular exponentiation.

The sequence modulo m is periodic, and its period is the Pisano period pi(m). For example pi(10) = 60, so the last digit of F(n) is the last digit of F(n mod 60). The period is at most 6m, so for small m you can find it by searching for the pair (0, 1) to reappear:

def fib_mod(n, m):
    a, b = 0, 1
    for bit in bin(n)[2:]:
        c = a * ((2 * b - a) % m) % m
        d = (a * a + b * b) % m
        a, b = (d, (c + d) % m) if bit == "1" else (c, d)
    return a

def pisano(m):
    a, b = 0, 1
    for i in range(1, 6 * m + 1):
        a, b = b, (a + b) % m
        if (a, b) == (0, 1):
            return i

assert pisano(10) == 60
assert fib_mod(10**18, 10) == fib_mod(10**18 % 60, 10)

Floating point and fixed-width limits

Binet's formula, F(n) = (phiⁿ - psiⁿ)/sqrt 5 with psi = (1 - sqrt 5)/2, is exact in real arithmetic and tempting in code as round(phi**n / sqrt(5)). In IEEE-754 double precision it is first wrong at n = 71 (measured on CPython 3.13; the script below finds the point on your platform), because doubles carry about 53 bits and the error in phiⁿ grows with n. Never use it for exact values.

from math import sqrt
phi = (1 + sqrt(5)) / 2
a, b, n = 0, 1, 0
while round(phi**n / sqrt(5)) == a:
    a, b, n = b, a + b, n + 1
print("Binet with doubles first fails at n =", n)

The last row surprises people. Since Python 3.11 (and security releases of some earlier versions), converting an integer of more than 4,300 decimal digits to a string raises ValueError by default, so print(fib(10**6)) fails even though the computation succeeded. Raise the limit deliberately with sys.set_int_max_str_digits(), or keep the number in binary and write it with hex() or to_bytes() if decimal output is not required.

Failure modes

• Overflow without an error. Java and C wrap silently. In Java use Math.addExact and Math.multiplyExact to make overflow throw, or switch to BigInteger.
• Recursion depth. Naive recursion and memoised recursion hit Python's default limit near n = 1,000. Recursive fast doubling only recurses log₂(n) deep, so it is safe.
• Negative intermediate under a modulus. 2*b - a can be negative in C, Java and Go before reduction.
• Benchmarking the printer. Timing print(fib(n)) measures decimal conversion as much as Fibonacci.
• Off-by-one indexing. Some sources start at F(1) = F(2) = 1 and some at F(0) = 0. Pin the convention in a test with known values.
• Trusting floats. Binet and logarithm tricks are fine for estimates such as digit counts, never for exact answers.

Trade-offs

If you need every value up to n, the loop is optimal, because you must produce them all anyway. If you need one value, use fast doubling. If the recurrence is not Fibonacci but another linear recurrence, use the matrix method. Reason about growth rates the way big-O analysis teaches, but count bit operations once numbers stop fitting in a register.

What to do next

1. Implement fib and fib_pair from this article and check them against the loop for n up to 1,000.
2. Run the timing script for n = 10⁵, 10⁶ and 10⁷; then repeat with gmpy2.mpz operands and compare.
3. Print F(10⁶) and handle the ValueError by raising the digit limit explicitly.
4. Write fib_mod in Java or C with 64-bit arithmetic and m = 10⁹ + 7, guarding negative intermediates.
5. Verify pisano(10) = 60 and use it to compute the last digit of F(10¹⁸) two ways.
6. Generalise: rewrite the doubling step as a 2-by-2 matrix squaring and apply it to a recurrence such as tribonacci with a 3-by-3 matrix.