Faster Fibonacci Doubling?

Binet's Formula

In chapter 8 of his excellent book Love & Math, Edward Frenkel introduces the concept of a generating function for Fibonacci numbers as an alternative to the usual recursive rule \[F_n = F_{n - 1} + F_{n - 2},n > 2.\] In the notes on page 260, he derives the closed formula which has become known as Binet's formula \[F_n = \frac{1}{\sqrt{5}}\left(\left(\frac{1 + \sqrt{5}}{2}\right)^n - \left(\frac{1 - \sqrt{5}}{2}\right)^n\right).\] To code this formula in Python, for example, one can't use sqrt(5) as this gives an approximation to \(\sqrt{5}\) when an exact result is required. If \(\left(\frac{1 + \sqrt{5}}{2}\right)^n\) is represented as \(\frac{a + b\sqrt{5}}{2}\), \(\left(\frac{1 - \sqrt{5}}{2}\right)^n\) will be \(\frac{a - b\sqrt{5}}{2}\), so their difference is \(b\sqrt{5}\) and \(F_n = b\).

Such a representation can be squared by \[\left(\frac{a + b\sqrt{5}}{2}\right)^2 = \frac{\frac{a^2 + 5b^2}{2} + ab\sqrt{5}}{2}\] The power of this representation can be incremented by \[\left(\frac{a + b\sqrt{5}}{2}\right)\left(\frac{1 + \sqrt{5}}{2}\right) = \frac{\frac{a + 5b}{2} + \frac{a + b}{2}\sqrt{5}}{2}\] Python code using squaring and incrementing with this representation is 

def Binet(n):
    a, b = 2, 0
    bits = bin(n)[2:]
    for bit in bits:  # Iterate over binary representation of n
        a, b = (a * a + 5 * b * b) // 2, a * b # square the representation
        if bit == '1':  # If the current binary bit of n is 1
            a, b = (a + 5 * b) // 2, (a + b) // 2 # increment the exponent
    return b
This code works, but is not as fast as the Fibonacci Fast Doubling method.

Fibonacci Fast Doubling

It can be shown that \[F_{2k} = F_k(2F_{k + 1}- F_k)\] \[F_{2k + 1} = F_{k + 1}^2 + F_k^2\] Python code using these identities is 

def fast_doubling1(n):
    a = 0 # a is F subscript k (k is zero initially)
    b = 1 # b is F subscript k + 1
    bits = bin(n)[2:]
    for bit in bits:  # Iterate over binary representation of n
        a, b = a * (2 * b - a), b * b + a * a # double k
        if bit == '1':  # If the current binary bit of n is 1
            a, b = b, a + b # increment k
    return a

It can also be shown that \[F_{2k} = F_k(2F_{k - 1} + F_k)\] \[F_{2k - 1} = F_{k - 1}^2 + F_k^2\] Python code using these identities is

def fast_doubling2(n):
    a = 0 # a is F subscript k (k is zero initially)
    b = 1 # b is F subscript k - 1
    bits = bin(n)[2:]
    for bit in bits:  # Iterate over binary representation of n
        a, b = a * (2 * b + a), b * b + a * a # double k
        if bit == '1':  # If the current binary bit of n is 1
            a, b = a + b, a # increment k
    return a

It can be assumed that the timing of these functions is due almost entirely to the multiplication of numbers which can become very large. The only difference in timing between fast_doubling1(n) and fast_doubling2(n) is the multiplication b * b, as b is \(F_{2k + 1}\) in fast_doubling1(n) and \(F_{2k - 1}\) in fast_doubling2(n). This does not necessarily imply, however, that fast_doubling2(n) will always run faster than fast_doubling1(n). The time taken to multiply two numbers depends not only on their size but on the pattern of bits in their binary representation. So it's not possible to predict which function would run faster for a given Fibonacci number, and in any case there is not much difference beteen them. Both functions calculate a superfluous value in the final iteration which is not returned, and this is easier to avoid in fast_doubling1(n). It is also highly desirable to use gmpy2 as it is much more efficient for multiple-precision arithmetic than native Python. The recommended Python code is then

def fast_doubling3(n):
    a = gmpy2.mpz(0) # a is F subscript k (k is zero initially)
    b = gmpy2.mpz(1) # b is F subscript k + 1
    bits = bin(n)[2:]
    for bit in bits[:-1]:  # Iterate over binary representation of n
        a, b = a * (2 * b - a), b * b + a * a # double k
        if bit == '1':  # If the current binary bit of n is 1
            a, b = b, a + b # increment k
    return b * b + a * a if bits[-1] == '1' else a * (2 * b - a)

The biggest Fibonacci number that this function can calculate is fast_doubling3(6186179585) which has 1,292,835,074 decimal digits and runs in about 93 seconds on my Windows 10 PC. For higher values, gmpy2 crashes, presumably due to memory management issues when multiplying the huge integers. Higher values can be calculated using the DecInt module by Case Van Horsen instead of gmpy2.

def fast_doubling4(n):
    a = DecInt.DecInt(0) # a is F subscript k (k is zero initially)
    b = DecInt.DecInt(1) # b is F subscript k + 1
    bits = bin(n)[2:]
    for bit in bits[:-1]:  # Iterate over binary representation of n
        a, b = a * (2 * b - a), b * b + a * a # double k
        if bit == '1':  # If the current binary bit of n is 1
            a, b = b, a + b # increment k
    return b * b + a * a if bits[-1] == '1' else a * (2 * b - a)

It's a bit slower, but can give results for much higher values of n. fast_doubling4(10 ** 10) has 2,089,876,403 decimal digits and runs in about 25 minutes. The numbers of digits for powers of 10 agree with OEIS A068070. fast_doubling4(20069628979) has 4,194,304,401 decimal digits and runs in about 3 hours. This is the maximum value of n before it fails on my Windows 10 PC with a memory fault. On my Windows 11 PC, however, it can run without failure for higher values: fast_doubling4(7 * 10 ** 10) has 14,629,134,818 decimal digits and runs in about 15 hours.

How can the value of a very large Fibonacci number be communicated? Its decimal value can be written to a text file in blocks of 100 digits separated with commas and line feeds, like this for \(F_{1000}\) 

                                                                                           434665576,
8693745643568852767504062580256466051737178040248172908953655541794905189040387984007925516929592259,
3080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875

Very large text files can be viewed using apps like the Large Text Viewer.

How can we be sure that the DecInt module has calculated a Fibonacci number correctly? A number N is a Fibonacci number if and only if \(5N^2 + 4\) or \(5N^2 - 4\) is a perfect square (or both). If these expressions are not too big to be held in a gmpy2.mpz integer, they can be tested with (5 * N ** 2 + 4).is_square() and (5 * N ** 2 - 4).is_square(). This test showed that N = fast_doubling4(10 ** 9) was indeed a Fibonacci number but N = fast_doubling4(10 ** 10) was too big to test. For n < 2 ** 32, however, it is easier to just compare the result with gmpy2.fib(n).

Christopher B. Jones 2025-01-20