I’ve written before about computing big Fibonacci numbers, and about creating a certificate to verify a Fibonacci number has been calculated correctly. This post will revisit both, giving a different approach to computing big Fibonacci numbers that produces a certificate along the way.
As I’ve said before, I’m not aware of any practical reason to compute large Fibonacci numbers. However, the process illustrates techniques that are practical for other problems.
The fastest way to compute the nth Fibonacci number for sufficiently large n is Binet’s formula:
Fn = round( φn / √5 )
where φ is the golden ratio. The point where n is “sufficiently large” depends on your hardware and software, but in my experience I found the crossover to be somewhere 1,000 and 10,000.
The problem with Binet’s formula is that it requires extended precision floating point math. You need extra guard digits to make sure the integer part of your result is entirely correct. How many guard digits you’ll need isn’t obvious a priori. This post gave a way of detecting errors, and it could be turned into a method for correcting errors.
But how do we know an error didn’t slip by undetected? This question brings us back to the idea of certifying a Fibonacci number. A number f Fibonacci number if and only if one of 5f2 ± 4 is a perfect square.
Binet’s formula, implemented in finite precision, takes a positive integer n and gives us a number f that is approximately the nth Fibonacci number. Even in low-precision arithmetic, the relative error in the computation is small. And because the ratio of consecutive Fibonacci numbers is approximately φ, an approximation to Fn is far from the Fibonacci numbers Fn − 1 and Fn + 1. So the closest Fibonacci number to an approximation of Fn is exactly Fn.
Now if f is approximately Fn, then 5f2 is approximately a square. Find the integer r minimizing |5f2 − r²|. In Python you could do this with the isqrt function. Then either r² + 4 or r² − 4 is divisible by 5. You can know which one by looking at r² mod 5. Whichever one is is, divide it by 5 and you have the square of Fn. You can take the square root exactly, in integer arithmetic, and you have Fn. Furthermore, the number r that you computed along the way is the certificate of the calculation of Fn.
With the algorithm described above, we don’t need the initial approximation of Fn to be very accurate. We could use ordinary floating point arithmetic, as long as our numbers don’t overflow. If we need extended floating point, it would not be for accuracy but to avoid overflow. It in terms of floating point representation, you need lots of exponent bits but not many significand bits.
You say “The fastest way to compute the nth Fibonacci number for sufficiently large n is Binet’s formula” but I think this is simply not true.
The naive a,b -> b,a+b iteration is of course much slower than Binet, but those are not at all the only options on the table.
Here are some concrete timings. For n=10^3,4,5,6,7 using Binet’s formula (I did the exact same as you, but with fewer guard digits which may make it slightly faster) takes 13us, 65us, 1.8ms, 41ms, 709ms on my machine.
For the same n, using the most simple-minded possible version of the recursion I mentioned in my comment on your previous post takes 3us, 22us, 800us, 32ms, 1.3s on my machine. Aha! that’s slower for n=10^7.
Yeah, but mpmath uses better underlying algorithms for large-integer arithmetic than Python’s integer type does. Specifically, it uses GMP. So I made my most-simple-minded-possible recursive code use GMP too via the same Python library that mpmath uses. 3us, 8us, 130us, 2.9ms, 42ms.
That’s _much_ faster, for all these choices of n, and the ratio is _increasing_ as n increases.