矩阵快速幂
前言
神仙算法
前置知识：矩阵
1. 矩阵怎么快速幂
矩阵的幂与数字的幂一样，就是自己和自己相乘。求一个矩阵的 \(n\) 次幂同样可以用快速幂加速。
首先我们打上数字快速幂的模板：
ll ksm(ll n,ll k){
ll ans=1,f=n;
while(k){
if((k&1)==1)ans=ans*f;
f=f*f;
k=k>>1;
}
return ans;
}
在加入一个重载过 * 运算符的支持乘法运算的矩阵结构体：
struct matrix{
ll data[105][105];ll n,m;
matrix()memset(data,0,sizeof data);
matrix(ll line,ll row,ll k){
n=line;m=row;
if(m!=n){
cout<<"invaild initializing"<<'\n';
return;
}
for(int i=1;i<=n;i++)data[i][i]=k;
}
matrix operator*(const matrix &x)const{//重载 * 操作
if(m!=x.n){
cout<<"invaild multi operator"<<'\n';
matrix res(1,1,1);
return res;
}
matrix c;
c.n=n;
c.m=x.m;
for(int i=1;i<=n;i++)
for(int j=1;j<=m;j++)
for(int k=1;k<=x.m;k++)
c.data[i][k]+=data[i][j]*x.data[j][k];
return c;
}
};
将原来快速幂中的数据类型改成矩阵就可以了。
模板代码：
#define rd read()
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll read(){
ll x=0,f=1;
char c=getchar();
while(c>'9'||c<'0'){if(c=='-') f=-1;c=getchar();}
while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+(c^48);c=getchar();}
return x*f;
}
struct matrix{
ll data[105][105];ll n,m;
matrix()memset(data,0,sizeof data);
matrix(ll line,ll row,ll k){
n=line;m=row;
if(m!=n){
cout<<"invaild initializing"<<'\n';
return;
}
for(int i=1;i<=n;i++)data[i][i]=k;
}
matrix operator*(const matrix &x)const{//重载 * 操作
if(m!=x.n){
cout<<"invaild multi operator"<<'\n';
matrix res(1,1,1);
return res;
}
matrix c;
c.n=n;
c.m=x.m;
for(int i=1;i<=n;i++)
for(int j=1;j<=m;j++)
for(int k=1;k<=x.m;k++)
c.data[i][k]+=data[i][j]*x.data[j][k];
return c;
}
};
matrix ksm(matrix n,ll k){
matrix ans(n.n,n.m,1),f=n;
while(k){
if((k&1)==1)ans=ans*f;
f=f*f;
k=k>>1;
}
return ans;
}
int main(){
ll n;
cin>>n;
matrix bc(2,2,1);
bc.data[1][1]=1;bc.data[1][2]=1;
bc.data[2][1]=1;bc.data[2][2]=0;
matrix res=ksm(bc,n);//计算矩阵 bc 的 n 次方
return 0;
}
2. 矩阵快速幂有什么用
仅是计算一个矩阵的幂有点太简单了，我们常用矩阵快速幂来加速递推。
考虑一个经典的递推：斐波那契数列。
使用递推计算斐波那契数列的第 \(n\) 项的时间复杂度为 \(O(n)\)。
但是使用矩阵快速幂，我们可以把时间优化到 \(O(\log n)\)。
当然，一般这样的斐波那契数列的项数非常大，题目中会让对答案取模。
如何用矩阵递推呢？
我们知道，计算斐波那契数列的下一项需要前两项的值，因此先把 \(F_n\) 和 \(F_{n-1}\) 放到一个一行两列的矩阵里：
此时我们如果可以用某种操作让矩阵变成：
就可以完成一次递推。
我们想到矩阵乘法的一些特点：结果是矩阵的每一行列上的积的和。
因此，我们不妨让这个矩阵乘以：
为什么这样乘，我们一列一列考虑：
首先对于结果矩阵第一行第一列，它等于 \(F_n\times 1+F_{n-1}\times 1\)，即 \(F_n+F_{n-1}\)，也就是 \(F_{n+1}\)。
对于结果矩阵第一行第二列，它等于 \(F_n\times 1\)，也就是 \(F_n\)。
于是计算结果为：
变成了我们想要的样子！
而且我们发现，如果让结果矩阵再乘以那个特殊的矩阵，我们还能再完成一次递推。
也就是说，我们将一个加法递推转换为了矩阵乘法，且这个乘法是一个原始矩阵多次乘以同一个特殊矩阵的过程，即乘以这个特殊矩阵的多次幂。
求一个矩阵的多次幂，快速幂算法正好可以快速解决。
因此，我们要求斐波那契数列的第 \(n\) 项，可以先把对应特殊矩阵的多次幂求解出来，最后乘以一个给定的初始值矩阵，答案便在结果矩阵的特定位置。
我自己把那个特殊矩阵称为转移矩阵（有点像 DP 的转移方程）。
简单观察找规律可以发现，如果我们把初始矩阵定为
即：
求解斐波那契第 \(n\) 项，只需要用快速幂求解出转移矩阵的 \(n-2\) 次幂再与初始矩阵相乘，答案便在第一行，第一列。
于是可以写代码了，记得取模：
Luogu P1962 斐波那契数列
#define rd read()
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll read(){
ll x=0,f=1;
char c=getchar();
while(c>'9'||c<'0'){if(c=='-') f=-1;c=getchar();}
while(c>='0'&&c<='9'){x=(x<<3)+(x<<1)+(c^48);c=getchar();}
return x*f;
}
const int mod=1e9+7;
struct matrix{
ll data[6][6];ll n,m;
matrix()memset(data,0,sizeof data);n=0;m=0;
matrix(int n,int m){
memset(data,0,sizeof data);this->n=n;this->m=m;
}
matrix operator *(const matrix &x)const {
matrix res(n,x.m);
for(int i=1;i<=n;i++)
for(int j=1;j<=m;j++)
for(int k=1;k<=x.m;k++)
res.data[i][k]=(res.data[i][k]+data[i][j]*x.data[j][k])%mod;
return res;
}
};
matrix qp(matrix m,ll k){
matrix res=m;
while(k){
if(k&1)res=res*m;
m=m*m;
k>>=1;
}
return res;
}
int main(){
ll n;
cin>>n;
matrix bc(2,2);
bc.data[1][1]=1;bc.data[1][2]=1;
bc.data[2][1]=1;bc.data[2][2]=0;//这样排列可以在视觉上方便给矩阵赋初值
matrix m(1,2);
m.data[1][1]=1;m.data[1][2]=1;
if(n<=2){
cout<<1;
return 0;
}
matrix res=qp(bc,n-2);
matrix ans=res*m;
cout<<ans.data[1][1];
return 0;
}
回顾一下我们的解题过程，在这个转移矩阵中，只有 \(1\) 和 \(0\) 两种元素，因此我们可以看成结果矩阵 要 还是 不要 原矩阵中的某些元素。
在上面的例题中，我们对答案矩阵的每一个，对应转移矩阵的每一列应是什么样子，思考过程如下：
- 对于结果矩阵的第一个，由于我们要在这里存放 \(n+1\) 项斐波那契数列的值，因此我需要 \(F_{n}\) 和 \(f_{n-1}\)，也就是原矩阵的全部两个元素相加。所以，转移矩阵第一列的两个数都是 \(1\)，就代表着 “都要”。
- 对于结果矩阵的第二个，由于我们将来还要继续向下递推 \(n+2\) 项的值，因此还需要 \(n+1\) 项和 \(n\) 项的值。所以，转移矩阵第二列从上到下的两个数分别是 \(1,0\)，表示我要原矩阵中的 \(F_n\)，不要 \(F_{n-1}\)
当然，以上是转移矩阵中只有 \(0,1\) 的情况。如果遇到递推方程中有系数，还可以在转移矩阵上标任何的数，可以表示不取，取了还要变为几倍等更多含义。
某些求数列带系数前缀和的题目甚至不能递推，我们还可能要用一些数学方法间接推导出值。
例题：
LibreOJ - 10221 Fibonacci 前 n 项和
LibreOJ - 10222 佳佳的 Fibonacci
OI-wiki 上带系数的递推式的转移矩阵
迁移自洛谷
浙公网安备 33010602011771号