Implementation of different methods to calculate the _Fibonacci_ numbers by fast doubling #recursion #dynamic_programming #Fibonacci #math

README.md

Calculating Fibonacci Numbers by Fast Doubling

Implementation of different methods to calculate the Fibonacci numbers by fast doubling.

Please read my blog post for more detail.

Reference

• Fast Fibonacci Transform
• Fast Doubling Proof for Fibonacci
• Fast Fibonacci
• The Fastest Way to Compute the Nth Fibonacci Number: The Doubling Method

fibonacci.cpp

#include "fibonacci.h"
#include <cassert>  // assert
#include <stack>    // std::stack
#include <vector>   // std::vector

///////////////////////////////////////////////////////////////////////////////
// Fast doubling: O(log(n))
//   Using 2n to the Fibonacci matrix above, we can derive that:
//     F(2n)   = F(n) * [ 2 * F(n+1) – F(n) ]
//     F(2n+1) = F(n+1)^2 + F(n)^2
//     (and F(2n-1) = F(n)^2 + F(n-1)^2)
uint64_t fib0(unsigned int n)
{
  // When n = 2: k = 1 and we want to use F(k+1) to calculate F(2k),
  // However, F(2k) = F(k+1) = F(2) is unknown then.
  if (n < 2) {
    return n; // F(0) = 0, F(1) = 1.
  } else if (n == 2) {
    return 1; // F(2) = 1
  }

  unsigned int k = 0;
  if (n % 2) { // By F(n) = F(2k+1) = F(k+1)^2 + F(k)^2
    k = (n - 1) / 2;
    return fib0(k + 1) * fib0(k + 1) + fib0(k) * fib0(k);
  } else { // By F(n) = F(2k) = F(k) * [ 2 * F(k+1) – F(k) ]
    k = n / 2;
    return fib0(k) * (2 * fib0(k + 1) - fib0(k));
  }
}

uint64_t fib1(unsigned int n)
{
  // When n = 2: k = 1 and we want to use F(k+1) to calculate F(2k),
  // However, F(2k) = F(k+1) = F(2) is unknown then.
  if (n <= 2) {
    return n ? 1 : 0; // F(0) = 0, F(1) = F(2) = 1.
  }

  unsigned int k = n / 2; // k = n/2 if n is even. k = (n-1)/2 if n is odd.
  uint64_t a = fib1(k);
  uint64_t b = fib1(k + 1);

  if (n % 2) { // By F(n) = F(2k+1) = F(k+1)^2 + F(k)^2
    return a * a + b * b;
  }
  // By F(n) = F(2k) = F(k) * [ 2 * F(k+1) – F(k) ]
  return a * (2 * b - a);
}

const unsigned int SIZE = 1000;
// 4 is not a fibonacci number, so using it as initialized value.
const uint64_t INIT = 4;
// In this case, F is not calculated successively. For example,
// To get F(6), we only need F(4), F(3), F(2), F(1), F(0) (no F(5)),
// so the other elements in F is still INIT.
// Another way is to use hash map(std::unordered_map), however,
// it will be slower.
uint64_t MEM[SIZE] = { [0 ... SIZE-1] = INIT };
uint64_t fib2(unsigned int n)
{
  if (MEM[n] != INIT) { // Return the saved answer directly.
    return MEM[n];
  }

  if (n == 0) {
    return (MEM[n] = 0); // F(0) = 0.
  } else if (n <= 2) {
    return (MEM[n] = 1); // F(1) = F(2) = 1.
  }

  unsigned int k = n / 2; // k = n/2 if n is even. k = (n-1)/2 if n is odd.
  uint64_t a = fib2(k);
  uint64_t b = fib2(k + 1);

  // By F(n) = F(2k+1) = F(k+1)^2 + F(k)^2, if n is odd.
  //    F(n) = F(2k) = F(k) * [ 2 * F(k+1) – F(k) ], if n is even.
  return (MEM[n] = (n % 2) ? a * a + b * b : a * (2 * b - a));
}

// void fib3_helper(unsigned int n, uint64_t f[])
// {
//   if (n == 0) {
//     f[0] = 0; f[1] = 1;
//     return;
//   }

