## Recurrence relation

Suppose we have a difference equation defined by:

Then,

could be written into a matrix form(matrix difference equation):

Since $S$ can transit the state from $\vec{x_{n-1}}$ to $\vec{x_n}$, the state vector can be expanded by adding any pair $\vec{x_{t-1}}, \vec{x_t}$ to the above equation.

For example,

can be written into

In fact, this can be generalized. Given $y_n$ by:

, where $c_k$ is constant and $k \in [0, n-1]$ is a integer.

then we can rewritten the equations into:

and we could expand the matrix to

## Fibonacci Sequence

Fibonacci number is defined by:

Obviously, Fibonacci sequence is a difference equation ($a = b = 1$ in above example) and it could be written in:

### Matrix Form

If we expand the $\vec{F_n}$ by taking $t = n + 1$ in above example, then

#### Computing Fibonacci number by exponentiation

By the above formula, the Fibonacci number $F_n$ can be calculated in $O(\log n)$. The key is to compute the exponentiation by squaring.

I explained how to do it in my previous post. Please read it if you need.

As a result, we can compute Fibonacci number $F_n$ as follows:

##### By std::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 x (k^2) ^ ((k-1)/2)  , if n is odd
Matrix pow(unsigned int n)
{
Matrix x(*this); // Copy constructor = Matrix x(rows, cols, data);
Matrix r = Identity(rows);
while (n) {
if (/*n % 2*/n & 1) {
r = r * x;
}
x = x * x;
/*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;
};

// 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)
{
Matrix F { 2, 2, {
{ 1, 1 },
{ 1, 0 }
} };

// 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);
}

##### 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 x (k^2) ^ ((k-1)/2)  , if n is odd
Matrix pow(unsigned int n)
{
Matrix x(*this); // Copy constructor = Matrix x(rows, cols, data);
Matrix r = Identity(rows);
while (n) {
if (/*n % 2*/n & 1) {
r = r * x;
}
x = x * x;
/*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;
};

// 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)
{
Matrix F { 2, 2, new uint64_t*[2] {
new uint64_t[2] { 1, 1 },
new uint64_t[3] { 1, 0 }
} };

// 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);
}


### Fast Doubling

If we calculate $F_{2n}$ directly, we can get the equation follows:

Thus, $F_{2n}, F_{2n + 1}$ can be calculated by $F_n, F_{n+1}$ since

This two equations help us to calculate $F_{N}$ in $O(\log N)$ time, since $F_{N}$ can be derived from $F_{2N'}$ or $F_{2N' + 1}$.

#### implementation

///////////////////////////////////////////////////////////////////////////////
// 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 fib(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 fib(k + 1) * fib(k + 1) + fib(k) * fib(k);
} else { // By F(n) = F(2k) = F(k) * [ 2 * F(k+1) – F(k) ]
k = n / 2;
return fib(k) * (2 * fib(k + 1) - fib(k));
}
}


The implementation are easy to understand but it still has a lot of room to improve. We will discuss it in next post. Stay tuned!