Skip to content

Instantly share code, notes, and snippets.

@ChunMinChang
Last active January 31, 2018 14:16
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save ChunMinChang/b6325c148e8aff15b6e72dcac0aa904e to your computer and use it in GitHub Desktop.
Save ChunMinChang/b6325c148e8aff15b6e72dcac0aa904e to your computer and use it in GitHub Desktop.
Fibonacci numbers #Fibonacci #math #recursion #dynamic_programming

Fibonacci Tricks

Comparison of different calculations for Fibonacci numbers.

TODO

  • Saving all the pair of function name and function pointer and using it to check and print its results
    • e.g. Save a pair like ["Recursive", fibonacci_recursive]

Reference

#include "fibonacci.h"
#include <cassert> // assert
#include <cmath> // pow, sqrt
#include <cstdlib> // calloc, free
// #include <iostream> // std::cout, std::cin, std::endl
#include <utility> // std::swap
#include <vector> // std::vector
///////////////////////////////////////////////////////////////////////////////
// Recursive: O(2^n)
uint64_t fibonacci_recursive(unsigned int n)
{
// if (n <= 1) {
// return n; // F(0) = 0, F(1) = 1.
// }
// return fibonacci_recursive(n-1) + fibonacci_recursive(n-2);
return (n <= 1) ? n : fibonacci_recursive(n-1) + fibonacci_recursive(n-2);
}
///////////////////////////////////////////////////////////////////////////////
// Recursive with memoization: O(n)
std::vector<uint64_t> mem = { 0, 1 }; // F(k) = mem[k], F(0) = 0, F(1) = 1.
uint64_t fibonacci_memoization(unsigned int n)
{
// if (n + 1 <= mem.size()) { // Min size of mem[n] is n + 1
// return mem[n];
// }
// mem.push_back(fibonacci_memoization(n-1) + fibonacci_memoization(n-2));
// return mem[n];
if (n + 1 > mem.size()) { // If n is not calculated yet
mem.push_back(fibonacci_memoization(n-1) + fibonacci_memoization(n-2));
}
return mem[n];
}
///////////////////////////////////////////////////////////////////////////////
// Dynamic programming: O(n)
// uint64_t fibonacci_dynamic(unsigned int n)
// {
// uint64_t a = 0, b = 1, sum = 0; // a = F(0), b = F(1)
// for (unsigned int i = 1 ; i < n ; ++i) { // run if n >= 2
// sum = a + b; // sum = F(i+1)
// a = b; // a = F(i)
// b = sum; // b = F(i+1)
// }
// // Now, i = n, sum = F(n), a = F(n-1), b = F(n)
// return (n < 2) ? n : sum;
// }
uint64_t fibonacci_dynamic(unsigned int n)
{
uint64_t a = 0, b = 1; // a = F(k), b = F(k+1), k = 0 now.
for (unsigned int k = 1 ; k <= n ; ++k) { // loop k from 1 to n.
std::swap(a, b); // a = F(k+1), b = F(k)
b += a; // b = F(k) + F(k+1) = F(k+2)
}
return a;
}
///////////////////////////////////////////////////////////////////////////////
// closed-form: O(log(n))
// Theoretically, the power of n could be done in O(log(n)), but it's
// complicated to calculate the floating numbers.
uint64_t fibonacci_formula(unsigned int n)
{
// double sqrt5 = sqrt((double)5);
double sqrt5 = 2.2360679775;
return (pow((1 + sqrt5) / 2, n) - pow((1 - sqrt5) / 2, n)) / sqrt5;
}
///////////////////////////////////////////////////////////////////////////////
// Rounding: O(log(n))
uint64_t fibonacci_rounding(unsigned int n)
{
// double sqrt5 = sqrt((double)5);
double sqrt5 = 2.2360679775;
double phi = (1 + sqrt5) / 2;
return round(pow(phi, n) / sqrt5);
}
///////////////////////////////////////////////////////////////////////////////
// Power by matrix exponentiation: O(log(n))
// Matrix A:
// <--- cols: n --->
// +- -+
// | A11, A12, ... A1n | ^
// | A21, A22, ... A2n | |
// | ... | rows: m
// | ... | |
// | Am1, Am2, ... Amn | v
// +- -+
// #define MATRIX_BY_VECTOR // Comment this to run native-array version.
///////////////////////////////////////
// By std::vector
#if defined(MATRIX_BY_VECTOR)
class Matrix
{
public:
Matrix(unsigned int r, unsigned int c,
std::vector<std::vector<uint64_t>> d)
: rows(r)
, cols(c)
, data(d)
{
}
Matrix(unsigned int r, unsigned int c)
: rows(r)
, cols(c)
{
assert(rows && cols);
data.resize(rows);
for (unsigned int i = 0 ; i < rows ; ++i) {
data[i].resize(cols);
}
}
~Matrix()
{
}
uint64_t Read(unsigned int r, unsigned int c)
{
return data[r][c];
}
// friend std::ostream& operator<<(std::ostream& os, const Matrix& m)
// {
// for (unsigned int i = 0; i < m.rows; ++i) {
// for (unsigned int j = 0; j < m.cols; ++j) {
// os << m.data[i][j] << " ";
// }
// os << std::endl;
// }
// return os;
// }
Matrix operator*(const Matrix& other)
{
assert(cols == other.rows); // Check if they can be multiplied.
Matrix z(rows, other.cols);
for (unsigned int i = 0 ; i < rows ; ++i) {
for (unsigned int j = 0 ; j < other.cols; ++j) {
for (unsigned int k = 0 ; k < cols; ++k) {
z.data[i][j] += data[i][k] * other.data[k][j];
}
}
}
return z;
}
// Calculate the power by fast doubling:
// k ^ n = (k^2) ^ (n/2) , if n is even
// or k * (k^2) ^ ((k-1)/2) , if n is odd
Matrix pow(unsigned int n)
{
Matrix k(*this); // Copy constructor = Matrix x(rows, cols, data);
Matrix r = Identity(rows);
while (n) {
if (/*n % 2*/n & 1) {
r = r * k;
}
k = k * k;
/*n /= 2*/n >>= 1;
}
return r;
}
private:
Matrix Identity(unsigned int size)
{
Matrix z(size, size);
for (unsigned int i = 0 ; i < size ; ++i) {
z.data[i][i] = 1;
}
return z;
}
unsigned int rows;
unsigned int cols;
std::vector<std::vector<uint64_t>> data;
};
#else
///////////////////////////////////////
// By native array
class Matrix
{
public:
Matrix(unsigned int r, unsigned int c, uint64_t** d = nullptr)
: rows(r)
, cols(c)
, data(d)
{
if (!data) {
AllocateData();
}
}
Matrix(const Matrix& other) // copy ctor
: rows(other.rows)
, cols(other.cols)
{
assert(!data);
AllocateData();
for (unsigned int i = 0 ; i < rows ; ++i) {
for (unsigned int j = 0 ; j < cols ; ++j) {
data[i][j] = other.data[i][j];
}
}
}
~Matrix()
{
FreeData();
}
uint64_t Read(unsigned int r, unsigned int c)
{
return data[r][c];
}
// friend std::ostream& operator<<(std::ostream& os, const Matrix& m)
// {
// for (unsigned int i = 0; i < m.rows; ++i) {
// for (unsigned int j = 0; j < m.cols; ++j) {
// os << m.data[i][j] << " ";
// }
// os << std::endl;
// }
// return os;
// }
Matrix& operator=(Matrix&& other) noexcept // move assignment
{
if(this != &other) {
FreeData(); // Free this data if it exists.
// Move original other.data to data and set other.data to nullptr.
// data = std::exchange(other.data, nullptr); // C++14
data = other.data;
other.data = nullptr;
}
return *this;
}
Matrix operator*(const Matrix& other)
{
assert(cols == other.rows); // Check if they can be multiplied.
Matrix z(rows, other.cols);
for (unsigned int i = 0 ; i < rows ; ++i) {
for (unsigned int j = 0 ; j < other.cols; ++j) {
for (unsigned int k = 0 ; k < cols; ++k) {
z.data[i][j] += data[i][k] * other.data[k][j];
}
}
}
return z;
}
// Calculate the power by fast doubling:
// k ^ n = (k^2) ^ (n/2) , if n is even
// or k * (k^2) ^ ((k-1)/2) , if n is odd
Matrix pow(unsigned int n)
{
Matrix k(*this); // Copy constructor = Matrix x(rows, cols, data);
Matrix r = Identity(rows);
while (n) {
if (/*n % 2*/n & 1) {
r = r * k;
}
k = k * k;
/*n /= 2*/n >>= 1;
}
return r;
}
private:
Matrix Identity(unsigned int size)
{
Matrix z(size, size);
for (unsigned int i = 0 ; i < size ; ++i) {
z.data[i][i] = 1;
}
return z;
}
void AllocateData()
{
assert(!data);
data = (uint64_t**) calloc(rows, sizeof(uint64_t*));
for (unsigned int i = 0 ; i < rows ; ++i) {
data[i] = (uint64_t*) calloc(cols, sizeof(uint64_t));
}
}
void FreeData()
{
if (!data) {
return;
}
assert(rows);
for (unsigned int i = 0 ; i < rows ; ++i) {
free(data[i]);
}
free(data);
data = nullptr;
}
unsigned int rows;
unsigned int cols;
uint64_t** data = nullptr;
};
#endif
// The Fibonacci matrix can be written into the following equation:
// +- -+ +- -+^n
// | F(n+1) F(n) | | 1 1 |
// | | = | |
// | F(n) F(n-1) | | 1 0 |
// +- -+ +- -+
uint64_t fibonacci_matrix(unsigned int n)
{
#if defined(MATRIX_BY_VECTOR)
Matrix F { 2, 2, {
{ 1, 1 },
{ 1, 0 }
} };
#else
Matrix F { 2, 2, new uint64_t*[2] {
new uint64_t[2] { 1, 1 },
new uint64_t[3] { 1, 0 }
} };
#endif
// Using F.data[0][1] since n might be 0.
// (we need to power by n - 1 if we return F.data[0][0].)
F = F.pow(n);
return F.Read(0, 1);
}
///////////////////////////////////////////////////////////////////////////////
// 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 fibonacci_doubling(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 fibonacci_recursive(unsigned int n);
uint64_t fibonacci_memoization(unsigned int n);
uint64_t fibonacci_dynamic(unsigned int n);
uint64_t fibonacci_formula(unsigned int n);
uint64_t fibonacci_rounding(unsigned int n);
uint64_t fibonacci_matrix(unsigned int n);
uint64_t fibonacci_doubling(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 <cstdlib> // calloc, free
#include <cassert> // assert
#include <iostream> // std::cout, std::cin, std::endl
#include <utility> // std::exchange
#include <vector> // std::vector
// Matrix A:
// <--- cols: n --->
// +- -+
// | A11, A12, ... A1n | ^
// | A21, A22, ... A2n | |
// | ... | rows: m
// | ... | |
// | Am1, Am2, ... Amn | v
// +- -+
class Matrix
{
public:
Matrix(unsigned int r, unsigned int c,
std::vector<std::vector<uint64_t>> d)
: rows(r)
, cols(c)
, data(d)
{
}
Matrix(unsigned int r, unsigned int c)
: rows(r)
, cols(c)
{
assert(rows && cols);
data.resize(rows);
for (unsigned int i = 0 ; i < rows ; ++i) {
data[i].resize(cols);
}
}
~Matrix()
{
}
uint64_t Read(unsigned int r, unsigned int c)
{
return data[r][c];
}
friend std::ostream& operator<<(std::ostream& os, const Matrix& m)
{
for (unsigned int i = 0; i < m.rows; ++i) {
for (unsigned int j = 0; j < m.cols; ++j) {
os << m.data[i][j] << " ";
}
os << std::endl;
}
return os;
}
Matrix operator*(const Matrix& other)
{
assert(cols == other.rows); // Check if they can be multiplied.
Matrix z(rows, other.cols);
for (unsigned int i = 0 ; i < rows ; ++i) {
for (unsigned int j = 0 ; j < other.cols; ++j) {
for (unsigned int k = 0 ; k < cols; ++k) {
z.data[i][j] += data[i][k] * other.data[k][j];
}
}
}
return z;
}
// Calculate the power by fast doubling:
// k ^ n = (k^2) ^ (n/2) , if n is even
// or k * (k^2) ^ ((k-1)/2) , if n is odd
Matrix pow(unsigned int n)
{
Matrix k(*this); // Copy constructor = Matrix x(rows, cols, data);
Matrix r = Identity(rows);
while (n) {
if (/*n % 2*/n & 1) {
r = r * k;
}
k = k * k;
/*n /= 2*/n >>= 1;
}
return r;
}
private:
Matrix Identity(unsigned int size)
{
Matrix z(size, size);
for (unsigned int i = 0 ; i < size ; ++i) {
z.data[i][i] = 1;
}
return z;
}
unsigned int rows;
unsigned int cols;
std::vector<std::vector<uint64_t>> data;
};
// class Matrix
// {
// public:
// Matrix(unsigned int r, unsigned int c, uint64_t** d = nullptr)
// : rows(r)
// , cols(c)
// , data(d)
// {
// if (!data) {
// AllocateData();
// }
// }
// Matrix(const Matrix& other) // copy ctor
// : rows(other.rows)
// , cols(other.cols)
// {
// assert(!data);
// AllocateData();
// for (unsigned int i = 0 ; i < rows ; ++i) {
// for (unsigned int j = 0 ; j < cols ; ++j) {
// data[i][j] = other.data[i][j];
// }
// }
// }
// ~Matrix()
// {
// FreeData();
// }
// uint64_t Read(unsigned int r, unsigned int c)
// {
// return data[r][c];
// }
// friend std::ostream& operator<<(std::ostream& os, const Matrix& m)
// {
// for (unsigned int i = 0; i < m.rows; ++i) {
// for (unsigned int j = 0; j < m.cols; ++j) {
// os << m.data[i][j] << " ";
// }
// os << std::endl;
// }
// return os;
// }
// Matrix& operator=(Matrix&& other) noexcept // move assignment
// {
// if(this != &other) {
// FreeData(); // Free this data if it exists.
// // Move original other.data to data and set other.data to nullptr.
// // data = std::exchange(other.data, nullptr); // C++14
// data = other.data;
// other.data = nullptr;
// }
// return *this;
// }
// Matrix operator*(const Matrix& other)
// {
// assert(cols == other.rows); // Check if they can be multiplied.
// Matrix z(rows, other.cols);
// for (unsigned int i = 0 ; i < rows ; ++i) {
// for (unsigned int j = 0 ; j < other.cols; ++j) {
// for (unsigned int k = 0 ; k < cols; ++k) {
// z.data[i][j] += data[i][k] * other.data[k][j];
// }
// }
// }
// return z;
// }
// // Calculate the power by fast doubling:
// // k ^ n = (k^2) ^ (n/2) , if n is even
// // or k * (k^2) ^ ((k-1)/2) , if n is odd
// Matrix pow(unsigned int n)
// {
// Matrix k(*this); // Copy constructor = Matrix x(rows, cols, data);
// Matrix r = Identity(rows);
// while (n) {
// if (/*n % 2*/n & 1) {
// r = r * k;
// }
// k = k * k;
// /*n /= 2*/n >>= 1;
// }
// return r;
// }
// private:
// Matrix Identity(unsigned int size)
// {
// Matrix z(size, size);
// for (unsigned int i = 0 ; i < size ; ++i) {
// z.data[i][i] = 1;
// }
// return z;
// }
// void AllocateData()
// {
// assert(!data);
// data = (uint64_t**) calloc(rows, sizeof(uint64_t*));
// for (unsigned int i = 0 ; i < rows ; ++i) {
// data[i] = (uint64_t*) calloc(cols, sizeof(uint64_t));
// }
// }
// void FreeData()
// {
// if (!data) {
// return;
// }
// assert(rows);
// for (unsigned int i = 0 ; i < rows ; ++i) {
// free(data[i]);
// }
// free(data);
// data = nullptr;
// }
// unsigned int rows;
// unsigned int cols;
// uint64_t** data = nullptr;
// };
int main()
{
Matrix A { 2, 4, {
{ 1, 4, 6, 10 },
{ 2, 7, 5, 3 }
} };
std::cout << A << std::endl;
Matrix B { 4, 3, {
{ 1, 4, 6 },
{ 2, 7, 5 },
{ 9, 0, 11 },
{ 3, 1, 0 },
} };
std::cout << B << std::endl;
Matrix C = A * B;
std::cout << C << std::endl;
Matrix D { 3, 3, {
{ 1, 2, 3 },
{ 4, 5, 6 },
{ 7, 8, 9 }
} };
std::cout << D << std::endl;
Matrix E = D.pow(5);
std::cout << E << std::endl;
Matrix F { 2, 2, {
{ 1, 1 },
{ 1, 0 }
} };
std::cout << F << std::endl;
// Matrix A { 2, 4, new uint64_t*[2] {
// new uint64_t[4] { 1, 4, 6, 10 },
// new uint64_t[4] { 2, 7, 5, 3 }
// } };
// std::cout << A << std::endl;
// Matrix B { 4, 3, new uint64_t*[4] {
// new uint64_t[3] { 1, 4, 6 },
// new uint64_t[3] { 2, 7, 5 },
// new uint64_t[3] { 9, 0, 11 },
// new uint64_t[3] { 3, 1, 0 }
// } };
// std::cout << B << std::endl;
// Matrix C = A * B;
// std::cout << C << std::endl;
// Matrix D { 3, 3, new uint64_t*[3] {
// new uint64_t[3] { 1, 2, 3 },
// new uint64_t[3] { 4, 5, 6 },
// new uint64_t[3] { 7, 8, 9 }
// } };
// std::cout << D << std::endl;
// Matrix E = D.pow(5);
// std::cout << E << std::endl;
// Matrix F { 2, 2, new uint64_t*[2] {
// new uint64_t[2] { 1, 1 },
// new uint64_t[3] { 1, 0 }
// } };
// std::cout << F << std::endl;
Matrix FIB = F.pow(10);
std::cout << FIB << std::endl;
// Fibonacci F(10)
std::cout << FIB.Read(0, 1) << std::endl;
return 0;
}
#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 r = calculate(fibonacci_recursive, k);
Result m = calculate(fibonacci_memoization, k);
Result d = calculate(fibonacci_dynamic, k);
Result f = calculate(fibonacci_formula, k);
Result n = calculate(fibonacci_rounding, k);
Result x = calculate(fibonacci_matrix, k);
Result b = calculate(fibonacci_doubling, k);
assert(r.value == m.value &&
m.value == d.value &&
d.value == f.value &&
f.value == n.value &&
n.value == x.value &&
x.value == b.value);
std::cout << "Fibonacci(" << k << ") = " << r.value << std::endl;
std::cout << "\nTime\n-----" << std::endl;
std::cout << "Recursive \tcpu: " << r.time.cpu << "ms, wall: " << r.time.wall << " ms" << std::endl;
std::cout << "Memoization \tcpu: " << m.time.cpu << "ms, wall: " << m.time.wall << " ms" << std::endl;
std::cout << "Dynamic \tcpu: " << d.time.cpu << "ms, wall: " << d.time.wall << " ms" << std::endl;
std::cout << "Formula \tcpu: " << f.time.cpu << "ms, wall: " << f.time.wall << " ms" << std::endl;
std::cout << "Rounding \tcpu: " << n.time.cpu << "ms, wall: " << n.time.wall << " ms" << std::endl;
std::cout << "Matrix \t\tcpu: " << x.time.cpu << "ms, wall: " << x.time.wall << " ms" << std::endl;
std::cout << "Doubling \tcpu: " << b.time.cpu << "ms, wall: " << b.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