//   unsigned int k = 0;
//   if (n % 2) {
//     k = (n - 1) / 2;
//     fib3_helper(k, f);
//     uint64_t a = f[0];            // F(k) = F((n-1)/2)
//     uint64_t b = f[1];            // F(k + 1) = F((n- )/2 + 1)
//     uint64_t c = a * (2 * b - a); // F(n-1) = F(2k) = F(k) * [2 * F(k + 1) - F(k)]
//     uint64_t d = a * a + b * b;   // F(n) = F(2k + 1) = F(k)^2 + F(k+1)^2
//     f[0] = d;                     // F(n)
//     f[1] = c + d;                 // F(n+1) = F(n-1) + F(n)
//   } else {
//     k = n / 2;
//     fib3_helper(k, f);
//     uint64_t a = f[0];            // F(k) = F(n/2)
//     uint64_t b = f[1];            // F(k + 1) = F(n/2 + 1)
//     f[0] = a * (2 * b - a);       // F(n) = F(2k) = F(k) * [2 * F(k + 1) - F(k)]
//     f[1] = a * a + b * b;         // F(n + 1) = F(2k + 1) = F(k)^2 + F(k+1)^2
//   }
// }

// Set f[0], f[1] to F(n), F(n+1).
void fib3_helper(unsigned int n, uint64_t f[])
{
  if (!n) {
    f[0] = 0;
    f[1] = 1;
    return;
  }

  fib3_helper(n / 2, f);
  // k = floor(n/2), so k = n / 2 if n is even, k = (n - 1) / 2 if n is odd.
  uint64_t a = f[0]; // F(k)
  uint64_t b = f[1]; // F(k+1)

  uint64_t c = a * (2 * b - a); // F(2k) = F(k) * [ 2 * F(k+1) – F(k) ]
  uint64_t d = a * a + b * b;   // F(2k+1) = F(k+1)^2 + F(k)^2

  if (n % 2) {    // k = (n - 1) / 2, so F(2k) = F(n-1), F(2k+1) = F(n).
    f[0] = d;     // F(n) = F(2k+1).
    f[1] = c + d; // F(n+1) = F(n-1) + F(n) = F(2k) + F(2k+1).
  } else {        // k = n / 2, so F(2k) = F(n), F(2k+1) = F(n+1).
    f[0] = c;     // F(n) = F(2k).
    f[1] = d;     // F(n+1) = F(2k).
  }
}

uint64_t fib3(unsigned int n)
{
  uint64_t f[2] = { INIT, INIT };
  fib3_helper(n, f);
  return f[0];
}

std::vector<uint64_t> fib4_helper(unsigned int n)
{
  if (!n) {
    // [F(0), F(1)] = [0 , 1]
    return { 0 , 1 };
  }

  std::vector<uint64_t> f(fib4_helper(n / 2));
  // k = floor(n/2), so k = n / 2 if n is even, k = (n - 1) / 2 if n is odd.
  uint64_t a = f[0]; // F(k)
  uint64_t b = f[1]; // F(k+1)

  uint64_t c = a * (2 * b - a); // F(2k) = F(k) * [ 2 * F(k+1) – F(k) ]
  uint64_t d = a * a + b * b;   // F(2k+1) = F(k+1)^2 + F(k)^2

  if (n % 2) { // k = (n - 1) / 2, so F(2k) = F(n-1), F(2k+1) = F(n).
    // [F(n), F(n+1)] = [F(2k+1), F(2k+2)] = [F(2k+1), F(2k) + F(2k+1)]
    return { d, c + d };
  } else { // k = n / 2, so F(2k) = F(n), F(2k+1) = F(n+1).
    // [F(n), F(n+1)] = [F(2k), F(2k+1)].
    return { c, d };
  }
}

uint64_t fib4(unsigned int n)
{
  return fib4_helper(n)[0];
}

// uint64_t fib5(unsigned int n)
// {
//   std::stack<unsigned int> s;
//   while(n) {
//     s.push(n);
//     n /= 2; // n = floor(n/2)
//   }
//   s.push(n); // n = 0 now.

//   uint64_t a; // F(n)
//   uint64_t b; // F(n+1)
//   while (!s.empty()) {
//     unsigned int m = s.top(); s.pop();

//     if (m == 0) {
//       a = 0; // F(0) = 0
//       b = 1; // F(1) = 1
//       continue;
//     }

