矩阵乘法专题

基础？递推

斐波那契，$F_0=1,F_1=1,F_n=F_{n-2}+F_{n-1}(n\ge 2)$
所以通过$O(n)$解出

---

矩阵乘法

$$C_{i,j}=\sum A_{i,k}\times B_{k,j}$$
矩阵乘法是具有结合律的，所以通过这个特点可以用到快速幂。

• 状态矩阵，长度为$n$的一维向量
• 转移矩阵，与状态矩阵相乘的矩阵

---

洛谷 1962 斐波那契数列

分析

容易可得，状态矩阵就是$\left[\begin{array}{cc}{F}_{n-1}& F_n\end{array}\right]$
根据斐波那契数列，可以得到
$$\left[\begin{array}{cc}0& 1\\ 1& 1\end{array}\right]$$
显然，对于$F_{n-1}$在下一时刻会给$F_n$，而$F_n$会结合$F_{n-1}$在下一时刻变成$F_{n+1}$,所以得到了这个转移矩阵
那么$$F_n=F_0\times {\left[\begin{array}{cc}0& 1\\ 1& 1\end{array}\right]}^n$$
代码太丑，就不贴了

---

HDU 5015 233 Matrix

题目

给定$a_{1,0},a_{2,0},\dots ,a_{n,0}$，然后$a_{0,0}=0,a_{0,1}=233,a_{0,2}=2333,\dots$,规定$a_{i,j}=a_{i-1,j}+a_{i,j-1}$
问$a_{n,m}\mathrm{m}\mathrm{o}\mathrm{d}10000007$

---

分析

首先想要维护前缀和，但是233怎么办，可以假使$a_{0,0}=23$，那么之后也就是$a_{0,i}=a_{0,i-1}\times 10+3$
那么状态矩阵是
$$\left[\begin{array}{cccc}23& a_{1,0}& a_{2,0}& 3\end{array}\right]$$
想到转移矩阵是
$$\left[\begin{array}{cccc}10& 0& 0& 1\\ 10& 1& 0& 1\\ 10& 1& 1& 1\\ 0& 0& 0& 1\end{array}\right]$$
所以
$$F_m=F_0\times {\left[\begin{array}{cccc}10& 0& 0& 1\\ 10& 1& 0& 1\\ 10& 1& 1& 1\\ 0& 0& 0& 1\end{array}\right]}^m$$
但是如何表示矩阵中的1呢，也就是
$$\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

---

代码

#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)$，问期望释放的能量

---

分析

根据概率期望，可以得到答案也就是求
$$\frac{{\sum }_{i=1}^{n}i\times Fib(i)}{n(n+1)÷2}$$
下面用逆元就可以解决，但是上面尝试化简
首先要$$\sum _{i=1}^{n}fib(i)=fib(1)+fib(2)+\cdots +fib(n-1)+fib(n)$$
$$=1+fib(1)+fib(2)+\cdots +fib(n-1)+fib(n)-1$$
$$=2fib(2)+fib(1)+\cdots +fib(n-1)+fib(n)-1=\cdots =fib(n+2)-1$$
那么$$\sum _{i=1}^{n}i\times fib(i)=n\sum _{i=1}^{n}fib(i)-\sum _{i=1}^{n-1}\sum _{j=1}^{i}fib(j)$$
$$=nfib(n+2)-n-\sum _{i=1}^{n-1}fib(i+2)-1=nfib(n+2)-1-\sum _{i=1}^{n-1}fib(i+2)$$
$$=nfib(n+2)+1-\sum _{i=1}^{n+1}fib(i)=nfib(n+2)+1-(fib(n+3)-1)=nfib(n+2)-fib(n+3)+2$$,剩下的就是一个斐波那契数列解决的问题了，代码就不贴了

---

洛谷 5107 能量采集

题目

求${\sum }_{i=1}^n{\sum }_{j=1}^{m}gcd(i,j)$
好的并不是这样
给定一个$n$个点$m$条边的有向图，每个点有初始能量$a_i$。
每过一秒，每个点的能量便会等量地流向所有出边，另外，会有一份流向自己（你可以当做有一个自环）。
现在dkw有$q$次询问，每次询问会给你一个时间$t$ ，dkw想知道$t$秒时每个点的能量。
不保证图中没有重边和自环，答案对998244353取模。

---

分析

可以发现这就是一个矩阵乘法，按照题目建就可以了，时间复杂度$O(q{n}^{3}logt)$，然而会TLE，怎么办，考虑二进制拆分，也就是预处理$2^x$的答案，然后拼凑答案，那么时间复杂度是$O((n^3+q{n}^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，可以知道$$F_n=F_{n-1}+10^{\lfloor \log n\rfloor +1}+n$$
那么转移矩阵即为
$$\left[\begin{array}{ccc}10^k& 0& 0\\ 1& 1& 0\\ 1& 1& 1\end{array}\right]$$
然而$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\ast a_{n-1}+666\ast a_{n-2},a_0=0,a_1=1$,多组数据求$a_n$

---

分析

可以轻松列出转移矩阵是
$$\left[\begin{array}{cc}233& 666\\ 1& 0\end{array}\right]$$
但是显然是TLE的$O(T{n}^{3}logm)$
二进制优化后是$O(n^{3}logm+T{n}^{2}logm)$
但是$T$太大了，可能还是会T，至少我没打这种方法
因为$a^b=a^{b\mathrm{m}\mathrm{o}\mathrm{d}\phi (p)}(\mathrm{m}\mathrm{o}\mathrm{d}p)[gcd(a,p)==1]$,所以可以把时间复杂度优化到$O(n^{3}logmod+T{n}^{2}logmod)$
然而这样还是会TLE，尝试用矩阵分块，设$k=\lceil \sqrt{m}od\rceil$
那么其实$a^b=(a^k{)}^{\lfloor \frac{b}{k}\rfloor }+a^{b\mathrm{m}\mathrm{o}\mathrm{d}k}$，两个东西都可以预处理，从而把时间复杂度优化到$O(n^3\sqrt{k}+T{n}^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}^k{A}_k$

---

分析

这道题按照蓝书的版本应该是分治，但是貌似有点慢，考虑找到转移矩阵，根据dalao的解法
$$\left[\begin{array}{cc}A& E\\ 0& E\end{array}\right]$$
即为转移矩阵，虽然解释可能会有点问题，但是惊人的发现
$${\left[\begin{array}{cc}A& E\\ 0& E\end{array}\right]}^k=\left[\begin{array}{cc}A& E+{\sum }_{i=1}^k{A}_k\\ 0& E\end{array}\right]$$
所以最后还要减掉一个单位矩阵

---

代码

#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^{3}logm)$，空间复杂度$O(n^2)$
对于多组数据，时间复杂度显然是$O(T{n}^{3}logm)$，那么需要通过预处理减少询问时的时间
可以二进制拆分，也就是预处理$2^k$次幂的结果再拼凑，时间复杂度为$O(n^{3}logm+T{n}^{2}logm)$，空间复杂度$O(n^{2}logm)$
当然还有矩阵分块，大致就是$a^b=a^{b\mathrm{m}\mathrm{o}\mathrm{d}r}\times (a^r{)}^{\lfloor \frac{b}{r}\rfloor }$，所以时间复杂度应该是$O(n^{3}r+T{n}^3)$。这个$r$就非常关键了
然而空间相应会耗损较大,为$O(n^{2}r)$，所以要权衡利弊，二进制拆分会比较平衡，但也非常容易被卡，貌似还有$k$进制拆分这种操作，然而我啥也不会
完结，撒花