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).