//     // Let k = floor(m/2), so `a` is F(k) and `b` is F(k+1) now.
//     // k = m/2, if m is even. k = (m-1)/2, if m is odd.
//     uint64_t c = a * (2 * b - a); // F(2k) = F(k) * [ 2 * F(k+1) – F(k) ]
//     uint64_t d = a * a + b * b;   // F(2k+1) = F(k)^2 + F(k+1)^2

//     if (m % 2) {  // m = 2k+1:
//       a = d;      //  F(m) = F(2k+1)
//       b = c + d;  //  F(m+1) = F(m) + F(m-1) = F(2k+1) + F(2k)
//     } else {      // m = 2k:
//       a = c;      //  F(m) = F(2k)
//       b = d;      //  F(m+1) = F(2k+1)
//     }
//   }

//   return a;
// }

uint64_t fib5(unsigned int n)
{
  std::stack<unsigned int> s;
  while (n) {
    s.push(n);
    /*n /= 2*/ n >>= 1;
  }

  uint64_t a = 0; // F(0) = 0
  uint64_t b = 1; // F(1) = 1
  while (!s.empty()) {
    unsigned int m = s.top(); s.pop();

    // Let k = floor(m/2), so `a` is F(k) and `b` is F(k+1) now.
    // k = m/2, if m is even. k = (m-1)/2, if m is odd.
    uint64_t c = a * (2 * b - a); // F(2k) = F(k) * [ 2 * F(k+1) – F(k) ]
    uint64_t d = a * a + b * b;   // F(2k+1) = F(k)^2 + F(k+1)^2

    if (/*m % 2*/ m & 1) {  // m = 2k+1:
      a = d;                //  F(m) = F(2k+1)
      b = c + d;            //  F(m+1) = F(m) + F(m-1) = F(2k+1) + F(2k)
    } else {                // m = 2k:
      a = c;                //  F(m) = F(2k)
      b = d;                //  F(m+1) = F(2k+1)
    }
  }

  return a;
}

uint64_t fib6(unsigned int n)
{
  // The position of the highest bit of n.
  // So we need to loop `h` times to get the answer.
  // Example: n = (Dec)50 = (Bin)00110010, then h = 6.
  //                               ^ 6th bit from right side
  unsigned int h = 0;
  for (unsigned int i = n ; i ; ++h, i >>= 1);

  uint64_t a = 0; // F(0) = 0
  uint64_t b = 1; // F(1) = 1
  for (int j = h - 1 ; j >= 0 ; --j) {
    // n_j = floor(n / 2^j) = n >> j, k = floor(n_j / 2), (n_j = n when j = 0)
    // then a = F(k), b = F(k+1) now.
    uint64_t c = a * (2 * b - a); // F(2k) = F(k) * [ 2 * F(k+1) – F(k) ]
    uint64_t d = a * a + b * b;   // F(2k+1) = F(k)^2 + F(k+1)^2

    if ((n >> j) & 1) { // n_j is odd: k = (n_j-1)/2 => n_j = 2k + 1
      a = d;            //   F(n_j) = F(2k+1)
      b = c + d;        //   F(n_j + 1) = F(2k + 2) = F(2k) + F(2k+1)
    } else {            // n_j is even: k = n_j/2 => n_j = 2k
      a = c;            //   F(n_j) = F(2k)
      b = d;            //   F(n_j + 1) = F(2k + 1)
    }
  }

  return a;
}

// Rewrite the fib6 so the iterator could be `unsigned int` instead of `int`.
uint64_t fib7(unsigned int n)
{
  // The position of the highest bit of n.
  // So we need to loop `h` times to get the answer.
  // Example: n = (Dec)50 = (Bin)00110010, then h = 6.
  //                               ^ 6th bit from right side
  unsigned int h = 0;
  for (unsigned int i = n ; i ; ++h, i >>= 1);

  uint64_t a = 0; // F(0) = 0
  uint64_t b = 1; // F(1) = 1
  for (unsigned int i = 1 ; i <= h ; ++i) {
    // Let j = h-i, n_j = floor(n / 2^j) = n >> j (n_j = n when j = 0, i = h),
    // k = floor(n_j / 2), then a = F(k), b = F(k+1) now.
    uint64_t c = a * (2 * b - a); // F(2k) = F(k) * [ 2 * F(k+1) – F(k) ]
    uint64_t d = a * a + b * b;   // F(2k+1) = F(k)^2 + F(k+1)^2

    if ((n >> (h - i)) & 1) { // n_j is odd: k = (n_j-1)/2 => n_j = 2k + 1
      a = d;                  //   F(n_j) = F(2k+1)
      b = c + d;              //   F(n_j + 1) = F(2k + 2) = F(2k) + F(2k+1)
    } else {                  // n_j is even: k = n_j/2 => n_j = 2k
      a = c;                  //   F(n_j) = F(2k)
      b = d;                  //   F(n_j + 1) = F(2k + 1)
    }
  }

  return a;
}

// uint64_t fib8(unsigned int n)
// {
//   // The position of the highest bit of n.
//   // So we need to loop `h` times to get the answer.
//   // Example: n = (Dec)50 = (Bin)00110010, then h = 6.
//   //                               ^ 6th bit from right side
//   unsigned int h = 0;
//   for (unsigned int i = n ; i ; ++h, i >>= 1);

//   uint64_t a = 0; // F(0) = 0
//   uint64_t b = 1; // F(1) = 1

//   // There is only one `1` in the bits of `mask`. The `1`'s position is same as
//   // the highest bit of n(mask = 2^(h-1) at first), and it will be shifted right
//   // iteratively to do `AND` operation with `n` to check `n / 2^j` is odd
//   // or even.
//   unsigned int mask = 1 << (h - 1);
//   for (int j = h - 1 ; j >= 0 ; --j, mask >>= 1) {
//     // n_j = floor(n / 2^j) = n >> j, k = floor(n_j / 2), (n_j = n when j = 0)
//     // then a = F(k), b = F(k+1) now.
//     uint64_t c = a * (2 * b - a); // F(2k) = F(k) * [ 2 * F(k+1) – F(k) ]
//     uint64_t d = a * a + b * b;   // F(2k+1) = F(k)^2 + F(k+1)^2

//     if (n & mask) { // n_j is odd: k = (n_j-1)/2 => n_j = 2k + 1
//       a = d;        //   F(n_j) = F(2k+1)
//       b = c + d;    //   F(n_j + 1) = F(2k + 2) = F(2k) + F(2k+1)
//     } else {        // n_j is even: k = n_j/2 => n_j = 2k
//       a = c;        //   F(n_j) = F(2k)
//       b = d;        //   F(n_j + 1) = F(2k + 1)
//     }
//   }

//   return a;
// }

// uint64_t fib8(unsigned int n)
// {
//   // The position of the highest bit of n.
//   // So we need to loop `h` times to get the answer.
//   // Example: n = (Dec)50 = (Bin)00110010, then h = 6.
//   //                               ^ 6th bit from right side
//   unsigned int h = 0;
//   for (unsigned int i = n ; i ; ++h, i >>= 1);

//   uint64_t a = 0; // F(0) = 0
//   uint64_t b = 1; // F(1) = 1
//   // There is only one `1` in the bits of `mask`. The `1`'s position is same as
//   // the highest bit of n(mask = 2^(h-1) at first), and it will be shifted right
//   // iteratively to do `AND` operation with `n` to check `n / 2^(h-i)` is odd
//   // or even.
//   for (unsigned int i = 1, mask = 1 << (h - 1) ; i <= h ; ++i, mask >>= 1) {
//     // Let j = h-i, n_j = floor(n / 2^j) = n >> j (n_j = n when j = 0, i = h),
//     // k = floor(n_j / 2), then a = F(k), b = F(k+1) now.
//     uint64_t c = a * (2 * b - a); // F(2k) = F(k) * [ 2 * F(k+1) – F(k) ]
//     uint64_t d = a * a + b * b;   // F(2k+1) = F(k)^2 + F(k+1)^2

//     if (mask & n) { // n_j is odd: k = (n_j-1)/2 => n_j = 2k + 1
//       a = d;        //   F(n_j) = F(2k + 1)
//       b = c + d;    //   F(n_j + 1) = F(2k + 2) = F(2k) + F(2k + 1)
//     } else {        // n_j is even: k = n_j/2 => n_j = 2k
//       a = c;        //   F(n_j) = F(2k)
//       b = d;        //   F(n_j + 1) = F(2k + 1)
//     }
//   }

