矩阵乘法优化DP复习
前言
最近做毒瘤做多了……联赛难度的东西也该复习复习了。
Warning:本文较长，难度分界线在“中场休息”部分，如果只想看普及难度的可以从第五部分直接到注意事项qwq
文中用（比如现在这个文本）引用文本书写的部分为总结性内容，即使是跳过部分也建议阅读awa
没事，最难也就NOI2020的签到题，不怕（
0——P1962 斐波那契数列
题意
$$ n\leq 2,F(n)=1. \\ n>2,F(n)=F(n-1)+F(n-2). $$
对于上述递推式，求 $F(n)\mod 1e9+7$ ，给出的 $n\leq 2^{63}$ .
思路
递推式给你了。照理讲 $O(n)$ 是多么优秀的复杂度啊 然而并不。
于是，对于一个清晰明白的递推式，考虑矩乘。
当然你要上生成函数也没人拦你，不过有一说一斐波那契的生成函数不咋好看（
首先是答案矩阵：
$$
F(n)=
\begin{pmatrix}
f(n) & f(n-1)\\
\end{pmatrix}
$$
然后是转移矩阵。转移矩阵 $base$ 可以通过：
$$
F(n)\times base=F(n+1)
$$
来确定。保险一点就是设出每个矩阵中的数，然后根据每个项解方程即可 当然也可以看着递推式直接手推 。
但是我个人比较喜欢通过矩阵意义推理。
首先，根据矩阵乘法 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$ :
$$
\begin{pmatrix}
1 & 1\\\\
1 & 0\\\\
\end{pmatrix}
$$
代码
#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 【模板】矩阵加速（数列）
题意
$$ x\leq 3,a_x=1 \\ x>3,a_x=a_{x-1}+a_{x-3} $$
求 $a_n\mod 1e9+7$
思路
推导过程同上，把 $f(n+1)$ 拆成递推式即可。
答案矩阵：
$$
F(n)=
\begin{pmatrix}
f(n),f(n-1),f(n-2)
\end{pmatrix};
$$
下一个：
$$
F(n+1)=
\begin{pmatrix}
f(n)+f(n-2),f(n),f(n-1)
\end{pmatrix}
$$
转移矩阵：
$$
base=
\begin{pmatrix}
1 & 1 & 0\\\\
0 & 0 & 1\\\\
1 & 0 & 0\\\\
\end{pmatrix}
$$
代码
#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))$
显然此题属于第三种情况。为什么能递推呢，原因在于要让最后的和最大，又是个可重集，显然每次贪心选最大和次大就好了。
答案矩阵：
$$
F(n)=
\begin{pmatrix}
f(n),f(n-1),S(n)
\end{pmatrix}
$$
下一个：
$$
F(n+1)=
\begin{pmatrix}
f(n)+f(n-1),f(n),S(n)+f(n)+f(n-1)
\end{pmatrix}
$$
转移矩阵：
$$
base=
\begin{pmatrix}
1 & 1 & 1\\\\
1 & 0 & 1\\\\
0 & 0 & 1\\\\
\end{pmatrix}
$$
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}$ ，之后每一个城市的魔法值都是：
$$ f_{x,j}=f_{v_1,j-1} \oplus f_{v_2,j-1} \oplus … \oplus f_{v_k,j-1} $$
其中 $\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开始学习也推荐这篇博文。
orz
题简单了点，但还不错
简单？！全文一个不懂！
Orz
哇，这几天也看的这个
巧了ｑｗｑ
好，很有精神！
OTZ