[Python] Is it possible to use the closed-form of Fibonacci series to generate the Nth...

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)n, using the binomial theorem we get 1 + na + nC2a2 + nC3a3 + nC4a4 ... + an. Because the choose function is symmetric we only need to compute half of the terms.

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

Then we expand (1 - a)n: 1 - na + nC2a2 - nC3a3 + nC4a4 - nC5a5 + nC6a6...

The signs alternate, now if we evaluate (1 + a)n - (1 - a)n, 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 = 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 = 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 = power
power *= 5

fib = 0
a, b, c, d = n, 1, 2, -1
i = di = length - 1
for _ in range(length):
fib += c * powers
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?

Continue reading...