//   return a;
// }

// uint64_t fib8(unsigned int n)
// {
//   // The position of the highest bit of n.
//   // So we need to loop `h` times to get the answer.
//   // Example: n = (Dec)50 = (Bin)00110010, then h = 6.
//   //                               ^ 6th bit from right side
//   unsigned int h = 0;
//   for (unsigned int i = n ; i ; ++h, i >>= 1);

//   uint64_t a = 0; // F(0) = 0
//   uint64_t b = 1; // F(1) = 1
//   // There is only one `1` in the bits of `mask`. The `1`'s position is same as
//   // the highest bit of n(mask = 2^(h-1) at first), and it will be shifted right
//   // iteratively to do `AND` operation with `n` to check `n / 2^(h-1-i)` is odd
//   // or even.
//   for (unsigned int i = 0, mask = 1 << (h - 1) ; i < h ; ++i, mask >>= 1) {
//     // Let j = h-1-i, n_j = floor(n / 2^j) = n >> j (n_j = n when i = h-1),
//     // k = floor(n_j / 2), then a = F(k), b = F(k+1) now.
//     uint64_t c = a * (2 * b - a); // F(2k) = F(k) * [ 2 * F(k+1) – F(k) ]
//     uint64_t d = a * a + b * b;   // F(2k+1) = F(k)^2 + F(k+1)^2

//     if (mask & n) { // n_j is odd: k = (n_j-1)/2 => n_j = 2k + 1
//       a = d;        //   F(n_j) = F(2k + 1)
//       b = c + d;    //   F(n_j + 1) = F(2k + 2) = F(2k) + F(2k + 1)
//     } else {        // n_j is even: k = n_j/2 => n_j = 2k
//       a = c;        //   F(n_j) = F(2k)
//       b = d;        //   F(n_j + 1) = F(2k + 1)
//     }
//   }

//   return a;
// }

uint64_t fib8(unsigned int n)
{
  // The position of the highest bit of n.
  // So we need to loop `h` times to get the answer.
  // Example: n = (Dec)50 = (Bin)00110010, then h = 6.
  //                               ^ 6th bit from right side
  unsigned int h = 0;
  for (unsigned int i = n ; i ; ++h, i >>= 1);

  uint64_t a = 0; // F(0) = 0
  uint64_t b = 1; // F(1) = 1
  // There is only one `1` in the bits of `mask`. The `1`'s position is same as
  // the highest bit of n(mask = 2^(h-1) at first), and it will be shifted right
  // iteratively to do `AND` operation with `n` to check `n_j` is odd or even,
  // where n_j is defined below.
  for (unsigned int mask = 1 << (h - 1) ; mask ; mask >>= 1) { // Run h times!
    // Let j = h-i (looping from i = 1 to i = h), n_j = floor(n / 2^j) = n >> j
    // (n_j = n when j = 0), k = floor(n_j / 2), then a = F(k), b = F(k+1) now.
    uint64_t c = a * (2 * b - a); // F(2k) = F(k) * [ 2 * F(k+1) – F(k) ]
    uint64_t d = a * a + b * b;   // F(2k+1) = F(k)^2 + F(k+1)^2

    if (mask & n) { // n_j is odd: k = (n_j-1)/2 => n_j = 2k + 1
      a = d;        //   F(n_j) = F(2k + 1)
      b = c + d;    //   F(n_j + 1) = F(2k + 2) = F(2k) + F(2k + 1)
    } else {        // n_j is even: k = n_j/2 => n_j = 2k
      a = c;        //   F(n_j) = F(2k)
      b = d;        //   F(n_j + 1) = F(2k + 1)
    }
  }

  return a;
}

fibonacci.h

#ifndef FIBONACCI
#define FIBONACCI

#include <cstdint> // uint64_t

uint64_t fib0(unsigned int n);
uint64_t fib1(unsigned int n);
uint64_t fib2(unsigned int n);
uint64_t fib3(unsigned int n);
uint64_t fib4(unsigned int n);
uint64_t fib5(unsigned int n);
uint64_t fib6(unsigned int n);
uint64_t fib7(unsigned int n);
uint64_t fib8(unsigned int n);

