Matrix Difference Equation for Fibonacci Sequence

Recurrence relation

Suppose we have a difference equation defined by:

Then,

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

Since can transit the state from to , the state vector can be expanded by adding any pair to the above equation.

For example,

can be written into

In fact, this can be generalized. Given by:

, where is constant and 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 ( in above example) and it could be written in:

Matrix Form

If we expand the by taking in above example, then

Computing Fibonacci number by exponentiation

By the above formula, the Fibonacci number can be calculated in . 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 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);
  return F.Read(0, 1);
}
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);
  return F.Read(0, 1);
}

Fast Doubling

If we calculate directly, we can get the equation follows:

Thus, can be calculated by since

This two equations help us to calculate in time, since can be derived from or .

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!