矩阵乘法优化DP复习
前言
最近做毒瘤做多了……联赛难度的东西也该复习复习了。
Warning:本文较长，难度分界线在“中场休息”部分，如果只想看普及难度的可以从第五部分直接到注意事项qwq
文中用（比如现在这个文本）引用文本书写的部分为总结性内容，即使是跳过部分也建议阅读awa
没事，最难也就NOI2020的签到题，不怕（
0——P1962 斐波那契数列
题意
对于上述递推式，求 \(F(n)\mod 1e9+7\) ，给出的 \(n\leq 2^{63}\) .
思路
递推式给你了。照理讲 \(O(n)\) 是多么优秀的复杂度啊 然而并不。
于是，对于一个清晰明白的递推式，考虑矩乘。
首先是答案矩阵：
然后是转移矩阵。转移矩阵 \(base\) 可以通过：
来确定。保险一点就是设出每个矩阵中的数，然后根据每个项解方程即可 当然也可以看着递推式直接手推 。
但是我个人比较喜欢通过矩阵意义推理。
首先，根据矩阵乘法 c[i][j]=a[i][k]*b[k][j] ，可以知道, \(base\) 的第一行就是 \(F(n)\) 第一列也就是 \(f(n)\) 对下一个答案矩阵的两列的贡献。
具体地讲，\(base\) 的第一行就是 \(f(n)\) 分别对 \(F(n+1)\) 中 \(f(n+1)\) （第一列）和 \(f(n)\) （第二列）贡献了 \(1\times f(n)\) .\(f(n-1)\) 同理。
最后的 \(base\) :
代码
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const ll mod=1e9+7;
struct matrix
{
ll mt[5][5];
matrix() { memset(mt,0,sizeof(mt)); }
matrix operator * ( matrix b )
{
matrix c;
for ( int i=1; i<=2; i++ )
for ( int j=1; j<=2; j++ )
for ( int k=1; k<=2; k++ )
c.mt[i][j]=(c.mt[i][j]+(mt[i][k]*b.mt[k][j])%mod)%mod;
return c;
}
}ans,base;
matrix power( ll b )
{
matrix res=ans,a=base;
for ( ; b; b>>=1,a=a*a )
if ( b&1 ) res=res*a;
return res;
}
int main()
{
ll n; scanf( "%lld",&n );
if ( n<=2 ) { printf( "1\n"); return 0; }
ans.mt[1][1]=1; ans.mt[1][2]=1; ans.mt[2][1]=0; ans.mt[2][2]=0;
base.mt[1][1]=1; base.mt[1][2]=1; base.mt[2][1]=1; base.mt[2][2]=0;
ans=power(n-2);
printf( "%lld\n",ans.mt[1][1] );
}
1——P1939 【模板】矩阵加速（数列）
题意
求 \(a_n\mod 1e9+7\)
思路
推导过程同上，把 \(f(n+1)\) 拆成递推式即可。
答案矩阵：
下一个：
转移矩阵：
代码
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const ll mod=1e9+7;
struct matrix
{
ll mt[5][5];
matrix() { memset(mt,0,sizeof(mt)); }
matrix operator * ( matrix b )
{
matrix c;
for ( int i=1; i<=3; i++ )
for ( int j=1; j<=3; j++ )
for ( int k=1; k<=3; k++ )
c.mt[i][j]=(c.mt[i][j]+(mt[i][k]*b.mt[k][j])%mod)%mod;
return c;
}
}ans,base;
matrix power( ll b )
{
matrix res=ans,a=base;
for ( ; b; b>>=1,a=a*a )
if ( b&1 ) res=res*a;
return res;
}
int main()
{
int T; scanf( "%d",&T );
ans.mt[1][1]=1; ans.mt[1][2]=1; ans.mt[1][3]=1;
base.mt[1][1]=1; base.mt[1][2]=1; base.mt[2][3]=1; base.mt[3][1]=1;
while ( T-- )
{
ll n; scanf( "%lld",&n );
if ( n<=3 ) { printf( "1\n"); continue; }
matrix res=power( n-3 );
printf( "%lld\n",res.mt[1][1] );
}
}
2——HDU5171 GTY's birthday gift
题意
有一个大小为 \(n\) 的可重集 \(S\) ，小奇每次操作可以加入一个数 \(a+b(a,b\in S)\) ，求 \(k\) 次操作后可获得的 \(S\) 的和的最大值。
思路
关键在于求和。涉及到了：
矩阵加维
- 带常数项
\(f(n)=f(n−1)+f(n−2)+k\)
答案矩阵： \(F(n)=(f(n),f(n-1),k)\)
- 带未知项
$f(n)=f(n−1)+f(n−2)+n $
答案矩阵： \(F(n)=(f(n),f(n-1),n,1)\)
- 求和
\(f(n)=f(n-1)+f(n-2),S(n)=\sum_{i=1}^nf(i)\)
答案矩阵： \(F(n)=(f(n),f(n-1),S(n))\)
显然此题属于第三种情况。为什么能递推呢，原因在于要让最后的和最大，又是个可重集，显然每次贪心选最大和次大就好了。
答案矩阵：
下一个：
转移矩阵：
Tips
(负数这个据说 HDU上面没有，但是BZOJ是有的)
当最大和次大都是负数，这两个的和累加 k 次即可。
若次大是负数，那么先把最大和次大累加直到次大大于0（每次k--），然后正常矩乘即可。
（这道题还要注意不同寻常的模数……）
还有就是我代码在HDU上面死活过不去，TLE；然后就跑到（黑暗）BZOJ上去交了，去掉了多测，加上了负数之后不用快读 193ms 跑得飞快……要交哪个OJ自己定吧。（主要dark上面有数据下载，虽然我没用）
代码
//-----------矩乘部分省略，同上-------------------
ll val[100010];
int main()
{
ll n,k,tot=0;
scanf( "%lld%lld",&n,&k );
for ( ll i=1; i<=n; i++ )
scanf( "%lld",&val[i] ),tot+=val[i],tot%=mod;
sort( val+1,val+1+n );
if ( val[n]<=0 )
{
tot+=(val[n]+val[n-1])*k; printf( "%lld\n",tot%mod ); return 0;
}
while ( val[n-1]<0 && k ) val[n-1]=(val[n]+val[n-1])%mod,tot=(tot+val[n-1])%mod,k--;
ans.mt[1][1]=val[n]; ans.mt[1][2]=val[n-1]; ans.mt[1][3]=tot;
base.mt[1][1]=1; base.mt[1][2]=1; base.mt[1][3]=1; base.mt[2][1]=1; base.mt[2][3]=1; base.mt[3][3]=1;
matrix res=power( k );
printf( "%lld\n",res.mt[1][3] );
}
3——Power of Matrix(UVA11149)
题意
给定 \(n,k,A(n\times n)\) ，求 \(A+A^2+\dots+A^k\) 的矩阵个位（就是矩阵每个元素都是一位，也就是 \(\mod =10\) .
思路
终于不用推转移矩阵了（
采用分治的思想。记 \(P_k=A+A^2+\dots +A^k\) ，对 \(k\) 的奇偶性进行讨论：
\(k\mod 2==0\) ,那么 原式 \(=(1+P_{k/2} )\times P_{k/2}\);
\(k\mod 2==1\)，那么可以转化为 \(A+A\times P_{k-1}\) ，回到第一种情况。
代码
惊奇地发现我五个月前写过这道题 为了复习还是再写一遍。
有一说一，这道题练矩阵运算板子还是不错的qwq，但是要知道 UVA 和 LA 对输出的要求之高简直一言难尽……拼命搞多测和严格的格式要求……作为刷训练指南的人已经习惯了（各位加油
然后，结构体内函数如果 没有类型且和结构体同名 ，是一种奇妙的写法，就相当于你每建一个新的结构体会自动运行这部分代码，可以用来初始化。如果你传参数（比如 (int a=0) ），就是默认 a=0 ，如果你写了别的参数就是你写的那个值，这样可以自定义初始化（这题不用）。
//主要部分
struct matrix
{
ll mt[50][50];
matrix() { memset( mt,0,sizeof(mt) ); }
matrix operator + ( const matrix &b )
{
matrix c;
for ( int i=1; i<=n; i++ )
for ( int j=1; j<=n; j++ )
c.mt[i][j]=(mt[i][j]+b.mt[i][j])%mod;
return c;
}
matrix operator * ( const matrix &b )
{
matrix c;
for ( int i=1; i<=n; i++ )
for ( int j=1 ;j<=n; j++ )
for ( int k=1 ;k<=n; k++ )
c.mt[i][j]=(c.mt[i][j]+mt[i][k]*b.mt[k][j]%mod)%mod;
return c;
}
}f[50],ans;
void init( ll x )
{
ll cnt=1;
for ( int i=0; cnt<=x; i++,cnt<<=1 )
f[i+1]=f[i]*f[i];
}
matrix power( ll b )
{
matrix res; ll cnt=0;
for ( ll i=1; i<=n; i++ )
res.mt[i][i]=1;
for ( ; b; b>>=1,cnt++ )
if ( b&1 ) res=res*f[cnt];
return res;
}
matrix split( ll x )
{
if ( x==1 ) return f[0];
matrix res=split( x>>1 );
res=res+res*power(x>>1);
if ( x&1 ) res=res+power(x);
return res;
}
4——How many ways?（HDU2157）
题意
询问有向图上从 \(S\) 点恰好经过 \(k\) 个点到达 \(T\) 点的方案数对 1000 取模的余数。
思路
非常经典的一类题目，矩阵加速DP解决图上问题。这类题目相当于一种套路，在 \(NOI\ Online\) 和 \(NOI2020\) 里面都有出现，后面会讲到。
首先把给定的图转化为邻接矩阵（这注定了这类题目的点数不会很大，也是判断是不是用矩阵加速来做的一个标志）。
令 \(C=A\times A\) ，那么 \(C[i][j]=\sum A[i][k]\times A[k][j]\) ，实际上就等同于 \(i\) 到 \(j\) 恰好经过两条边的路径数量，\(k\) 步同理。
如何理解？一种方式是直接把矩阵乘法的转移看做是 \(Floyd\) 最短路，边权都是1；另一种方式是，考虑每次相乘的两个数，当且仅当都是 1，结果才会是1.放到整体去看就相当于 \((i,k),(k,j)\) 这两条边都是存在的，也就是经过两条边的路径。（相当于 Floyd 枚举中转点的思想）
代码
#include <bits/stdc++.h>
using namespace std;
const int N=110,mod=1000;
int n,m;
struct matrix
{
int mt[N][N];
matrix() { memset( mt,0,sizeof(mt) ); }
matrix operator * ( const matrix &b )
{
matrix c;
for ( int i=0; i<n; i++ )
for ( int j=0; j<n; j++ )
for ( int k=0; k<n; k++ )
c.mt[i][j]=(c.mt[i][j]+mt[i][k]*b.mt[k][j]%mod)%mod;
return c;
}
};
matrix power( matrix a,int b )
{
matrix res;
for ( int i=0; i<n; i++ )
for ( int j=0; j<n; j++ )
res.mt[i][j]=(i==j);
for ( ; b; b>>=1,a=a*a )
if ( b&1 ) res=res*a;
return res;
}
int main()
{
while ( cin>>n>>m && (n||m) )
{
matrix g;
for ( int i=0,x,y; i<m; i++ )
scanf( "%d%d",&x,&y ),g.mt[x][y]=1;
int q; scanf( "%d",&q );
while ( q-- )
{
int x,y,k; scanf( "%d%d%d",&x,&y,&k );
matrix res=power( g,k );
printf( "%d\n",res.mt[x][y] );
}
}
}
5——P2886 [USACO07NOV]Cow Relays G
题意
给出一张无向连通图，求 \(S\) 到 \(E\) 经过 \(k\) 条边的最短路。
思路
跟 4 一样，不过是把矩阵乘法里面的相乘换成了取min，多了一点 Floyd 的味道。注意矩阵初始化也得改改，毕竟是最短路。
然后还得离散化emmm。
代码
#include <bits/stdc++.h>
using namespace std;
const int N=110;
int n,m,k,S,T;
struct matrix
{
int mt[N][N];
matrix() { memset( mt,0x3f,sizeof(mt) ); }
matrix operator * ( const matrix &b )
{
matrix c;
for ( int i=1; i<=n; i++ )
for ( int j=1; j<=n; j++ )
for ( int k=1; k<=n; k++ )
c.mt[i][j]=min(c.mt[i][j],mt[i][k]+b.mt[k][j]);
return c;
}
}g;
matrix power( matrix a,int b )
{
matrix res=a; b--;
for ( ; b; b>>=1,a=a*a )
if ( b&1 ) res=res*a;
return res;
}
int id[1010];
int main()
{
scanf( "%d%d%d%d",&k,&m,&S,&T );
memset( id,0,sizeof(id) );
for ( int i=0,u,v,w; i<m; i++ )
{
scanf( "%d%d%d",&w,&u,&v );
u=id[u] ? id[u] : id[u]=++n;
v=id[v] ? id[v] : id[v]=++n;
g.mt[u][v]=g.mt[v][u]=min( g.mt[u][v],w );
}
S=id[S]; T=id[T];
matrix ans=power( g,k );
printf( "%d",ans.mt[S][T] );
}
中场休息
好的！那么这样我们就把 @[Xing_Ling] 神仙的作业写完了！bushi
现在，要不要再来点甜点呢——
6——[NOI Online3]提高组T3 魔法值
题目链接：Acwing1828 luogu
题意
\(n\) 座城市， \(m\) 条道路，任意两个城市之间之多一条道路。
在第 \(j\) 天，\(i\) 号城市的魔法值为 \(f_{i,j}\) 。给出所有的 \(f_{i,0}\) ，之后每一个城市的魔法值都是：
其中 \(\oplus\) 表示异或，\(j\ge 1,v_1,...,v_k\) 是所有与 \(x\) 直接相连的城市。
给出 \(q\) 个询问，每次回答 \(a_i\) 天 1 号城市的魔法值。
思路
难度陡增 倍增+矩乘优化DP
看到题目里的 \(n\leq 100\) 和这个形式，结合上面两题，你就该想到这个做法（当然还要优化）。
仔细观察，发现最终要求的总是 1 号城市（原题面里的首都）的魔法值。考虑一个 \(f_{i,0}\) 如果产生了贡献，首先 1 一定和这个点有一条路径，而且根据异或的性质，异或偶数次跟没异或一样，那么只需要考虑长度为 \(a_i\) 的路径个数的奇偶性即可。于是用矩阵快速幂求 \(a_i\) 次方即可。
但是还不够——复杂度 \(O(n^3qloga)\) ，TLE了！怎么办呢？再想想那个特殊的要求——只需要知道 1 到其他的点的情况，其余都不用考虑。那么可以预处理邻接矩阵的 \(2^k\) 次方，倍增即可。复杂度 \(O(n^3\log a+qn^2\log a)\)
代码
/*
ai的范围是2^32，直接做肯定炸掉。考虑是2的幂次，可以用倍增解决。
首先，用dis[i,j,k]表示i到j长度为2^k的路径数量。矩乘计算dis。
对于一个询问，把ai二进制拆分。f[i,j]表示前j个二进制里的1，从1到i的方案数。
*/
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N=110,M=35;
ll f[N][M],road[N][N][M],a[N],a1[M];
int n,m,q;
int main()
{
scanf( "%d%d%d",&n,&m,&q );
for ( int i=1; i<=n; i++ )
scanf( "%lld",&a[i] );
for ( int i=1,ta,b; i<=m; i++ )
scanf( "%d%d",&ta,&b ),road[ta][b][0]=road[b][ta][0]=1;
for ( int i=1; i<M; i++ )
for ( int l=1; l<=n; l++ )
for ( int j=1; j<=n; j++ )
for ( int k=1; k<=n; k++ )
road[j][k][i]+=road[j][l][i-1]*road[l][k][i-1];
for ( int i=1; i<=q; i++ )
{
ll tot=0; ll tmp,ans=0; memset( f,0,sizeof(f) );
scanf( "%lld",&tmp );
for ( int j=M-1; j>=0; j-- )
if ( (1ll<<j)<=tmp ) a1[++tot]=j,tmp-=1ll<<j;
f[1][0]=1;
for ( int j=1; j<=tot; j++ )
for ( int k=1; k<=n; k++ )
for ( int l=1; l<=n; l++ )
f[k][j]+=f[l][j-1]*road[l][k][a1[j]];
for ( int j=1; j<=n; j++ )
f[j][tot]&=1;
for ( int j=1; j<=n; j++ )
if ( f[j][tot] ) ans^=a[j];
printf( "%lld\n",ans );
}
}
7——NOI2020 Day1T1 美食家
（别问我为什么不放 AcWing 的链接，因为我这题 LOJ 上过了但是 AcWing 上 TLE了……）
题意
有 \(n\) 个城市，\(1\sim n\) 编号，\(i\) 能提供 \(c_i\) 的愉悦值。有 \(m\) 条单向道路，道路 \(i\) 起始为 \(u_i,v_i\) ，花费 \(w_i\) 天（经过之后天数+=\(w_i\)）。
计划进行为期 \(T\) 天的旅行，第 \(0\) 天从 \(1\) 出发，\(T\) 天后恰好回到 \(1\) 。每到达（包括起始点）一个城市会获得当地的愉悦值，且多次到达可以累加。中途不能停留。
有 \(k\) 个不同时间的节日，\(i\) 次在 \(t_i\) 天于 \(x_i\) 举办，如果当时在这个城市会额外得到 \(y_i\).
求愉悦值之和的最大值。
\(1\leq n\leq 50，n \leq m \leq 50,0 \leq k \leq 200,1\leq t_i \leq T \leq 10^9,1 \leq w_i \leq 5\)
思路
开始了.jpg （然而这已经是最后一题了！）
拆点，DP，矩乘，倍增优化
这类问题的解决方式（总结）
首先是描述：边有边权，问从 i 出发恰好过 k 条边到 j 的最长路
\(f[k,i,j]\) 表示上述意义，设邻接矩阵为 \(g\) ，得到转移：
\[f[k,i,j]=max_p\{f[k-1,i,p]+g_{p,j}\} \]定义一个广义矩阵乘法：
\[C_{i,j}=max_k\{A_{i,k}+B_{k,j}\} \]然后把 \(f[k]\) 看做一个矩阵，可以改写上式为：
\[f[k]=f[k-1]\times g \]根据结合律（证明略），有
\[f[k]=f[0]\times g^k \]于是得到了 \(O(n^3\log k)\) 的优良复杂度，结合倍增等方式能进一步优化。
回归题目。
首先，本题中是点权不是边权。那么把每个点的点权作为所有入边的边权，特判起始点点权。(\(f[0][1][1]=c_1\))
其次，每条边是当 \(w_i\) 条计时的，\(w_i\leq 5\) ，考虑拆点。边 \(e\) 拆分成若干点 \(e_i\) 表示在这条边上的 \(i\) 天。那么 \(e=(u,v,w)=> (u,e_1,0),...,(e_{w-2},e_{w-1},0),(e_{w-1},v,c_v)\)
最后，考虑如何处理 \(k\) 个美食节。在时间相邻的两个美食节之间，乘 \(g\) 转移，然后在美食节那天乘上一个举办城市入边边权加了 \(y\) 的转移矩阵即可。
但是这样子点数是 \(n+4m\) ，还是过不去。
考虑更换拆边为拆点，对于 \(u\to v=w\) 这个过程，看做 \(u\) 待了 \(w-1\) 天，第一天获得点权，在 \(w\) 天到达 \(v\) ，每个点 \(u\) 拆成 “待在 u 的第 i 天即可，点数 \(5n\).
朴素实现 \(O((5n)^3k\log t)\) ，由于最后只和 \(f[T,1,1]\) 有关， \(f[k]\) 保留第一行即可，再预处理 \(g^{2^k}\) 就可以通过本题。时间复杂度 \(O((5n)^3\log T+(5n)^2k\log T).\)
代码
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N=250+10,K=200+10;
int n,m,T,k,c[N],id[N][6],tot;
ll A[N];
struct node
{
int time,x,y;
bool operator < ( node tmp ) { return time<tmp.time; }
}t[K];
struct matrix
{
ll s[N][N];
matrix() { memset( s,0xc0,sizeof(s) ); }
ll * operator [] ( int x ) { return s[x]; }
matrix operator * ( matrix tmp )
{
matrix c;
for ( int i=1; i<=tot; i++ )
for ( int k=1; k<=tot; k++ )
{
if ( s[i][k]<0 ) continue;
for ( int j=1; j<=tot; j++ )
c[i][j]=max( c[i][j],s[i][k]+tmp[k][j] );
}
return c;
}
}D[35];
void mul( matrix T )
{
ll B[N];
for ( int i=1; i<=tot; i++ )
B[i]-=4e18;
for ( int k=1; k<=tot; k++ )
{
if ( A[k]<0 ) continue;
for ( int j=1; j<=tot; j++ )
B[j]=max( B[j],A[k]+T[k][j] );
}
for ( int i=1; i<=tot; i++ )
A[i]=B[i];
}
int main()
{
freopen( "delicacy.in","r",stdin ); freopen( "delicacy.out","w",stdout );
scanf( "%d%d%d%d",&n,&m,&T,&k ); tot=n;
for ( int i=1; i<=n; i++ )
id[i][0]=i;
for ( int i=1; i<=n; i++ )
scanf( "%d",&c[i] );
for ( int i=1,u,v,w; i<=m; i++ )
{
scanf( "%d%d%d",&u,&v,&w );
for ( int j=1; j<w; j++ )
if ( !id[u][j] ) id[u][j]=++tot;
for ( int j=1; j<w; j++ )
D[0][id[u][j-1]][id[u][j]]=0;
D[0][id[u][w-1]][v]=c[v];
}
for ( int i=1,t1,t2,t3; i<=k; i++ )
scanf( "%d%d%d",&t1,&t2,&t3 ),t[i]=(node){t1,t2,t3};
sort( t+1,t+1+k ); t[0]=(node){0,0,0}; t[k+1]=(node){T,0,0};
for ( int i=1; i<=30; i++ )
D[i]=D[i-1]*D[i-1];
for ( int i=2; i<=tot; i++ )
A[i]=-4e18;
A[1]=c[1];
for ( int i=1; i<=k+1; i++ )
{
int del=t[i].time-t[i-1].time;
for ( int j=30; ~j; j-- )
if ( del & (1<<j) ) mul( D[j] );
A[t[i].x]+=t[i].y;
}
if ( A[1]<0 ) printf( "-1" );
else printf( "%lld",A[1] );
fclose( stdin ); fclose( stdout );
}
8——注意事项
- 矩阵初始化，尤其是多测；最短路和普通矩乘的初始化不一样。
- 矩乘快速幂的 res 初始值，赋成单位矩阵或者递推式的初始答案矩阵
- 矩乘和Floyd不一样的！！！Floyd的 k 循环一定要在外面（枚举中间点，不然你得跑很多遍），但是矩阵就是单纯乘法，没差。
Last——鸣谢
以上的题单（前面部分）均来源于：
的附录习题，另加了一些补充题目（NOI Online和NOI2020）。
要是想系统从0开始学习也推荐这篇博文。
