矩阵快速幂详解（以斐波那契数列为例）
前言
刷题时正好遇到这方面的知识，以前学过，但没写过博文，忘得差不多了，就重新学下。
找了个基础题：https://www.luogu.com.cn/problem/P1962
以求斐波那契数列为例，正常操作是直接循环，时间复杂度\(O(n)\)，然而使用矩阵快速幂时间复杂度为\(O(logn)\)
快速幂
这部分较为简单，重点为下面的公式
\[a^n=\begin{cases}
a^{n/2}*a^{n/2} & \quad \text{if } n \text{ is even}\\
a^{(n-1)/2}*a^{(n-1)/2}*a & \quad \text{if } n \text{ is odd}
\end{cases}
\]
例如求\(2^{18} = 2^9*2^9\)，然后\(2^9 = 2^4*2^4*2\)，接下来\(2^4 = 2^2*2^2\)，最后求\(2^2\)，求一个18次方也仅仅需要4步即可，依次求\(2->2^2->2^4->2^9->2^{18}\)，所以时间复杂度仅为\(O(logn)\)
代码如下：
int quick_pow(int a,int b){
int ans = 1;
while(b){
if(b&1) ans *= a;
a *= a;
b >>= 1;
}
return ans;
}
矩阵快速幂
首先假设\(F_{n}\)为斐波那契数列第n项，矩阵\(Fib(n) = \begin{bmatrix}F_{n} &F_{n-1}\end{bmatrix}\)，矩阵\(base = \begin{bmatrix}a & b\\
c&d
\end{bmatrix}\)
给出一个矩阵公式
\[Fib(n)=Fib(n-1)*base
\]
即
\[\begin{bmatrix}F_{n} &F_{n-1}\end{bmatrix} = \begin{bmatrix}F_{n-1} &F_{n-2}\end{bmatrix} * \begin{bmatrix}a & b\\
c&d
\end{bmatrix}
\]
可得
\[\begin{cases}
a*F_{n-1}+c*F_{n-2} = F_{n}\\
b*F_{n-1}+d*F_{n-2} = F_{n-1}
\end{cases}
\]
显然由上式可计算出\(base = \begin{bmatrix}1 & 1\\
1&0
\end{bmatrix}\)
最终获得公式
\[\begin{bmatrix}F_{n} &F_{n-1}\end{bmatrix}
= \begin{bmatrix}F_{n-1} &F_{n-2}\end{bmatrix} *
\begin{bmatrix}1 & 1\\1&0\end{bmatrix}\]
基本到尾声了，对于斐波那契数列来说，\(F_{1}=F_{2}=1\)，则对于\(Fib(3)\)，得
\[\begin{aligned}
\begin{bmatrix}F_{3} &F_{2}\end{bmatrix} &= \begin{bmatrix}F_{2} &F_{1}\end{bmatrix} * \begin{bmatrix}1 & 1\\
1&0
\end{bmatrix}\\
&= \begin{bmatrix}1 &1\end{bmatrix} * \begin{bmatrix}1 & 1\\
1&0
\end{bmatrix}
\end{aligned}
\]
所以，当\(n>2\)时，可得
\[Fib(n) =
\begin{bmatrix}F_{n} &F_{n-1}\end{bmatrix} = \begin{bmatrix}1 &1\end{bmatrix} * \begin{bmatrix}1 & 1\\
1&0
\end{bmatrix}^{n-2}
\]
最后，\(\begin{bmatrix}1 &1\end{bmatrix}\)与\(\begin{bmatrix}1 & 1\\ 1&0 \end{bmatrix}\)的第一行正好相同，我们也只需要第一个数字，所以最终正好能简化成以下公式
\[\begin{bmatrix}F(n)&F(n-1)\\F(n-1)&F(n-2)\end{bmatrix}=
\begin{bmatrix}1 & 1\\
1&0
\end{bmatrix}^{n-1}
\]
最终代码如下
#include <iostream>
#include <cstring>
#define Max_rank 3
#define mod 1000000007
struct Matrix {
long long a[Max_rank][Max_rank];
Matrix() {
memset(a, 0, sizeof(a));
}
void init(){
a[1][1] = a[1][2] = a[2][1] = 1;
a[2][2] = 0;
}
Matrix operator*(const Matrix b) {
Matrix res;
for (int i = 1; i <= 2; i++)
for (int j = 1; j <= 2; j++)
for (int u = 1; u <= 2; u++)
res.a[i][j] = (res.a[i][j] + a[i][u]*b.a[u][j])%mod;
return res;
}
};
long long q_pow(long long n){
Matrix ans,base;
ans.init();
base.init();
while(n > 0){
if(n&1) ans =ans *base;
base = base *base;
n >>= 1;
}
return ans.a[1][1];
}
int main() {
long long n;
while(std::cin >> n){
std::cout << q_pow(n-2) << std::endl;
}
return 0;
}
