前言
刷题时正好遇到这方面的知识，以前学过，但没写过博文，忘得差不多了，就重新学下。
找了个基础题：https://www.luogu.com.cn/problem/P1962
以求斐波那契数列为例，正常操作是直接循环，时间复杂度(O(n))，然而使用矩阵快速幂时间复杂度为(O(logn))
快速幂
这部分较为简单，重点为下面的公式
[a^n=egin{cases}
a^{n/2}*a^{n/2} & quad ext{if } n ext{ is even}\
a^{(n-1)/2}*a^{(n-1)/2}*a & quad ext{if } n ext{ 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 = 0;
while(b){
if(b&1) ans *= a;
a *= a;
b >>= 1;
}
return ans;
}
矩阵快速幂
首先假设(F_{n})为斐波那契数列第n项，矩阵(Fib(n) = egin{bmatrix}F_{n} &F_{n-1}end{bmatrix})，矩阵(base = egin{bmatrix}a & b\
c&d
end{bmatrix})
给出一个矩阵公式
[Fib(n)=Fib(n-1)*base
]
即
[egin{bmatrix}F_{n} &F_{n-1}end{bmatrix} = egin{bmatrix}F_{n-1} &F_{n-2}end{bmatrix} * egin{bmatrix}a & b\
c&d
end{bmatrix}
]
可得
[egin{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 = egin{bmatrix}1 & 1\
1&0
end{bmatrix})
最终获得公式
[egin{bmatrix}F_{n} &F_{n-1}end{bmatrix}
= egin{bmatrix}F_{n-1} &F_{n-2}end{bmatrix} *
egin{bmatrix}1 & 1\1&0end{bmatrix}]
基本到尾声了，对于斐波那契数列来说，(F_{1}=F_{2}=1)，则对于(Fib(3))，得
[egin{aligned}
egin{bmatrix}F_{3} &F_{2}end{bmatrix} &= egin{bmatrix}F_{2} &F_{1}end{bmatrix} * egin{bmatrix}1 & 1\
1&0
end{bmatrix}\
&= egin{bmatrix}1 &1end{bmatrix} * egin{bmatrix}1 & 1\
1&0
end{bmatrix}
end{aligned}
]
所以，当(n>2)时，可得
[Fib(n) =
egin{bmatrix}F_{n} &F_{n-1}end{bmatrix} = egin{bmatrix}1 &1end{bmatrix} * egin{bmatrix}1 & 1\
1&0
end{bmatrix}^{n-2}
]
最后，(egin{bmatrix}1 &1end{bmatrix})与(egin{bmatrix}1 & 1\ 1&0 end{bmatrix})的第一行正好相同，我们也只需要第一个数字，所以最终正好能简化成以下公式
[egin{bmatrix}F(n)&F(n-1)\F(n-1)&F(n-2)end{bmatrix}=
egin{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;
}