Calculating Fibonacci number exactly in C ++?

To accurately calculate Fibonacci numbers using Binet's formula, you need an accurate interpretation of & radic; Since & radic; 5 is irrational, it cannot be accurately represented with double

or float

, so Binet's formula does not work with these types (however, rounding in the calculation produces accurate results for some small inputs). Since Fibonacci numbers are integers, you can get accurate results from Binet's formula using double

or float

for more arguments, rounding afterwards.

double binet(unsigned int n)
{
    static const double phi = (1 + sqrt(5))*0.5;
    double fib = (pow(phi,n) - pow(1-phi,n))/sqrt(5);
    return round(fib);
}

      

        
        
        
      

    

This will return the correct result for almost all n

small enough so that the result can be represented exactly as double

. However, they are few. A double

usually only has 53 bits of precision, so only Fibonacci numbers smaller than 2 ⁵³ can be accurately represented as double

(plus a few larger ones divisible by powers of 2 high enough), the last Fibonacci number less than 2 ⁵³ is F (77), but F (78) is divisible by 8, so it is also exactly represented as a double

with 53 bits of precision. However, the above only gives correct results for n <= 70

here 71 through, the rounding error is too large (by the way, the result of the Binet formula using doubles

here is always too large, so use floor

insteadround

will give the correct result also for F (71), but not further).

With standard data types, not many Fibonacci numbers are accurately represented, the latter to match (unsigned) 64-bit type is F (93); for 128 bits, the last is F (186). For such small indices, there is practically nothing compared to a simple iterative algorithm

unsigned long long fibonacci(unsigned int n)
{
    unsigned long long a = 0, b = 1;
    for(; n > 0; --n)
    {
        b += a;
        a = b-a;
    }
    return a;
}

      

        
        
        
      

    

unless you are using lookup table

static const unsigned long long fibs[94] = { 0, 1, 1, 2, ... , 12200160415121876738ull };

      

        
        
        
      

    

For accurate results, consider & radic; 5 (and / or .phis;) as a symbolic constant and evaluate the formula using that. It boils down to evaluating the formula in the ring

ℤ[φ] = { a + b*φ : a, b ∈ ℤ }

      

        
        
        
      

    

algebraic integers in ℚ(√5)

, using the fact that φ² = 1 + φ

. Equivalent to Binet's formula:

φ^n = F(n-1) + φ*F(n)

      

        
        
        
      

    

which can be used to efficiently compute Fibonacci numbers by re-squaring in O (log n) steps (but note that F (n) has an & Theta; (n) bit, so the number of bit operations cannot be lower than Na )). A slightly more efficient version than vanilla re-squaring uses

φ^(2n) = (φ^n)² = (F(n-1) + φ*F(n))² = F(n-1)² + φ*2*F(n-1)*F(n) + φ²*F(n)²
       = (F(n-1)² + F(n)²) + φ*(2*F(n-1)*F(n) + F(n)²)

      

        
        
        
      

    

find F(2n) = 2*F(n)*F(n-1) + F(n)² = 2*F(n)*F(n+1) - F(n)² = F(n)*(F(n+1) + F(n-1))

and F(2n+1) = F(n)² + F(n+1)²

using φ² = 1 + φ

. These formulas allow calculating F (2n), F (2n + 1) and F (2n + 2) from F (n) and F (n + 1) with at most two multiplications and two additions / subtractions by a number, which gives an algorithm to compute a pair (F(n),F(n+1))

in O (log n) steps with only two numbers as a state (vanilla re-squaring uses four numbers as a state and requires a few more multiplications).

Left to right iterative algorithm -

unsigned long long fib(unsigned int n){
    if (n == 0) return 0;
    unsigned int h = n/2, mask = 1;
    // find highest set bit in n, can be done better
    while(mask <= h) mask <<= 1;
    mask >>= 1;
    unsigned long long a = 1, b = 1, c; // a = F(k), b = F(k+1), k = 1 initially
    while(mask)
    {
        c = a*a+b*b;        // F(2k+1)
        if (n&mask)
        {
            b = b*(b+2*a);  // F(2k+2)
            a = c;          // F(2k+1)
        } else {
            a = a*(2*b-a);  // F(2k)
            b = c;          // F(2k+1)
        }
        mask >>= 1;
    }
    return a;
}

      

        
        
        
      

    

With an arbitrary type of precision instead unsigned long long

, allowing you to quickly calculate large Fibonacci numbers. But, of course, arbitrary-precision libraries often have their own optimized Fibonacci functions, so it might make sense to implement them.