YbtOJ 「数学基础」第1章 矩阵快速幂
例题1.序列的第 k 个数
等差数列可以直接用乘法求，等比数列写一个普通快速幂。
注意数据范围是 \(10^9\)，而不是 \(109\)。\(k\le 3\) 的情况记得特判。
code
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int mod=200907;
int T,a,b,c,k;
int qpow(int n,int k)
{
int ans=1;
while(k)
{
if(k&1) ans=ans*n%mod;
n=n*n%mod;k>>=1;
}
return ans;
}
signed main()
{
scanf("%lld",&T);
while(T--)
{
scanf("%lld%lld%lld%lld",&a,&b,&c,&k);
if(k==1) {cout<<a<<endl;continue;}
if(k==2) {cout<<b<<endl;continue;}
if(b-a==c-b)
{
int qwq=b-a;
c=(c+qwq*(k-3)%mod)%mod;
cout<<c<<endl;
}
else
{
int qwq=b/a;
c=c*qpow(qwq,k-3)%mod;
cout<<c<<endl;
}
}
return 0;
}
例题2.斐波那契数列
转移矩阵为[1,1,1,0] 脑补换行吧，懒得打 LaTeX (
code
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int mod=1e9+7;
int n;
struct matrix{
int n,m,f[3][3];
matrix(){memset(f,0,sizeof(f));}
};
matrix operator *(const matrix &x,const matrix &y)
{
matrix res;
res.n=x.n,res.m=y.m;
for(int i=1;i<=res.n;i++)
for(int j=1;j<=res.m;j++)
for(int k=1;k<=x.m;k++)
res.f[i][j]=(res.f[i][j]+x.f[i][k]*y.f[k][j]%mod)%mod;
return res;
}
matrix qpow(matrix n,int k)
{
matrix res;
res.n=res.m=2;res.f[1][1]=res.f[2][2]=1;
while(k>0)
{
if(k&1) res=res*n;
n=n*n;k>>=1;
}
return res;
}
signed main()
{
scanf("%lld",&n);
matrix a,b;
a.n=1,a.m=2,a.f[1][1]=a.f[1][2]=1;
b.n=b.m=2;b.f[1][1]=b.f[1][2]=b.f[2][1]=1;
matrix ans=a*qpow(b,n-2);
cout<<ans.f[1][1]<<endl;
return 0;
}
例题3.行为方案
邻接矩阵的 \(k\) 次方代表走 \(k\) 步后的方案数。
对于自爆，将所有点向 \(0\) 连一条边。
code
#include<bits/stdc++.h>
using namespace std;
const int mod=2017,N=105;
int n,m,t;
struct matrix {
int n,m,f[N][N];
matrix(){memset(f,0,sizeof(f));}
}mp;
matrix operator * (const matrix &x,const matrix &y)
{
matrix res;
res.n=x.n,res.m=y.m;
for(int i=0;i<=x.n;i++)
for(int j=0;j<=y.m;j++)
for(int k=0;k<=x.m;k++)
res.f[i][j]=(res.f[i][j]+x.f[i][k]*y.f[k][j]%mod)%mod;
return res;
}
matrix qpow(matrix a,int k)
{
matrix res;res.n=res.m=n;
for(int i=0;i<=n;i++) res.f[i][i]=1;
while(k)
{
if(k&1) res=res*a;
a=a*a;k>>=1;
}
return res;
}
int main()
{
scanf("%d%d",&n,&m);
for(int i=1,u,v;i<=m;i++)
{
scanf("%d%d",&u,&v);
mp.f[u][v]=mp.f[v][u]=1;
}
mp.n=mp.m=n;
for(int i=0;i<=n;i++) mp.f[i][0]=1,mp.f[i][i]=1;
scanf("%d",&t);
matrix ans=qpow(mp,t);
int sum=0;
for(int i=0;i<=n;i++) sum=(sum+ans.f[1][i])%mod;
cout<<sum<<endl;
return 0;
}
例题4.矩阵求和
因为求的是整个矩阵，考虑构造一个包含 \(A\) 的转移矩阵。则矩阵 \(B\) 为：\(\begin{bmatrix}
A & I\\
0&I
\end{bmatrix}\)，其中 \(I\) 为 \(n\times n\) 的单位矩阵。
则 \(B^k\) 的右上角矩阵为 \(\Sigma^{k-1}_{i=0}A^i\)。故答案为 \(B^{k+1}\) 的右上角减 \(I\)。
注意这题坑人，std 取模没判负数所以答案里有 \(-1\)。也就是说你样例右下角那个数跑出来是 \(-1\) 才能 AC。
code
#include<bits/stdc++.h>
using namespace std;
const int N=65;
int n,k,mod;
struct matrix {
int n,m,f[N][N];
matrix(){memset(f,0,sizeof(f));}
};
matrix operator * (const matrix &x,const matrix &y)
{
matrix res;
res.n=x.n,res.m=y.m;
for(int i=1;i<=x.n;i++)
for(int j=1;j<=y.m;j++)
for(int k=1;k<=x.m;k++)
res.f[i][j]=(res.f[i][j]+x.f[i][k]*y.f[k][j]%mod)%mod;
return res;
}
matrix qpow(matrix a,int k)
{
matrix res;res.n=res.m=2*n;
for(int i=1;i<=2*n;i++) res.f[i][i]=1;
while(k)
{
if(k&1) res=res*a;
a=a*a;k>>=1;
}
return res;
}
int main()
{
matrix a;
scanf("%d%d%d",&n,&k,&mod);//mod*=2;
for(int i=1;i<=n;i++)
{
for(int j=1;j<=n;j++)
{
scanf("%d",&a.f[i][j]);
}
}
a.n=a.m=2*n;
for(int i=n+1;i<=2*n;i++)
a.f[i][i]=1;
for(int i=1;i<=n;i++)
a.f[i][i+n]=1;
matrix ans=qpow(a,k+1);
for(int i=1;i<=n;i++) ans.f[i][i+n]=(ans.f[i][i+n]-1)%mod;
for(int i=1;i<=n;i++)
{
for(int j=n+1;j<=2*n;j++)
{
cout<<ans.f[i][j]<<" ";
}
cout<<endl;
}
return 0;
}
例题5.最短路径
书上告诉我们可以把 \(f[i][j][k]=\min\{f[i][t][k-1]+f[t][j][1]\}\) 转化为 \(f[i][j]=\min\{f[i][t]+f[t][j]\}\)。这东西乍一看不太正确的样子，但写完代码发现他想表达的其实是 \(f[i][j][k]=\min\{f[i][t][k-l]+f[t][j][l]\}\) emmm...
这东西长得很像矩阵乘法公式，所以矩阵乘法。res 初始值不能赋成单位矩阵，或许可以直接用原矩阵开始乘。（也可能还有别的方法但窝tcl不会/kk
code
#include<bits/stdc++.h>
using namespace std;
const int N=205;
int n,t,s,e;
int cnt,num[1005];
struct matrix{
int f[N][N];
matrix(){memset(f,0x3f,sizeof(f));}
}a;
matrix operator *(const matrix &x,const matrix &y)
{
matrix res;
for(int k=1;k<=cnt;k++)
for(int i=1;i<=cnt;i++)
for(int j=1;j<=cnt;j++)
res.f[i][j]=min(res.f[i][j],x.f[i][k]+y.f[k][j]);
return res;
}
matrix qpow(matrix n,int k)
{
matrix res=n;
while(k)
{
if(k&1) res=res*n;
n=n*n;k>>=1;
}
return res;
}
int main()
{
scanf("%d%d%d%d",&t,&n,&s,&e);
for(int i=1,u,v,w;i<=n;i++)
{
scanf("%d%d%d",&w,&u,&v);
if(!num[u]) num[u]=++cnt;
if(!num[v]) num[v]=++cnt;
a.f[num[u]][num[v]]=a.f[num[v]][num[u]]=w;
}
a=qpow(a,t-1);
cout<<a.f[num[s]][num[e]]<<endl;
return 0;
}
1.数列求解
构造转移矩阵\(\begin{bmatrix} p & 1\\ q & 0 \end{bmatrix}\)。
code
#include<bits/stdc++.h>
#define int long long
using namespace std;
int p,q,a1,a2,n,m;
struct matrix{
int n,m;
int f[4][4];
matrix(){memset(f,0,sizeof(f));}
};
matrix operator *(const matrix &x,const matrix &y)
{
matrix res;
res.n=x.n,res.m=y.m;
for(int i=1;i<=x.n;i++)
for(int j=1;j<=y.m;j++)
for(int k=1;k<=x.m;k++)
res.f[i][j]=(res.f[i][j]+x.f[i][k]*y.f[k][j])%m;
return res;
}
matrix qpow(matrix n,int k)
{
matrix res;
res.n=res.m=2;res.f[1][1]=res.f[2][2]=1;
while(k)
{
if(k&1) res=res*n;
n=n*n;k>>=1;
}
return res;
}
signed main()
{
scanf("%lld%lld%lld%lld%lld%lld",&p,&q,&a1,&a2,&n,&m);
matrix a,b;
a.n=1,a.m=2;a.f[1][1]=a2,a.f[1][2]=a1;
b.n=b.m=2;
b.f[1][1]=p,b.f[1][2]=1,b.f[2][1]=q,b.f[2][2]=0;
matrix ans=a*qpow(b,n-2);
cout<<ans.f[1][1]<<endl;
return 0;
}
2.填充棋盘
经找规律，答案为 \(6\times (2^n+2^m)-24\)。
code
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int mod=1e9+7;
int n,m;
int qpow(int n,int k)
{
int ans=1;
while(k)
{
if(k&1) ans=ans*n%mod;
n=n*n%mod;k>>=1;
}
return ans;
}
signed main()
{
scanf("%lld%lld",&n,&m);
cout<<(6*(qpow(2,n)+qpow(2,m))%mod-24+mod)%mod<<endl;
return 0;
}
3.集合统计
为什么这题没题解啊/qd
分别考虑集合中的每一个元素 \(a_i\)。对于我们选择的子集，要保证至少有一个不含 \(a_i\)。则对于这个元素，选择方案数为 \(2^k-1\)。共有 \(n\) 个元素。故最后的答案为 \((2^k-1)^n\)。
这里数据范围是 \(10^{63}\) 超出了 long long 的范围，用欧拉定理处理一下就好了。
hbc学长太强辣！！！
code
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int mod=1e9+7;
int qpow(int n,int k)
{
int ans=1;
while(k)
{
if(k&1) ans=ans*n%mod;
n=n*n%mod;k>>=1;
}
return ans;
}
unsigned long long n,k;
int read()
{
int xr=0;char cr;
cr=getchar();
while(cr<'0'||cr>'9') cr=getchar();
while(cr>='0'&&cr<='9') xr=(xr*10%(mod-1)+(cr-'0'))%(mod-1),cr=getchar();
return xr;
}
signed main()
{
n=read();k=read();
//cout<<n<<" "<<k<<endl;
n%=(mod-1),k%=(mod-1);
//cout<<n<<" "<<k<<endl;
cout<<qpow((qpow(2,k)-1+mod)%mod,n)<<endl;
return 0;
}
4.公式推导
两两配对，组合数常见套路。
然后就可以得到答案为 \((2s+dn)2^{n-1}\)。
这题又会爆 long long，要写龟速乘。
code
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int mod=998244353;
int n,s,d;
int qpow(int n,int k)
{
int ans=1;
while(k)
{
if(k&1) ans=ans*n%mod;
n=n*n%mod;k>>=1;
}
return ans;
}
int mul(int n,int k)
{
int ans=0;
while(k)
{
if(k&1) ans=(ans+n)%mod;
n=(n+n)%mod;k>>=1;
}
return ans;
}
signed main()
{
scanf("%lld%lld%lld",&n,&s,&d);
cout<<(2*s+mul(n,d)%mod)%mod*qpow(2,n-1)%mod<<endl;
return 0;
}
5.斐波拉契
写一下做这题的心路历程吧。
根据题目和题里的式子联想到斐波那契的通项公式，移项得答案为斐波那契数列对应位置的数乘一个 \(\sqrt5\)。提交，后面一半WA了。精度会炸。
用前面小数据打表。发现答案序列也是一个类似斐波那契的东西（？
然后又WA了。然后发现偶数答案要减1。于是过了。
code
#include<bits/stdc++.h>
#define int long long
using namespace std;
int n,mod;
struct matrix{
int n,m,f[3][3];
matrix(){memset(f,0,sizeof(f));}
};
matrix mul(matrix x,matrix y)
{
matrix res;
res.n=x.n,res.m=y.m;
for(int i=1;i<=res.n;i++)
for(int j=1;j<=res.m;j++)
for(int k=1;k<=x.m;k++)
res.f[i][j]=(res.f[i][j]+x.f[i][k]*y.f[k][j]%mod)%mod;
return res;
}
matrix qpow(matrix n,int k)
{
matrix res;
res.n=res.m=2;res.f[1][1]=res.f[2][2]=1;
while(k>0)
{
if(k&1) res=mul(res,n);
n=mul(n,n);k>>=1;
}
return res;
}
signed main()
{
scanf("%lld%lld",&n,&mod);
matrix a,b;
a.n=1,a.m=2,a.f[1][1]=3,a.f[1][2]=1;
b.n=b.m=2;b.f[1][1]=b.f[1][2]=b.f[2][1]=1;
if(n==0) {cout<<1<<endl;return 0;}
matrix ans=mul(a,qpow(b,n-2));
double qwq=sqrt(5.0);
cout<<ans.f[1][1]-(n%2==0)<<endl;
//cout<<(int)(ans.f[1][1]*qwq)%mod<<endl;
return 0;
}
本文来自博客园，作者：樱雪喵，转载请注明原文链接：https://www.cnblogs.com/ying-xue/p/16583395.html
