Is it possible to use the closed-form of Fibonacci series to generate the Nth Fibonacci number exactly and efficiently?

The closed-form of the Fibonacci series is the following:

As you can see the expression contains square roots so we cannot use it directly to generate Nth Fibonacci number exactly, as sqrt(5) is irrational, we can only approximate it even in pure mathematics, and floating point inaccuracies introduces a whole other lot of complications that ensure if we implement the formula naively we are bound to lose accuracy.

I want to get the Nth Fibonacci number exactly using this formula, assume N is very large (in the thousands range), and since Python's integer is unbounded and Fibonacci sequence grows extremely quickly out of the range of uint64 (at the 94th term to be exact, we start at 0th), I cannot use Numba or NumPy to vectorize this.

Now I have implemented a function that generates Nth Fibonacci number exactly using the closed form, but it is very inefficient.

First I will just quickly describe what I have done, we start with (1 + a)ⁿ, using the binomial theorem we get 1 + na + _nC₂a² + _nC₃a³ + _nC₄a⁴ ... + aⁿ. Because the choose function is symmetric we only need to compute half of the terms.

Now we can compute _nC_i from the previous term _nC_i-1: _nC_i = _nC_i-1 * (n + 1 - i) / i, this is much faster than using factorials.

Then we expand (1 - a)ⁿ: 1 - na + _nC₂a² - _nC₃a³ + _nC₄a⁴ - _nC₅a⁵ + _nC₆a⁶...

The signs alternate, now if we evaluate (1 + a)ⁿ - (1 - a)ⁿ, we find the even powers cancel out, leaving only odd powers, and since we divide by sqrt(5) we are left with only powers of 5.

---

Code:

def fibonacci(lim: int) -> list[int]:
    a, b = 1, 1
    result = [0] * lim
    for i in range(1, lim):
        result[i] = a
        a, b = b, a + b

    return result


def Fibonacci_phi_even(n: int) -> int:
    fib = 0
    power = 1
    a, b, c = n - 1, 2, 2 * n
    length = int(n / 4 + 0.5)
    coeffs = [1] * length
    for i in range(length):
        coeffs[i] = c
        fib += c * power
        power *= 5
        c = c * (a - 1) * a // ((b + 1) * b)
        a -= 2
        b += 2

    for coeff in coeffs[length - 1 - (n // 2 & 1) :: -1]:
        fib += coeff * power
        power *= 5

    return fib >> n


def Fibonacci_phi_odd(n):
    length = n // 2 + 1
    powers = [1] * length
    power = 1
    for i in range(length):
        powers[i] = power
        power *= 5

    fib = 0
    a, b, c, d = n, 1, 2, -1
    i = di = length - 1
    for _ in range(length):
        fib += c * powers[i]
        c = c * a // b
        a -= 1
        b += 1
        i += d * di
        d *= -1
        di -= 1

    return fib >> n


def Fibonacci_phi(n: int) -> int:
    if n <= 2:
        return int(n > 0)

    return Fibonacci_phi_odd(n) if n & 1 else Fibonacci_phi_even(n)

It works:

In [377]: print(fibonacci(25))
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368]

In [378]: print([Fibonacci_phi(i) for i in range(25)])
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368]

But it is slow:

In [379]: %timeit fibonacci(1025)
87.9 μs ± 633 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [380]: %timeit Fibonacci_phi(1024)
699 μs ± 8.42 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

How can we make it faster without losing accuracy?

I have done it, I have made the closed form spit out Nth Fibonacci number exactly efficiently.

Now instead of doing all that inefficient binomial expansion, we can reduce the problem to the most basic part. Since we are doing exponentiation, we are really just doing multiplications repeatedly, we first need to find a way to multiply numbers like (a + b√5) exactly.

Now using simple algebra, we have the following relationship:

(a + b√5) * (c + d√5) = ac + bc√5 + ad√5 + 5bd

Now what do we do? We group the rational terms and irrational terms, and then we get rid of the radical, we simplify the above to this:

(a, b) * (c, d) = (ac + 5bd, ad + bc)

The above expression doesn't use floating points anywhere and is guaranteed to be exact.

Now what do we do next? We just multiply the numbers by itself N times to do exponentiation. But of course that is inefficient, since we have defined multiplication, we can use exponentiation by squaring.

Now we are computing (1, 1)ⁿ - (1, -1)ⁿ, the rational part will get cancelled out, the irrational part will have opposite signs and so we only need to compute (1, 1)ⁿ and double the irrational part, and we then right shift by n to get nth Fibonacci number.

def poly_mult(t1: tuple[int, int], t2: tuple[int, int]) -> tuple[int, int]:
    (a, b), (c, d) = t1, t2
    return (a * c + 5 * b * d, a * d + b * c)


def exp_by_sqr(base: tuple[int, int], exp: int) -> tuple[int, int]:
    prod = (1, 0)
    if not exp:
        return prod
    
    while exp > 1:
        if exp & 1:
            prod = poly_mult(base, prod)
            exp -= 1
        
        base = poly_mult(base, base)
        exp >>= 1
    
    return poly_mult(base, prod)


def Fibonacci_phi_poly(n: int) -> int:
    return (2 * exp_by_sqr((1, 1), n)[1]) >> n

For comparison, the fast doubling approach:

def fibonacci_fast_doubling(n: int) -> tuple[int, int]:
    if not n:
        return (0, 1)

    a, b = fibonacci_fast_doubling(n >> 1)
    c = (a2 := a * a) + (b2 := b * b)
    d = 2 * a * b
    return (c, b2 + d) if n & 1 else (d - a2, c)

In [10]: print([Fibonacci_phi_poly(i) for i in range(26)])
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025]

In [11]: %timeit fibonacci(1024)
90 μs ± 664 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [12]: %timeit Fibonacci_phi(1023)
862 μs ± 5.74 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [13]: %timeit fibonacci_fast_doubling(1023)
5 μs ± 38.8 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

In [14]: %timeit Fibonacci_phi_poly(1023)
19 μs ± 84.3 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

It is only slightly less efficient than fast doubling and much more efficient than simple iteration, and tremendously more efficient than the binomial expansion approach.