Skip to content

Instantly share code, notes, and snippets.

@ChunMinChang
Last active March 6, 2022 13:23
Show Gist options
  • Save ChunMinChang/f80ef4decca23b88df16f2f7846049b6 to your computer and use it in GitHub Desktop.
Save ChunMinChang/f80ef4decca23b88df16f2f7846049b6 to your computer and use it in GitHub Desktop.
Implementation of different methods to calculate the _Fibonacci_ numbers by fast doubling #recursion #dynamic_programming #Fibonacci #math
#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;
}
#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
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)
#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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment