算法总结—矩阵加速
矩阵
矩阵类似于二维数组：
\(A=\begin{bmatrix}a&b\\ c&d\end{bmatrix}\)
矩阵的用途
-
解方程。
-
加速线性 DP。
-
加速 Floyd。
常数乘法
\(Ak=\begin{bmatrix}ak&bk\\ ck&dk\end{bmatrix}\)
矩阵乘法
若 \(A=\begin{bmatrix}a&b\\ c&d\end{bmatrix}\)，\(B=\begin{bmatrix}x&y&z\\ u&v&w\end{bmatrix}\)
则 \(C=A\times B=\begin{bmatrix}ax+bu&ay+bv&az+bw\\ cx+du&cy+dv&cz+dw\end{bmatrix}\)
也就是 \(C_{i,j}=\sum_{i=1}^{A\tt{的列数}} A_{i,k}B_{k,j}\)，所以在矩阵乘法中，\(A\) 的列数要等与 \(B\) 的行数。
单位矩阵
单位矩阵就是形如 \(\begin{bmatrix}1&0&\cdots&0\\0&1&\cdots&0\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&1\end{bmatrix}\) 的矩阵，任何矩阵乘单位矩阵都等于自己，就相当于整数中的 \(1\)。
矩阵快速幂
矩阵快速幂与整数快速幂类似，唯一的不同是需要重载一下矩阵的乘法。
【模板】矩阵快速幂
代码
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int mod=1e9+7;
int n,p;
struct matrix{
int n,m;
int a[105][105];
matrix():n(0),m(0){
memset(a,0,sizeof(a));
return;
}
matrix(int x,int y):n(x),m(y){
memset(a,0,sizeof(a));
return;
}
matrix operator*(const matrix&b){
matrix c(n,b.m);
for(int i=1;i<=n;i++){
for(int j=1;j<=b.m;j++){
for(int k=1;k<=m;k++){
c.a[i][j]=(c.a[i][j]+a[i][k]*b.a[k][j])%mod;
}
}
}
return c;
}
};
matrix qpow(matrix a,int p){
matrix ans(n,n);
for(int i=1;i<=n;i++){
ans.a[i][i]=1;
}
while(p){
if(p&1){
ans=ans*a;
}
a=a*a;
p/=2;
}
return ans;
}
signed main(){
cin>>n>>p;
matrix A(n,n);
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
cin>>A.a[i][j];
}
}
matrix ans=qpow(A,p);
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
cout<<ans.a[i][j]<<" ";
}
cout<<"\n";
}
return 0;
}
广义斐波那契数列
可以把 \(\begin{bmatrix}a_1&a_2\end{bmatrix}\) 看作初始的矩阵，想 \(\begin{bmatrix}a_1&a_2\end{bmatrix}\) 乘上什么矩阵等于 \(\begin{bmatrix}a_2&a_3\end{bmatrix}\)。
首先看 \(a_1\times?+a_2\times?=a_2\)，\(a_1\times 0+a_2\times 1=a_2\)。
再看 \(a_1\times?+a_2\times?=a_3\)，\(a_1\times q+a_2\times p=a_3\)。
所以 \(\begin{bmatrix}a_1&a_2\end{bmatrix}\times\begin{bmatrix}0&q\\1&p\end{bmatrix}=\begin{bmatrix}a_2&a_3\end{bmatrix}\)。
那么 \(\begin{bmatrix}a_n&a_{n+1}\end{bmatrix}=\begin{bmatrix}a_1&a_2\end{bmatrix}\times\begin{bmatrix}0&q\\1&p\end{bmatrix}^{n-1}\)。
而矩阵又可以进行快速幂，所以 \(a_n\) 就可以在 \(\mathcal{O}(2^3\times \log n)\) 的时间内算出来了。
代码
#include<bits/stdc++.h>
#define int long long
using namespace std;
int p,q,x,y,n,mod;
struct matrix{
int n,m;
int a[5][5];
matrix():n(0),m(0){
memset(a,0,sizeof(a));
return;
}
matrix(int x,int y):n(x),m(y){
memset(a,0,sizeof(a));
return;
}
matrix operator*(const matrix&b){
matrix c(n,b.m);
for(int i=1;i<=n;i++){
for(int j=1;j<=b.m;j++){
for(int k=1;k<=m;k++){
c.a[i][j]=(c.a[i][j]+a[i][k]*b.a[k][j])%mod;
}
}
}
return c;
}
};
matrix qpow(matrix a,int p){
matrix ans(2,2);
for(int i=1;i<=2;i++){
ans.a[i][i]=1;
}
while(p){
if(p&1){
ans=ans*a;
}
a=a*a;
p/=2;
}
return ans;
}
signed main(){
cin>>p>>q>>x>>y>>n>>mod;
if(n==1){
cout<<x;
return 0;
}
if(n==2){
cout<<y;
return 0;
}
matrix A(1,2);
A.a[1][1]=x,A.a[1][2]=y;
matrix B(2,2);
B.a[1][1]=0,B.a[1][2]=q;
B.a[2][1]=1,B.a[2][2]=p;
matrix C=qpow(B,n-1);
A=A*C;
cout<<A.a[1][1];
return 0;
}
骨牌
这道题显然是 DP。
\(dp_1=1\)
\(dp_2=2\)
\(dp_3=4\)
\(dp_i=dp_{i-1}+dp_{i-2}+dp_{i-3}(i\geq 4)\)
但是 \(n\) 可以达到 \(10^{18}\)，所以不能 \(\mathcal{O}(n)\) 写，但是可以用矩阵加速。
可以把 \(\begin{bmatrix}dp_1&dp_2&dp_3\end{bmatrix}\) 看作初始的矩阵，想 \(\begin{bmatrix}dp_1&dp_2&dp_3\end{bmatrix}\) 乘上什么矩阵等于 \(\begin{bmatrix}dp_2&dp_3&dp_4\end{bmatrix}\)。
首先看 \(dp_1\times?+dp_2\times?+dp_3\times?=dp_2\)，\(dp_1\times0+dp_2\times1+dp_3\times0=dp_2\)。
再看 \(dp_1\times?+dp_2\times?+dp_3\times?=dp_3\)，\(dp_1\times0+dp_2\times0+dp_3\times1=dp_3\)。
最后看 \(dp_1\times?+dp_2\times?+dp_3\times?=dp_4\)，\(dp_1\times1+dp_2\times1+dp_3\times1=dp_4\)。
所以 \(\begin{bmatrix}dp_1&dp_2&dp_3\end{bmatrix}\times\begin{bmatrix}0&0&1\\1&0&1\\0&1&1\end{bmatrix}=\begin{bmatrix}dp_2&dp_3&dp_4\end{bmatrix}\)。
那么 \(\begin{bmatrix}dp_n&dp_{n+1}&dp_{n+2}\end{bmatrix}=\begin{bmatrix}dp_1&dp_2&dp_3\end{bmatrix}\times\begin{bmatrix}0&0&1\\1&0&1\\0&1&1\end{bmatrix}^{n-1}\)。
这样，所以 \(dp_n\) 就可以在 \(\mathcal{O}(3^3\times \log n)\) 的时间内算出来了。
代码
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int mod=1e9+7;
struct matrix{
int n,m;
int a[5][5];
matrix():n(0),m(0){
memset(a,0,sizeof(a));
return;
}
matrix(int x,int y):n(x),m(y){
memset(a,0,sizeof(a));
return;
}
matrix operator*(const matrix&b){
matrix c(n,b.m);
for(int i=1;i<=n;i++){
for(int j=1;j<=b.m;j++){
for(int k=1;k<=m;k++){
c.a[i][j]=(c.a[i][j]+a[i][k]*b.a[k][j])%mod;
}
}
}
return c;
}
};
matrix qpow(matrix a,int p){
matrix ans(3,3);
for(int i=1;i<=3;i++){
ans.a[i][i]=1;
}
while(p){
if(p&1){
ans=ans*a;
}
a=a*a;
p/=2;
}
return ans;
}
signed main(){
int n;
cin>>n;
matrix A(1,3);
A.a[1][1]=1,A.a[1][2]=2,A.a[1][3]=4;
matrix B(3,3);
B.a[1][3]=B.a[2][1]=B.a[2][3]=B.a[3][2]=B.a[3][3]=1;
matrix C=qpow(B,n-1);
A=A*C;
cout<<A.a[1][1];
return 0;
}
[NOI2012] 随机数生成器
初始矩阵 \(\begin{bmatrix}X_0&c\end{bmatrix}\)，根据上面的方法 \(\begin{bmatrix}X_0&c\end{bmatrix}\times\begin{bmatrix}a&0\\1&1\end{bmatrix}=\begin{bmatrix}X_1&c\end{bmatrix}\)，\(\begin{bmatrix}X_n&c\end{bmatrix}=\begin{bmatrix}X_0&c\end{bmatrix}\times\begin{bmatrix}a&0\\1&1\end{bmatrix}^n\)。
但这道题很坑，\(X_0,m\leq 10^{18}\)，矩阵相乘时，会爆 long long ，所以要用到龟速乘，边乘边取模，龟速乘的时间复杂度是 \(\log\) 级别的，代码如下：
int qmul(int a,int p){
int ans=0;
while(p){
if(p&1){
ans=(ans+a)%mod;
}
a=(a+a)%mod;
p/=2;
}
return ans;
}
代码
#include<bits/stdc++.h>
#define int long long
using namespace std;
int mod,a,c,x0,n,g;
int qmul(int a,int p){
int ans=0;
while(p){
if(p&1){
ans=(ans+a)%mod;
}
a=(a+a)%mod;
p/=2;
}
return ans;
}
struct matrix{
int n,m;
int a[5][5];
matrix():n(0),m(0){
memset(a,0,sizeof(a));
return;
}
matrix(int x,int y):n(x),m(y){
memset(a,0,sizeof(a));
return;
}
matrix operator*(const matrix&b){
matrix c(n,b.m);
for(int i=1;i<=n;i++){
for(int j=1;j<=b.m;j++){
for(int k=1;k<=m;k++){
c.a[i][j]=(c.a[i][j]+qmul(a[i][k],b.a[k][j]))%mod;
}
}
}
return c;
}
};
matrix qpow(matrix a,int p){
matrix ans(2,2);
for(int i=1;i<=2;i++){
ans.a[i][i]=1;
}
while(p){
if(p&1){
ans=ans*a;
}
a=a*a;
p/=2;
}
return ans;
}
signed main(){
cin>>mod>>a>>c>>x0>>n>>g;
matrix A(1,2);
A.a[1][1]=x0,A.a[1][2]=c;
matrix B(2,2);
B.a[1][1]=a,B.a[1][2]=0;
B.a[2][1]=1,B.a[2][2]=1;
matrix C=qpow(B,n);
A=A*C;
cout<<A.a[1][1]%g;
return 0;
}
浙公网安备 33010602011771号