矩阵乘法小记
移植于原csdn博客
基础？递推
斐波那契，\(F_0=1,F_1=1,F_n=F_{n-2}+F_{n-1}(n\geq 2)\)
所以通过\(O(n)\)解出
矩阵乘法
矩阵乘法是具有结合律的，所以通过这个特点可以用到快速幂。
- 状态矩阵，长度为\(n\)的一维向量
- 转移矩阵，与状态矩阵相乘的矩阵
洛谷 1962 斐波那契数列
分析
容易可得，状态矩阵就是\(\left[\begin{matrix}F_{n-1}&F_n\end{matrix}\right]\)
根据斐波那契数列，可以得到
显然，对于\(F_{n-1}\)在下一时刻会给\(F_{n}\)，而\(F_{n}\)会结合\(F_{n-1}\)在下一时刻变成\(F_{n+1}\),所以得到了这个转移矩阵
那么
HDU 5015 233 Matrix
题目
给定\(a_{1,0},a_{2,0},\dots,a_{n,0}\)，然后
,规定\(a_{i,j}=a_{i-1,j}+a_{i,j-1}\)
询问\(a_{n,m} \bmod 10000007\)
分析
首先想要维护前缀和，但是233怎么办，可以假使\(a_{0,0}=23\)，那么之后也就是\(a_{0,i}=a_{0,i-1}\times 10+3\)
那么状态矩阵是
想到转移矩阵是
所以
但是如何表示矩阵中的1呢，也就是
代码
#include <cstdio>
#include <cctype>
#include <cstring>
#define rr register
using namespace std;
struct maix{int p[13][13];}A,ANS;
const int mod=10000007; int n;
inline signed iut(){
rr int ans=0; rr char c=getchar();
while (!isdigit(c)) c=getchar();
while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();
return ans;
}
inline maix mul(maix A,maix B){
rr maix C;
for (rr int i=0;i<n+2;++i)
for (rr int j=0;j<n+2;++j) C.p[i][j]=0;
for (rr int i=0;i<n+2;++i)
for (rr int j=0;j<n+2;++j)
for (rr int k=0;k<n+2;++k)
C.p[i][j]=(1ll*A.p[i][k]*B.p[k][j]+C.p[i][j])%mod;
return C;
}
signed main(){
while (scanf("%d",&n)==1){
for (rr int i=0;i<n+2;++i)
for (rr int j=0;j<n+2;++j)
ANS.p[i][j]=i==j;
for (rr int i=0;i<n+2;++i)
for (rr int j=0;j<n+2;++j){
if (!j&&i!=n+1) A.p[i][j]=10;
else if ((i>=j&&i!=n+1)||j==n+1) A.p[i][j]=1;
else A.p[i][j]=0;
}
for (rr int m=iut();m;m>>=1,A=mul(A,A))
if (m&1) ANS=mul(ANS,A);
rr int ans=23ll*ANS.p[n][0];
for (rr int i=1;i<=n;++i)
ans=(1ll*iut()*ANS.p[n][i]+ans)%mod;
ans=(3ll*ANS.p[n][n+1]+ans)%mod;
printf("%d\n",ans);
}
return 0;
}
进阶？？？
洛谷 4159 迷路
题目
windy在有向图中迷路了。 该有向图有\(N\)个节点，windy从节点1出发，他必须恰好在\(T\)时刻到达节点\(N\)。 现在给出该有向图，你能告诉windy总共有多少种不同的路径吗？ 注意：windy不能在某个节点逗留，且通过某有向边的时间严格为给定的时间。
分析
首先边权是0或者1的时候，很容易想到矩阵乘法，然而边权是\(0\sim 9\)，考虑拆点，把一个点拆成9个点，再建立转移矩阵，对于每个拆的点，连向前一个拆的点，对于连边，用原来的点连向边权为\(w\)的拆出的点，嗯，就是这个样子
代码
#include <cstdio>
#define rr register
using namespace std;
struct maix{int p[91][91];}A,ANS;
int n,m,t;
inline maix mul(maix A,maix B){
rr maix C;
for (rr int i=1;i<=t;++i)
for (rr int j=1;j<=t;++j) C.p[i][j]=0;
for (rr int i=1;i<=t;++i)
for (rr int j=1;j<=t;++j)
for (rr int k=1;k<=t;++k)
C.p[i][j]=(A.p[i][k]*B.p[k][j]+C.p[i][j])%2009;
return C;
}
signed main(){
scanf("%d%d",&n,&m); t=n*9;
for (rr int i=1;i<=n;++i){
for (rr int j=1;j<9;++j) A.p[i+j*n][i+j*n-n]=1;
for (rr int j=1,x;j<=n;++j){
scanf("%1d",&x);
if (x) A.p[i][j+x*n-n]=1;
}
}
for (rr int i=1;i<=t;++i) ANS.p[i][i]=1;
for (;m;m>>=1,A=mul(A,A))
if (m&1) ANS=mul(ANS,A);
return !printf("%d",ANS.p[1][n]);
}
洛谷 4834 萨塔尼亚的期末考试
题目
有\(n\)个点，假设第1个点被选的概率为\(P\),那么第\(i\)个点就是\(iP\),释放的能量为\(Fib(i)\)，问期望释放的能量
分析
根据概率期望，可以得到答案也就是求
下面用逆元就可以解决，但是上面尝试化简
首先要
那么
剩下的就是一个斐波那契数列解决的问题了，代码就不贴了
洛谷 5107 能量采集
题目
求\(\sum_{i=1}^n\sum_{j=1}^mgcd(i,j)\)
好的并不是这样
给定一个\(n\)个点\(m\)条边的有向图，每个点有初始能量\(a_i\)。
每过一秒，每个点的能量便会等量地流向所有出边，另外，会有一份流向自己（你可以当做有一个自环）。
现在dkw有\(q\)次询问，每次询问会给你一个时间\(t\) ，dkw想知道\(t\)秒时每个点的能量。
不保证图中没有重边和自环，答案对998244353取模。
分析
可以发现这就是一个矩阵乘法，按照题目建就可以了，时间复杂度\(O(qn^3logt)\)，
然而会TLE，怎么办，考虑二进制拆分，也就是预处理\(2^x\)的答案，
然后拼凑答案，那么时间复杂度是\(O((n^3+qn^2)logt)\)，虽然这个时间复杂度好像很难办，但至少可以过。
代码
50分警告还是不可以过？？发现矩阵乘法对于询问时显然是浪费了，于是就有了
还是T掉了，考虑取模优化
还是T掉了，考虑把long long改成int
#include <cstdio>
#include <cctype>
#define rr register
using namespace std;
const int mod=998244353; int n,m,q;
struct maix{int p[51][51];}A[31],T;
inline signed iut(){
rr int ans=0; rr char c=getchar();
while (!isdigit(c)) c=getchar();
while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();
return ans;
}
inline signed ksm(int x,int y){
rr int ans=1;
for (;y;y>>=1,x=1ll*x*x%mod)
if (y&1) ans=1ll*ans*x%mod;
return ans;
}
inline signed mo(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline maix mul(maix A,maix B,int t){
rr maix C;
for (rr int i=0;i<t;++i)
for (rr int j=0;j<n;++j) C.p[i][j]=0;
for (rr int i=0;i<t;++i)
for (rr int j=0;j<n;++j)
for (rr int k=0;k<n;++k)
C.p[i][j]=mo(C.p[i][j],1ll*A.p[i][k]*B.p[k][j]%mod);
return C;
}
signed main(){
n=iut(); m=iut(); q=iut();
for (rr int i=0;i<n;++i) T.p[0][i]=iut(),A[0].p[i][i]=1;
while (m--) ++A[0].p[iut()-1][iut()-1];
for (rr int i=0;i<n;++i){
rr int sum=0;
for (rr int j=0;j<n;++j) sum+=A[0].p[i][j];
sum=ksm(sum,mod-2);
for (rr int j=0;j<n;++j) A[0].p[i][j]=1ll*A[0].p[i][j]*sum%mod;
}
for (rr int i=1;i<30;++i) A[i]=mul(A[i-1],A[i-1],n);
while (q--){
rr int t=iut(); rr maix ANS=T;
for (rr int i=0;i<30;++i)
if ((t>>i)&1) ANS=mul(ANS,A[i],1);
rr int ans=0;
for (rr int i=0;i<n;++i) ans^=ANS.p[0][i];
printf("%d\n",ans>=mod?ans-mod:ans);
}
return 0;
}
洛谷 3216 数学作业
分析
根据dp，可以知道
那么转移矩阵即为
然而\(10^k\)比较难搞，考虑分块，枚举\(k\)把它分成\(k\)个部分依次相乘
代码
#include <cstdio>
#include <cstring>
#define rr register
using namespace std;
typedef long long ll;
struct maix{ll p[3][3];}A,ANS;
ll n; int mod;
inline maix mul(maix A,maix B){
rr maix C; memset(C.p,0,sizeof(C.p));
for (rr int i=0;i<3;++i)
for (rr int j=0;j<3;++j)
for (rr int k=0;k<3;++k)
C.p[i][j]=(A.p[i][k]*B.p[k][j]+C.p[i][j])%mod;
return C;
}
signed main(){
scanf("%lld%d",&n,&mod);
ANS.p[0][0]=ANS.p[1][1]=ANS.p[2][2]=1;
for (rr ll t=10,k;;t=t*10){
if (t>n) k=n-t/10+1; else k=t-t/10;
memset(A.p,0,sizeof(A.p));
A.p[0][0]=t%mod; A.p[1][0]=1; A.p[2][0]=1; A.p[1][1]=1; A.p[2][1]=1; A.p[2][2]=1;
for (;k;k>>=1,A=mul(A,A))
if (k&1) ANS=mul(ANS,A);
if (t>n) break;
}
return !printf("%lld\n",ANS.p[2][0]);
}
洛谷 5110 块速递推
题目
\(a_n=233*a_{n-1}+666*a_{n-2},a_0=0,a_1=1\),多组数据求\(a_n\)
分析
特征方程做法
可以轻松列出转移矩阵是
但是显然是TLE的\(O(Tn^3logm)\)
二进制优化后是\(O(n^3logm+Tn^2logm)\)
但是\(T\)太大了，可能还是会T，至少我没打这种方法
因为\(a^b=a^{b\bmod \varphi(p)}(\bmod p)[gcd(a,p)==1]\),所以可以把时间复杂度优化到\(O(n^3logmod+Tn^2logmod)\)
PS：这题是因为出题人模数设置的好，不然不能这么乱用欧拉定理
然而这样还是会TLE，尝试用矩阵分块，设\(k=\lceil\sqrt{mod}\rceil\)
那么其实\(a^b=(a^k)^{\lfloor\frac{b}{k}\rfloor}+a^{b\bmod k}\)，两个东西都可以预处理，从而把时间复杂度优化到\(O(n^3\sqrt k+Tn^3)\)，然后循环展开就可以草率的理解为\(O(\sqrt k+T)\)
代码
#include <cstdio>
#include <cstring>
#define rr register
using namespace std;
const int N=31623,mod=1000000007;
typedef unsigned long long ull;
struct maix{int p[2][2];}A[N+1],B[N+1];
inline signed mo(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline maix mul(maix A,maix B){
rr maix C;
C.p[0][0]=mo(1ll*A.p[0][0]*B.p[0][0]%mod,1ll*A.p[0][1]*B.p[1][0]%mod);
C.p[0][1]=mo(1ll*A.p[0][0]*B.p[0][1]%mod,1ll*A.p[0][1]*B.p[1][1]%mod);
C.p[1][0]=mo(1ll*A.p[1][0]*B.p[0][0]%mod,1ll*A.p[1][1]*B.p[1][0]%mod);
C.p[1][1]=mo(1ll*A.p[1][0]*B.p[0][1]%mod,1ll*A.p[1][1]*B.p[1][1]%mod);
return C;
}
signed main(){
A[1].p[0][0]=233,A[1].p[0][1]=666,A[0].p[0][0]=A[0].p[1][1]=A[1].p[1][0]=1,B[0]=A[0];
for (rr int i=2;i<=N;++i) A[i]=mul(A[i-1],A[1]);
for (rr int i=1;i<=N;++i) B[i]=mul(B[i-1],A[N]);
rr int t,n,ans=0; rr ull sa,sb,sc,temp; scanf("%d",&t);
for (scanf("%llu %llu %llu",&sa,&sb,&sc);t;--t){
sa^=sa<<32,sa^=sa>>13,sa^=sa<<1,temp=sa,sa=sb,sb=sc,sc^=temp^sa;
n=(sc-1)%(mod-1); rr maix C=mul(B[n/N],A[n%N]); ans^=C.p[0][0];
}
return !printf("%d",ans);
}
POJ 3233 Matrix Power Series
题目
给定一个矩阵\(A\)，求\(\sum_{i=1}^kA_k\)
分析
这道题按照蓝书的版本应该是分治，但是貌似有点慢，考虑找到转移矩阵，根据dalao的解法
即为转移矩阵，虽然解释可能会有点问题，但是惊人的发现
所以最后还要减掉一个单位矩阵
代码
#include <cstdio>
#include <cctype>
#include <cstring>
#define rr register
using namespace std;
struct maix{int p[61][61];}A,ANS; int n,mod,k;
inline signed iut(){
rr int ans=0; rr char c=getchar();
while (!isdigit(c)) c=getchar();
while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();
return ans;
}
inline maix mul(maix A,maix B){
rr maix C; memset(C.p,0,sizeof(C.p));
for (rr int i=1;i<=n;++i)
for (rr int j=1;j<=n;++j)
for (rr int k=1;k<=n;++k)
C.p[i][j]=(A.p[i][k]*B.p[k][j]+C.p[i][j])%mod;
return C;
}
signed main(){
n=iut(); k=iut()+1; mod=iut();
for (rr int i=1;i<=n;++i)
for (rr int j=1;j<=n;++j) A.p[i][j]=iut();
for (rr int i=1;i<=n;++i) A.p[i+n][i+n]=A.p[i][i+n]=1;
n<<=1; for (rr int i=1;i<=n;++i) ANS.p[i][i]=1;
for (;k;k>>=1,A=mul(A,A)) if (k&1) ANS=mul(ANS,A);
n>>=1; for (rr int i=1;i<=n;++i) --ANS.p[i][i+n];
for (rr int i=1;i<=n;++i)
for (rr int j=1;j<=n;++j)
printf("%d%c",ANS.p[i][j+n]<0?ANS.p[i][j+n]+mod:ANS.p[i][j+n],j==n?10:32);
return 0;
}
后续
如何优化卡常
减少取模次数，少用long long，减少不必要的循环，减少矩阵大小,循环展开
对于只有单个数据，时间复杂度显然是\(O(n^3logm)\)，空间复杂度\(O(n^2)\)
对于多组数据，时间复杂度显然是\(O(Tn^3logm)\)，那么需要通过预处理减少询问时的时间
可以二进制拆分，也就是预处理\(2^k\)次幂的结果再拼凑，时间复杂度为\(O(n^3logm+Tn^2logm)\)，空间复杂度\(O(n^2logm)\)
当然还有矩阵分块，大致就是\(a^b=a^{b \bmod r}\times (a^r)^{\lfloor\frac{b}{r}\rfloor}\)，所以时间复杂度应该是\(O(n^3r+Tn^3)\)。这个\(r\)就非常关键了
然而空间相应会耗损较大,为\(O(n^2r)\)，所以要权衡利弊，二进制拆分会比较平衡，但也非常容易被卡，貌似还有\(k\)进制拆分这种操作，然而我啥也不会
完结，撒花