#endif // FIBONACCI

makefile

CC=g++
CFLAGS=-Wall --std=c++11

SOURCES=test.cpp\
        fibonacci.cpp
OBJECTS=$(SOURCES:.cpp=.o)

# Name of the executable program
EXECUTABLE=run_test

all: $(EXECUTABLE)

# $@ is same as $(EXECUTABLE)
$(EXECUTABLE): $(OBJECTS)
	$(CC) $(CFLAGS) $(OBJECTS) -o $@

# ".cpp.o" is a suffix rule telling make how to turn file.cpp into file.o
# for an arbitrary file.
#
# $< is an automatic variable referencing the source file,
# file.cpp in the case of the suffix rule.
#
# $@ is an automatic variable referencing the target file, file.o.
# file.o in the case of the suffix rule.
#
# Use -c to generatethe .o file
.cpp.o:
	$(CC) -c $(CFLAGS) $< -o $@

clean:
	rm *.o $(EXECUTABLE)

test.cpp

#include "fibonacci.h"
#include <cstdint>  // uint64_t
#include <cassert>  // assert
#include <ctime>    // std::clock
#include <iostream> // std::cout, std::cin, std::endl

typedef struct Time
{
  double cpu;
  double wall;
} Time;

typedef struct Result
{
  uint64_t value;
  Time time;
} Result;

template <typename Function, typename... Args>
Result calculate(Function aFunction, Args&&... aArgs)
{
  std::clock_t cpuStart = std::clock();
  auto wallStart = std::chrono::high_resolution_clock::now();

  uint64_t r = aFunction(std::forward<Args>(aArgs)...);

  std::clock_t cpuEnd = std::clock();
  auto wallEnd = std::chrono::high_resolution_clock::now();
  return Result { r, Time { 1000.0 * (cpuEnd - cpuStart) / CLOCKS_PER_SEC,
                            std::chrono::duration<double, std::milli>
                              (wallEnd - wallStart).count()
                          }
                };
}

int main()
{
  unsigned int k = 0;
  std::cout << "Fibonacci number at: ";
  std::cin >> k;
  Result f0 = calculate(fib0, k);
  Result f1 = calculate(fib1, k);
  Result f2 = calculate(fib2, k);
  Result f3 = calculate(fib3, k);
  Result f4 = calculate(fib4, k);
  Result f5 = calculate(fib5, k);
  Result f6 = calculate(fib6, k);
  Result f7 = calculate(fib7, k);
  Result f8 = calculate(fib8, k);
  assert(f0.value == f1.value &&
         f1.value == f2.value &&
         f2.value == f3.value &&
         f3.value == f4.value &&
         f4.value == f5.value &&
         f5.value == f6.value &&
         f6.value == f7.value &&
         f7.value == f8.value);
  std::cout << "Fibonacci(" << k << ") = " << f0.value << std::endl;

  std::cout << "\nTime\n-----" << std::endl;
  std::cout << "0\tcpu: " << f0.time.cpu << "ms, wall: " << f0.time.wall << " ms" << std::endl;
  std::cout << "1\tcpu: " << f1.time.cpu << "ms, wall: " << f1.time.wall << " ms" << std::endl;
  std::cout << "2\tcpu: " << f2.time.cpu << "ms, wall: " << f2.time.wall << " ms" << std::endl;
  std::cout << "3\tcpu: " << f3.time.cpu << "ms, wall: " << f3.time.wall << " ms" << std::endl;
  std::cout << "4\tcpu: " << f4.time.cpu << "ms, wall: " << f4.time.wall << " ms" << std::endl;
  std::cout << "5\tcpu: " << f5.time.cpu << "ms, wall: " << f5.time.wall << " ms" << std::endl;
  std::cout << "6\tcpu: " << f6.time.cpu << "ms, wall: " << f6.time.wall << " ms" << std::endl;
  std::cout << "7\tcpu: " << f7.time.cpu << "ms, wall: " << f7.time.wall << " ms" << std::endl;
  std::cout << "8\tcpu: " << f8.time.cpu << "ms, wall: " << f8.time.wall << " ms" << std::endl;

  return 0;
}