前言
刷題時正好遇到這方面的知識，以前學過，但沒寫過博文，忘得差不多了，就重新學下。
找了個基礎題：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;
}