矩阵乘法小记

目录

移植于原csdn博客

基础？递推

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

矩阵乘法

[C_{i,j}=sum A_{i,k} imes B_{k,j} ]

矩阵乘法是具有结合律的，所以通过这个特点可以用到快速幂。

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

洛谷 1962 斐波那契数列

分析

容易可得，状态矩阵就是(left[egin{matrix}F_{n-1}&F_nend{matrix} ight])
根据斐波那契数列，可以得到

[left[egin{matrix}0&1\ 1&1end{matrix} ight] ]

显然，对于(F_{n-1})在下一时刻会给(F_{n})，而(F_{n})会结合(F_{n-1})在下一时刻变成(F_{n+1}),所以得到了这个转移矩阵
那么

[F_{n}=F_{0} imes left[egin{matrix}0&1\ 1&1end{matrix} ight]^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} mod 10000007)

分析

首先想要维护前缀和，但是233怎么办，可以假使(a_{0,0}=23)，那么之后也就是(a_{0,i}=a_{0,i-1} imes 10+3)
那么状态矩阵是

[left[ egin{matrix} 23&a_{1,0}&a_{2,0}&3 end{matrix} ight] ]

想到转移矩阵是

[left[ egin{matrix} 10&0&0&1\ 10&1&0&1\ 10&1&1&1\ 0&0&0&1 end{matrix} ight] ]

所以

[F_{m}=F_{0} imes left[egin{matrix}10&0&0&1\10&1&0&1\10&1&1&1\0&0&0&1end{matrix} ight]^m ]

但是如何表示矩阵中的1呢，也就是

[left[ egin{matrix} 1&0&0&0\ 0&1&0&0\ 0&0&1&0\ 0&0&0&1 end{matrix} ight] ]

代码

#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
",ans);
    }
    return 0;
}

进阶？？？

洛谷 4159 迷路

题目

windy在有向图中迷路了。 该有向图有(N)个节点，windy从节点1出发，他必须恰好在(T)时刻到达节点(N)。 现在给出该有向图，你能告诉windy总共有多少种不同的路径吗？ 注意：windy不能在某个节点逗留，且通过某有向边的时间严格为给定的时间。

分析

首先边权是0或者1的时候，很容易想到矩阵乘法，然而边权是(0sim 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 imes Fib(i)}{n(n+1)div2} ]

下面用逆元就可以解决，但是上面尝试化简
首先要

[sum_{i=1}^n fib(i)=fib(1)+fib(2)+dots +fib(n-1)+fib(n) ]

[=1+fib(1)+fib(2)+dots+fib(n-1)+fib(n)-1 ]

[=2fib(2)+fib(1)+dots+fib(n-1)+fib(n)-1=dots=fib(n+2)-1 ]

那么

[sum_{i=1}^n i imes fib(i)=nsum_{i=1}^nfib(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}^nsum_{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
",ans>=mod?ans-mod:ans);
    }
    return 0;
}

洛谷 3216 数学作业

分析

根据dp，可以知道

[F_n=F_{n-1}+10^{lfloor log n floor+1}+n ]

那么转移矩阵即为

[left[egin{matrix}10^k&0&0\1&1&0\1&1&1end{matrix} ight] ]

然而(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
",ANS.p[2][0]);
}

洛谷 5110 块速递推

题目

(a_n=233*a_{n-1}+666*a_{n-2},a_0=0,a_1=1),多组数据求(a_n)

分析

特征方程做法
可以轻松列出转移矩阵是

[left[egin{matrix}233&666\1&0end{matrix} ight] ]

但是显然是TLE的(O(Tn^3logm))
二进制优化后是(O(n^3logm+Tn^2logm))
但是(T)太大了，可能还是会T，至少我没打这种方法
因为(a^b=a^{bmod varphi(p)}(mod p)[gcd(a,p)==1]),所以可以把时间复杂度优化到(O(n^3logmod+Tn^2logmod))
PS：这题是因为出题人模数设置的好，不然不能这么乱用欧拉定理
然而这样还是会TLE，尝试用矩阵分块，设(k=lceilsqrt{mod} ceil)
那么其实(a^b=(a^k)^{lfloorfrac{b}{k} floor}+a^{bmod k})，两个东西都可以预处理，从而把时间复杂度优化到(O(n^3sqrt 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的解法

[left[egin{matrix}A&E\0&Eend{matrix} ight] ]

即为转移矩阵，虽然解释可能会有点问题，但是惊人的发现

[left[egin{matrix}A&E\0&Eend{matrix} ight]^k=left[egin{matrix}A&E+sum_{i=1}^kA_k\0&Eend{matrix} ight] ]

所以最后还要减掉一个单位矩阵

代码

#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 mod r} imes (a^r)^{lfloorfrac{b}{r} floor})，所以时间复杂度应该是(O(n^3r+Tn^3))。这个(r)就非常关键了
然而空间相应会耗损较大,为(O(n^2r))，所以要权衡利弊，二进制拆分会比较平衡，但也非常容易被卡，貌似还有(k)进制拆分这种操作，然而我啥也不会
完结，撒花

【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/Spare-No-Effort/p/13854990.html