Return to Answer

n is treated as the sum of powers of 2: 2^a + 2^b + ... for each iteration e, let p = 2^i, then

n is treated as the sum of powers of 2: 2^a + 2^b + ... for each iteration i (starting from 0), let p = 2^i, then

Note: older compilers may be missing <inttypes.h> or <stdint.h>. If <inttypes.h> is not present, use compiler specific format string for uint64_t. If <stdint.h> is not present, use typedef ... uint64_t (usually unsigned long long) and the appropriate format string.

During the calculation of a and b, let c = current cumulative sum of bits of n:

b = F(c)
a = F(c+1)

a' = F(c+1+p) = F(c+2) F(p) + F(c+1) F(p-1)
   = (F(c+1)+F(c)) F(p) + F(c+1) F(p-1)
   = F(c+1) F(p) + F(c) F(p) + F(c) F(p-1)
   = a d + b d + b c

b' = F(c+p) = F(c+1) F(p) + F(c) F(p-1)
   = a d + b c

Note: older compilers may be missing <inttypes.h> or <stdint.h>. If <inttypes.h> is not present (Visual Studio 2010), use compiler specific format string for uint64_t. If <stdint.h> is not present (Visual Studio 2005), use typedef ... uint64_t (usually unsigned long long) and the appropriate format string.

During the calculation of a and b, let m = current cumulative sum of bits of n:

b = F(m)
a = F(m+1)

a' = F(m+1+p) = F(m+2) F(p) + F(m+1) F(p-1)
   = (F(m+1)+F(m)) F(p) + F(m+1) F(p-1)
   = F(m+1) F(p) + F(m) F(p) + F(m) F(p-1)
   = a d + b d + b c

b' = F(m+p) = F(m+1) F(p) + F(m) F(p-1)
   = a d + b c

Note that if b' is the max value for uint64_t, a' will overflow, but it's not an issue. However, the algorithm can be modified so that when completed, a = fib(n-1):

a = fib(-1) = 1
b = fib(0)  = 0

a = fib(m-1)
b = fib(m)

b' = fib(m+p)
   = fib(m+1)fib(p) + fib(m)fib(p-1)
   = (fib(m) + fib(m-1))fib(p) + fib(m)fib(p-1)
   = fib(m)fib(p) + fib(m-1)fib(p) + fib(m)fib(p-1)
   = bd           + ad             + bc

a' = fib(m-1+p)
   = fib(m)fib(p) + fib(m-1)fib(p-1)
   = bd           + ac

uint64_t fib(uint64_t n)
{
    uint64_t a, b, c, d;
    a = d = 1;
    b = c = 0;
    while (1) {
        if (n & 1) {
            uint64_t bd = b*d;
            b = bd + a*d + b*c;
            a = bd + a*c;
        }
        n >>= 1;
        if (n == 0)
            break;
        {
            uint64_t dd = d*d;
            d = dd + 2*d*c;
            c = dd + c*c;
        }
    }
    return b;
}

Example code; c and d are used for the repeated squaring, while a and b are the cumulative results and end up as a = fib(n+1), b = fib(n). Older compilers may not have include files for inttypes.h or stdint.h. If inttypes.h is not present, use compiler specific format string for uint64_t. If stdint.h is not present, use typedef ... uint64_t (usually unsigned long long).

Example code; c and d are used for the repeated squaring, while a and b are the cumulative results and end up as a = fib(n+1), b = fib(n).

Note: older compilers may be missing <inttypes.h> or <stdint.h>. If <inttypes.h> is not present, use compiler specific format string for uint64_t. If <stdint.h> is not present, use typedef ... uint64_t (usually unsigned long long) and the appropriate format string.