Yeah, I know, it's been done, it's boring, nothing new to see here. I understand, this is probably all known, but I just wanted to show this.
As I'm sure you are all aware, the Fibonacci sequence can be defined recursively via the equation:
\[ F(n) = F(n-1) + F(n-2)~~~~~~~F(n) = n~~\forall~~n \in \{0, 1\} \]
Now, this is fine and good, but if you translate this formula into a recursive algorithm in the most naive way possible, it is also extremely computationally expensive. The fact that each number requires the knowledge of the two previous numbers means that a naive algorithm to compute this function runs in exponential time. The Lispy way (IMHO) to deal with this exponential explosion is to use memoization as a dynamic programming method. This will cache the previous results and turn this function into a linear time algorithm. Another common way to compute this is to run through the entire sequence until you get the value you are looking for. This also runs in linear time but (usually) with a smaller constant factor and has constant space requirements. Unfortunately, it requires a substantial change from the intuitive description of the problem, making the code much harder to understand, which is why I don't like it as much as the memoized recursion. However, in the name of clarity when exploring the complexity, I will be using the "sum up from the base case method".
But, anyway, my point was that you can easily write a linear time Fibonacci number calculator. That is what we would say under normal circumstances, but this isn't really true, or rather, it is only true if you assume constant time multiplication. For large \(n\), it becomes clear that this algorithm is not, actually, linear. This is because it takes \(O(N)\) time to add two \(N\) digit numbers, and because at least some of the numbers that we wish to add are of order \(F(n)\) (which have \(O(\log F(n)) \sim O(n)\) digits), we see now that the actual time is more like \(O(n\log F(n)) \sim O(n^2)\). This might seem a bit surprising, but look and see for yourself.
A quadratic fit to the time scaling of the Fibonacci sequence generator function
But, if we are willing to throw clarity out the window all together, we can also use a different technique to compute the number. We can use the closed form solution to the Fibonacci sequence, Binet's formula (all of this is from Wikipedia):
\[ F(n) = \frac {\psi^n - \phi^n}{\psi - \phi} \]
…where…
\[ \psi = \frac{1+\sqrt 5}{2} \\ \phi = \frac{1-\sqrt 5}{2} \]
And, because \(\psi = -\phi^{-1}\), we can simply say:
\[ F(n) = \frac {\phi^n - (- \phi)^{-n}}{\sqrt 5} \]
Now, I hear you thinking, these aren't going to work if you actually want the
exact number, because these have nasty things like \(\sqrt 5\) and and \(\phi^{-n}\)
in them. And you would be correct that this is not something that you can
simply type into any programming language and expect a valid number back. But
the fact of the matter is, these formulas are actually correct, it is the
language and the language's insistence that computation on general real numbers
be performed using approximate floating point arithmetic that causes such an
issue in the first place. If, however, we were to turn to the venerable bc,
an arbitrary precision calculator that is already on your computer if you have a
UNIX like operating system, then we can try it out for various \(n\) by setting
the precision to a "high enough" value (the scale variable sets how many
decimal digits of precision you have).
scale=100 phi = (1+sqrt(5))/2 define fib (n) { (phi^n - (-phi)^(-n))/sqrt(5) } fib(5) 4.999999999999999999999999999999999999999999999999999999999999999999\ 9999999999999999999999999999999988 fib(10) 54.99999999999999999999999999999999999999999999999999999999999999999\ 99999999999999999999999999999999728 fib(100) 354224848179261915074.9999999999999999999999999999999999999999999999\ 999999999999999999999999999999981555491928230697552907 fib(200) 280571172992510140037611932413038677189524.9999999999999999999999999\ 99999999999999999999999999999997069407044903016136967720432076973572\ 3021244 fib(400) 17602368064501396646822694539241125077038438330449219188672599289657\ 5345044216019674.999999999999996317359669914941458065726733504366353\ 6782163342769596482054432795143577733011851910134 fib(500) 13942322456169788013972438287040728395007025658769730726410896294832\ 5571622863290691557658876222517646901.483330817850463056986072289669\ 50027087266539483120658105031310816022996182603489404860217143670694\ 09
We see that as the numbers are computed in the most significant digits, the error accumulates in the least significant. Once the error meets the important numbers, which has happened at the 500th Fibonacci number, we have reached the point where we actually don't know if the last digits are correct or not. At that point we need to increase the scale and try again.
Also, if you try this at home, you will find that the computation takes longer for larger values of \(n\), despite this being a closed form expression. This isn't really surprising, the closed for expression may seem magical, but it still needs to be implemented with algorithms that multiply, divide, exponentiate, and "square root" very large or very precise numbers. More on this later.
Incidentally, the CLISP implementation of Common Lisp also has arbitrary
precision arithmetic and can be used to perform the same calculation much faster
(also, you can also do this in other imps with some success using rational
numbers if you take the time to write a sqrt function that will return a
rational approximation of the square root to within the needed accuracy). If
you try this at home using CLISP, note that the scale value in bc designates
decimal digits while (long-float-digits) in CLISP designates binary digits
(i.e. the CLISP one needs to be larger for the same precision. For example,
CLISP needs a value of 3322 to give around 1000 decimal digits of precision).
(setf (long-float-digits) 350) (let ((phi (/ (+ 1 (sqrt 5L0)) 2L0))) (defun fib (n) (/ (- (expt phi n) (expt (- phi) (- n))) (sqrt 5L0)))) (fib 100) 3.542248481792619150750000000000000000000000000000000000000000000000000\ 000000000000000000000000000000000005L20 (fib 200) 2.805711729925101400376119324130386771895250000000000000000000000000000\ 0000000000000000000000000000000000084L41 (fib 400) 1.760236806450139664682269453924112507703843833044921918867259928965753\ 4504421601967500000000000000000000109L83 (fib 500) 1.394232245616978801397243828704072839500702565876973072641089629483255\ 7162286329069155765887622252129412605L104
This works okay, but we have to constantly check to make sure that we haven't run out of precision in our number. We can fix this by computing the number of significant digits needed to accurately compute \(F(n)\) by first computing an approximate \(F(n)\) with floating point numbers and using that approximate value as a guess of how many \(b\)-ary digits you need in have by taking the \(\log_b\).
(defun fib-approx (n) (let ((phi (/ (+ 1 (sqrt 5d0)) 2d0))) (/ (- (expt (float phi 0d0) n) (expt (- (float phi 0d0)) (- n))) (sqrt 5d0)))) (defun fib (n) (setf (long-float-digits) (max 64 (round (* 1.5 (log (fib-approx n) 2))))) (let ((phi (/ (+ 1 (sqrt 5L0)) 2L0))) (nth-value 0 (round (/ (- (expt phi n) (expt (- phi) (- n))) (sqrt 5L0)))))) (fib 5) 5 (fib 100) 354224848179261916032 (fib 1000) 43466557686937456435688527675040625802564660517371780402481729089536\ 55541794905189040387984007925516929592259308032263477520968962323987\ 33224711616429964409065331879382989696499285160037044761377951668492\ 28875
However, this only really works up to the point when the double floats start to overflow (Fibonacci numbers less than around 300 digits). After that our fast approximation no longer works as the Fibonacci numbers are too big to fit into any machine sized number. Is there a way to just get an approximate ballpark figure for the size of a Fibonacci number without actually calculating it? Yes. All we need to do is look at the what is important in the equation when \(n\) is large. If \(n\) is large, \((-\phi)^{-n}\) is extremely small since \(\phi^{-1}\) is smaller than 1. This means that we can approximate the large \(n\) Fibonacci numbers as:
\[ F(n\gg 1) = \frac{\phi^n}{\sqrt 5} \]
Now, we know that the number of \(b\)-ary digits in a number is roughly equal to the \(\log_{b}\) of the number. This means that the number of digits in a large Fibonacci number is:
\[ \log F(n \gg 1) = \log \phi^n - \log \sqrt 5 = n\log\phi - \log\sqrt 5 \]
Adding this to our Fibonacci number calculator:
(defun log-fib-approx (n &optional (base 2)) (let ((phi (/ (+ 1 (sqrt 5d0)) 2d0))) (- (* n (log phi base)) (log (sqrt 5) base)))) (defun fib (n) (setf (long-float-digits) (max 64 (round (* 1.5 (log-fib-approx n 2))))) (let ((phi (/ (+ 1 (sqrt 5L0)) 2L0))) (nth-value 0 (round (/ (- (expt phi n) (expt (- phi) (- n))) (sqrt 5L0)))))) (fib 5) 5 (fib 100) 354224848179261915075 (fib 1000) 43466557686937456435688527675040625802564660517371780402481729089536555417\ 94905189040387984007925516929592259308032263477520968962323987332247116164\ 2996440906533187938298969649928516003704476137795166849228875 (fib 10000) 33644764876431.......84058580903081179711592936853976188104504236323110912
At this point, we might ask how all this compares to calculating the number via the recurrence equation? On the one hand, using intelligent algorithms we have found an \(O(n^2)\) time algorithm based on the recurrence relation. On the other, it isn't totally clear what the scaling behavior of the proposed closed form algorithm is. In the Binet's formula based method, we have (roughly in order carried out in the calculation):
- An arbitrary precision square root
- An fixed to arbitrary precision addition (adding one to a value)
- An arbitrary precision to fixed division (divide by two)
- Two arbitrary precision exponentiations
- An arbitrary precision subtraction
- An arbitrary precision division
Looking at this list, the obvious prime suspects for the bottleneck are the square root and the exponentiation, but let's look at each. Wikipedia tells me that multiplication and division have the same complexity lower bounded with \(O(N\log N)\) (probably) and that addition and subtraction run in \(O(N)\) time.
Using simple divide and conquer exponentiation gives an \(O(N\log N \log n)\) time scaling for exponent \(n\). The only thing that is left is the scaling of the square root. Looking through the GMP docs tells me that this can be done in \(O(N\log N)\) time (again assuming the Schönhage-Strassen conjecture for multiplication). This means that the most expensive part of the calculation is the exponentiation, at \(O(\log F(n)\log \log F(n) \log n)\) time. What does that mean for our algorithm?
The Fibonacci sequence grows in an exponential manner. This is clear from recurrence, \(F(n+1) = F(n) + F(n-1) \sim 2F(n)\). Roughly speaking, increasing \(n\) by a constant factor has the effect of multiplying \(F\) by a constant factor, which is basically the definition of "exponential". But, in our complexity, We are taking the \(\log\) of that number, which "undoes" the exponential function. More concretely, we actually calculated a large \(n\) approximation of \(\log F(n)\) above, which turned out to scale linearly with \(n\). All of this boils down to the relationship I mentioned before: \(O(\log F(n)) \sim O(n)\). Thus, we expect this calculation to be \(O(n (\log n)^2)\), whereas the complexity for the recurrence based results took \(O(n^2)\). This result demonstrates that Binet's formula is, in fact, an asymptotic winner in this computational contest. When you plot the timing data, it seems very linear (I certainly cannot seem to resolve the \((\log n)^2\) portion of the scaling behavior).
If you are actually interested in calculating Fibonacci numbers very quickly and are using CLISP, then it is probably a safe bet to use the Binet's formula method rather than the more intuitive algorithm even for small numbers. On CLISP, the differences in run times are around two orders of magnitude faster than the intuitive method (although it seems to be much less predictable in terms of run time in any particular run, my prime suspect on this is the GC). Using SBCL (a much more sophisticated compiler), the intuitive algorithm runs significantly faster than on CLISP, but still is no match for Binet's method (on CLISP).
