Simplified closed form for Fibonacci numbers and O(1) implementation

Question

A friend of mine got the task to implement Fibonacci such that it will take less than 10 seconds for the 2000000th number. This was an interesting task so I set myself the task to make a super fast implementation for any $n$.

The trivial recursion algorithm will take too much time ($O(F_n)$ operations), and using dynamic programming also won't work ($O(n)$ operations). Here even the closed form solution might fail as it take $O(\log n)$ operations using fast exponentiation. Another problem is that $F_n$ becomes huge and operations become more and more expensive.

To avoid these problems I decided to compute $\log(F_n)$ instead, The intuition is that $\log a^n = n\log a$ so we reduce the number of operations to 1. First let me present the math behind my code: $$a = \frac{1+\sqrt{5}}{2}, b = \frac{1-\sqrt{5}}{2}, c=\frac{1}{\sqrt{5}}, d=\frac{b}{a}$$ $$F_n = c(a^n - b^n) = ca^n(1-d^n)$$ I use tilde to denote numbers in the log domain (e.g $\tilde{a} = \log a)$: $$\tilde{F}_n = \tilde{c} + n\tilde{a} + \log(1-d^n)$$ This almost solves the exponentiation problem, but we still have the $d^n$, which can also be solved by: $$= \tilde{c} + n\tilde{a} + \log(1-(-1)^ne^{n\log (-d)})$$ This is not so clean as we could do the exact same thing in the $F_n$ formula, ideally I want to remove the $1-d^n$ completely - Note this does help numerically.

Now for the more technical part, I implemented this exact algorithm to python using python (numpy):

log_d = log(-(1 - sqrt(5))/(1 + sqrt(5)))
sign = -1 if n % 2 else 1
return log(1 / sqrt(5)) + n*log((1 + sqrt(5)) / 2) + log(1 - sign * exp(n * log_d))

This code works well, with less than 0.0001 seconds for n=2 million, and I haven't found an instance where round(exp(log_fibonacci(n))) != fibonacci(n). Something I noticed about my code is that round(exp(log_fibonacci(n))) != fibonacci(n) is 0 for $n>15$ is 0 because of numerical issues. This is actually very interesting, because it means my function computes $$\tilde{F}_n = \tilde{c} + n\tilde{a}$$ which consists of only 2 operations! This raised the following question: Is there a constant $N$ such that any $n>N$ satisfies $F_n=\text{round}\left(\frac{1}{\sqrt{5}}\left(\frac{1+\sqrt{5}}{2}\right)^n\right)$? If not, can we find the $n$'s that does not satisfy this formula?

Answer 1

It is indeed easy to verify that the rounding formula works, since $b^n$ approaches $0$ very fast.

Numerical Issues

As you have noted, however, there are severe numerical issues with this approach. It is clear from the relationships that you have written that $\log F_n$ is most nearly $n\log a$. Since you've stored this as a floating point number, you are essentially storing the mantissa and exponent simultaneously:

$$\log F_n=\rm\log(mantissa\times2^{exponent})=\underbrace{exponent\log2+\log mantissa}$$

In order to store the exponent with the mantissa, you lose significant digits in the mantissa.

To offset this, one requires increasing precision as $n$ increases. This means one of two things:

1. Either we must restrict the algorithm to small $n$ or

2. We need to use more precision as $n$ increases and find a way to compute the golden ratio further.

As you saw, double precision only works up until $n=15$, and when one considers how to handle larger $n$, all of the additional computations make this more than simply $2$ operations.

Exact Computation

Of course, if a rough approximation is all one desires, then this approach is fine.

For exact results that do not rely on floating point, common methods use the matrix form or the identity

$$\varphi^n=F_n\varphi+F_{n-1}$$

with exponentiation by squaring, or their derived identities

$$F_{2n-1}=F_n^2+F_{n-1}^2\\F_{2n}=(2F_{n-1}+F_n)F_n$$

where one computes pairs of Fibonacci numbers at a time or uses dynamic programming. These examples are all pointed out on Wikipedia.

Here is a python implementation which computes $\varphi^n$ using exponentiation by squaring. It can be seen that these types of identities compute the Fibonacci numbers exactly, and though they require large integer arithmetic, computing the $2$-millionth Fibonacci number takes no more than a second.

Answer 2

You just need $$|cb^n| \lt 0.5\\ \frac 1{\sqrt 5}\left(\frac {\sqrt 5-1}2\right)^n\lt 0.5\\ n \log \left(\frac {\sqrt 5-1}2\right)\lt \log \frac {\sqrt 5}2\\ n \gt \frac{\log \frac {\sqrt 5}2}{\log \left(\frac {\sqrt 5-1}2\right)}$$
Which is already true at $n=0$, so you can replace the $b^n$ term with rounding any time you like.