|
#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; |
|
} |