[模板]矩阵快速幂(斐波那契数列)

首先回顾矩阵乘法的定义

\[c_{ij}=\sum_{i=1}^{k}a_{ik} \times b_{kj} \]

显然，对\(F_{n}=F_{n-1}+F_{n-2}\)这样的柿子，我们可以用待定系数法求得递推矩阵：

设\(f_{n}=\begin{bmatrix} F_{n} & F_{n-1} \end{bmatrix}(n > 1)\),解\(f_{n}=f_{n-1} \begin{bmatrix} a & b \\ c & d \end{bmatrix}\)
得

\[\begin{bmatrix} a & b \\ c & d \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \]

于是

\[f_{n} = f_{2} \times\begin{bmatrix} 1 & 1\\ 1 &0 \end{bmatrix}^{n-2} \]

利用矩阵快速幂即可\(O(\log_{}{n})\)求解\(F_{n}\)
我们定义一个结构体来进行矩阵的表示，初始化，乘法运算等

struct Matrix {
  ll a[maxn][maxn];
  ll r, c;

  Matrix() {
    memset(a, 0, sizeof(a));//构造时进行初始化
  }

  void print() const { // for debug
    if (r < 0 || c < 0) return;
    for (int i = 1; i <= r; i++) {
      for (int j = 1; j <= c; j++) {
        printf("%lld ", a[i][j]);
      }
      puts("");
    }
  }

  Matrix operator * (const Matrix b) {
    Matrix res;
    if (c != b.r) {
      puts("error");
      exit(0);
    }
    res.r = r; res.c = b.c;
    //print(); b.print();
    for (int i = 1; i <= r; i++)
      for (int j = 1; j <= c; j++)
        for (int k = 1; k <= c; k++)
          res.a[i][j] = (res.a[i][j] + a[i][k] * b.a[k][j]) % mod;
    //res.print();
    return res;
  }
};

然后快速幂自然也是很简单的啦

Matrix matrix_pow(Matrix m, ll n) {
  Matrix res;
  res.r = res.c = m.r;
  for (int i = 1; i <= res.r; i++) res.a[i][i] = 1;//单位矩阵
  while (n) {
    if (n & 1)  res = res * m;
    m = m * m;
    n >>= 1;
  }
  return res;
}

附上完整代码

点击查看代码

#include <iostream>
#include <cstring>
#define ll long long
#define mod 1000000007
#define maxn 3
using namespace std;

struct Matrix {
  ll a[maxn][maxn];
  ll r, c;

  Matrix() {
    memset(a, 0, sizeof(a));
  }

  void print() const { // for debug
    if (r < 0 || c < 0) return;
    for (int i = 1; i <= r; i++) {
      for (int j = 1; j <= c; j++) {
        printf("%lld ", a[i][j]);
      }
      puts("");
    }
  }

  Matrix operator * (const Matrix b) {
    Matrix res;
    if (c != b.r) {
      puts("error");
      exit(0);
    }
    res.r = r; res.c = b.c;
    //print(); b.print();
    for (int i = 1; i <= r; i++)
      for (int j = 1; j <= c; j++)
        for (int k = 1; k <= c; k++)
          res.a[i][j] = (res.a[i][j] + a[i][k] * b.a[k][j]) % mod;
    //res.print();
    return res;
  }
};

Matrix matrix_pow(Matrix m, ll n) {
  Matrix res;
  res.r = res.c = m.r;
  for (int i = 1; i <= res.r; i++) res.a[i][i] = 1;
  while (n) {
    if (n & 1)  res = res * m;
    m = m * m;
    n >>= 1;
  }
  return res;
}

ll n;
Matrix f, base, I;
void init() {
  f.r = 1; f.c = 2;
  f.a[1][1] = f.a[1][2] = 1;
  base.r = base.c = 2;
  base.a[1][1] = base.a[1][2] = base.a[2][1] = 1;
  base.a[2][2] = 0;
}


int main() {
  while (cin >> n) {
    init();
    //(f * matrix_pow(base, n - 2)).print();
    if (n > 2) cout << (f * matrix_pow(base, n - 2)).a[1][1] << endl;
    else cout << 1 << endl;
  }
  return 0;
}