[数据结构与算法]矩阵乘法大作战

今天给我搞的要吐了

Part 1 常规矩阵乘法

[HNOI2008]GT考试

题目描述

阿申准备报名参加 GT 考试，准考证号为 $n$ 位数 $X_1{X}_2?X_n(0\le X_i\le 9)$，他不希望准考证号上出现不吉利的数字。

他的不吉利数字 $A_1{A}_2?A_m(0\le A_i\le 9)$ 有 $m$ 位，不出现是指 $X_1{X}_2?X_n$ 中没有恰好一段等于 $A_1{A}_2?A_m$，$A_1$ 和 $X_1$ 可以为 $0$。

$1\le n\le 10^9,1\le m\le 20,2\le K\le 1000$

解题思路

先不考虑$n$的范围，那么考虑一个普通DP
设$F[i][j]$表示$i$位数，与$A$匹配到了第$j$位，的合法准考证号个数
枚举下一位填的字符$c$，求出此时匹配的长度$k$，让$F[i+1][k]+=F[i][j]$
可以用$KMP$优化匹配过程，但是$n$是$10^9$，不能一位一位转移，我们发现对于每一位来说没有本质区别，所以设$G[i]$表示长度为$i$的合法个数，构造转移矩阵，进行矩阵快速幂就行了
最后答案$Ans={\sum }_{i=0}^{m?1}F[0][i]$

#include<bits/stdc++.h>
using namespace std;
const int N = 40;
struct mat
{
	int a[N][N];
	mat()
	{
		memset(a,0,sizeof(a));
	}
}F,T;
int n,m,K;
mat mul(mat x,mat y)
{
	mat res;
	for(int i=0;i<m;i++)
	{
		for(int j=0;j<m;j++)
		{
			for(int k=0;k<m;k++)
			{
				res.a[i][j]=(res.a[i][j]+x.a[i][k]*y.a[k][j]%K)%K;
			}
		}
	}
	return res;
}
char s[N];
int Next[N];
void KMP()
{
	Next[1]=0;
	for(int i=2,j=0;i<=m;i++)
	{
		while(j&&s[j+1]!=s[i]) j=Next[j];
		if(s[j+1]==s[i]) j++;
		Next[i]=j;
	}
}
void Pow(int b)
{
	while(b)
	{
		if(b&1) F=mul(F,T);
		T=mul(T,T);
		b>>=1;
	}
}
int main()
{
	cin>>n>>m>>K;
	scanf("%s",s+1);
	KMP();
	for(int j=0;j<m;j++)
	{
		for(int ch='0';ch<='9';ch++)
		{
			int k=j;
			while(k&&s[k+1]!=ch) k=Next[k];
			if(s[k+1]==ch) k++;
			if(k<m) T.a[j][k]++;
		}
	}
	F.a[0][0]=1;
	Pow(n);
	int ans=0;
	for(int i=0;i<m;i++)
	ans=(ans+F.a[0][i])%K;
	cout<<ans;
	return 0;
}

Matrix Power Series

题目描述

给定n $\times$ n矩阵A和正整数k，找到总和 $S=A+A^2+A^3+...+A^k$。
n $\le$ 30，$k\le 10^9$，$m<10^4$

解题思路

这个$k$如此大就能让我们想到矩阵快速幂
构造一种类似分治的做法

· k %2==0

$$A+A^2+A^3+\dots +A^k=(A^{\frac{k}{2}}+1)(A+A^2+\dots +A^{\frac{k}{2}})$$
其中$1$代表单位矩阵，也就是对角线是$1$的矩阵，其他矩阵乘单位矩阵还是原矩阵

·k %2==1

$$A+A^2+A^3+\dots +A^k=(A^{\frac{k?1}{2}}+1)(A+A^2+\dots +A^{\frac{k?1}{2}})+A^k$$

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 40;
struct mat
{
	LL a[N][N];
	mat()
	{
		memset(a,0,sizeof(a));
	}
}A,E,ans;
int n,k,mod;
mat mul(mat x,mat y)
{
	mat res;
	for(int i=1;i<=n;i++)
	for(int j=1;j<=n;j++)
	for(int k=1;k<=n;k++)
	res.a[i][j]=(res.a[i][j]+x.a[i][k]*y.a[k][j]%mod)%mod;
	return res;
}
mat add(mat x,mat y)
{
	mat res;
	for(int i=1;i<=n;i++)
	for(int j=1;j<=n;j++)
	res.a[i][j]=(x.a[i][j]+y.a[i][j])%mod;
	return res;
}
mat Pow(int b)
{
	mat res,trans=A;
	for(int i=1;i<=n;i++)
	res.a[i][i]=1;
	while(b)
	{
		if(b&1) res=mul(res,trans);
		trans=mul(trans,trans);
		b>>=1;
	}
	return res;
}
void solve(int k)
{
	if(k==1)
	{
		ans=A;
	}
	else
	{
		solve(k/2);	
		ans=mul(ans,add(Pow(k/2),E));
		if(k&1) ans=add(ans,Pow(k));
	}

}
int main()
{ 
	cin>>n>>k>>mod;
	for(int i=1;i<=n;i++)
	for(int j=1;j<=n;j++)
	scanf("%lld",&A.a[i][j]);
	for(int i=1;i<=n;i++)
	E.a[i][i]=1;
	solve(k);
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
		{
			cout<<ans.a[i][j]<<' ';
		}
		cout<<endl;
	}
	return 0;
}

守望者的烦恼

题目描述

守望者-warden，长期在暗夜精灵的的首都艾萨琳内担任视察监狱的任务，监狱呈长条形，守望者warden拥有一个技能名叫“闪烁”，这个技能可以把她传送到后面的监狱内查看，她比较懒，一般不查看完所有的监狱，只是从入口进入，然后再从出口出来就算完成任务了。

头脑并不发达的warden最近在思考一个问题，她的闪烁技能是可以升级的，k级的闪烁技能最多可以向前移动k个监狱，一共有n个监狱要视察，她从入口进去，一路上有n个监狱，而且不会往回走，当然她并不用每个监狱都视察，但是她最后一定要到第n个监狱里去，因为监狱的出口在那里，但是她并不一定要到第1个监狱。
守望者warden现在想知道，她在拥有k级闪烁技能时视察n个监狱一共有多少种方案？

解题思路

很容易构造出$DP:$$F[i]$表示到第i的监狱的方案数
$$F[i]=\sum _{j=1}^{k}F[i?j]$$
转移就像K阶斐波那契数列，矩阵快速幂转移

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int K =20;
const int mod = 7777777;
struct mat
{
	LL a[K][K];
	mat()
	{
		memset(a,0,sizeof(a));
	}
}F,T;
int n,k;
mat mul(mat x,mat y)
{
	mat res;
	for(int i=1;i<=k;i++)
	for(int j=1;j<=k;j++)
	for(int l=1;l<=k;l++)
	res.a[i][j]=(res.a[i][j]+x.a[i][l]*y.a[l][j]%mod)%mod;
	return res;
}
void put(mat x)
{
	for(int i=1;i<=k;i++)
	{
		for(int j=1;j<=k;j++)
		cout<<x.a[i][j];
		cout<<endl;
	 } 
}
void Pow(int b)
{
	while(b)
	{
		if(b&1) F=mul(F,T);
		T=mul(T,T);
		b>>=1;
	}
}
int main()
{
	cin>>k>>n;
	for(int i=1;i<k;i++)
	T.a[i+1][i]=1;
	for(int i=1;i<=k;i++)
	{
		T.a[i][k]=1;
		F.a[1][k]=1;
	}
	Pow(n);
	cout<<F.a[1][k];
	return 0;
}

[SDOI2008] 递归数列

题目描述

一个由自然数组成的数列按下式定义：

$对于i<=k:a[i]=b[i]$

$对于i>k:a[i]=c[1]a[i?1]+c[2]a[i?2]+...+c[k]a[i?k]$

其中$b[j]$ 和 $c[j]（1<=j<=k）$是给定的自然数。写一个程序，给定自然数$m<=n$,

计算$a[m]+a[m+1]+a[m+2]+...+a[n]$, 并输出它除以给定自然数p的余数的值。

解题思路

这个题明示矩阵乘法，关键是构造矩阵
这里盗一下其他大佬的图片

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL N =20;
struct mat
{
	LL a[N][N];
	mat()
	{
		memset(a,0,sizeof(a));
	}
	void clear()
	{
		memset(a,0,sizeof(a));
	}
}F,T;
LL b[N],c[N],sum[N];
LL n,m,K,p;
mat mul(mat a,mat b)
{
	mat res;
	for(LL i=1;i<=K+1;i++)
	{
		for(LL j=1;j<=K+1;j++)
		{
			for(LL k=1;k<=K+1;k++)
			{
				res.a[i][j]=(res.a[i][j]+a.a[i][k]*b.a[k][j]%p)%p;
			}
		}
	}
	return res;
}
void init()
{
	T.clear();
	for(LL i=2;i<=K;i++)
	T.a[i][i-1]=1;
	for(LL i=1;i<=K;i++)
	T.a[i][K]=T.a[i][K+1]=c[K-i+1]%p;
	T.a[K+1][K+1]=1;
}
void pre()
{
	F.clear();
	for(LL i=1;i<=K;i++)
	F.a[1][i]=b[i]%p;
	F.a[1][K+1]=sum[K];	
}
void Pow(LL b)
{
	while(b)
	{
		if(b&1) F=mul(F,T);
		T=mul(T,T);
		b>>=1;
	}
}
LL ans1=0,ans2=0;
int main()
{
	freopen("spp.in","r",stdin);
	freopen("spp.out","w",stdout);
	cin>>K;
	for(LL i=1;i<=K;i++)
	scanf("%lld",&b[i]);
	for(LL i=1;i<=K;i++)
	scanf("%lld",&c[i]);
	cin>>m>>n>>p;
	for(LL i=1;i<=K;i++)
	sum[i]=(sum[i-1]+b[i]%p)%p;
	if(m>K)
	{	
		init();
		pre();
		Pow(m-1-K);
		ans1=F.a[1][K+1];
	}
	else ans1=sum[m-1];
	if(n>=K)
	{
		init();
		pre(); 
		Pow(n-K);
		ans2=F.a[1][K+1];
	}
	else ans2=sum[n];
	cout<<((ans2-ans1+p)%p+p)%p;
	return 0;
}

Part 2 矩阵乘法加速图上问题

[HNOI2002]公交车路线

题目描述

在长沙城新建的环城公路上一共有 $8$ 个公交站，分别为 A、B、C、D、E、F、G、H。公共汽车只能够在相邻的两个公交站之间运行，因此你从某一个公交站到另外一个公交站往往要换几次车，例如从公交站 A 到公交站 D，你就至少需要换 $3$ 次车。

Tiger 的方向感极其糟糕，我们知道从公交站 A 到公交 E 只需要换 $4$ 次车就可以到达，可是 tiger 却总共换了 $n$ 次车，注意 tiger 一旦到达公交站 E，他不会愚蠢到再去换车。现在希望你计算一下 tiger 有多少种可能的乘车方案。
$4\le n\le 10^7$。

解题思路

$n$如果小一点可以DP搞出来，但是现在$n$太大啦
我们要知道一个非常有用的小trick

一个图的邻接矩阵的k次方中两点的值代表经过k条边从一点到另一点的方案数

注：这里的$k$次方指矩阵对自己做矩阵乘法
所以我们构造邻接矩阵，转移，因为题目说到达E之后不能再走，所以E没有出边

#include<bits/stdc++.h>
using namespace std;
const int N = 12;
const int mod = 2009;
int a[N][N],b[9*N][9*N];
int n,m,t;
struct node
{
	int a[9*N][9*N];
	node()
	{
		memset(a,0,sizeof(a));
	}
}ans,base;
int pos[N][N];
node mul(node x,node y)
{
	node res;
	for(int i=1;i<=m;i++)
	{
		for(int j=1;j<=m;j++)
		{
			for(int k=1;k<=m;k++)
			res.a[i][j]=(res.a[i][j]+x.a[i][k]*y.a[k][j])%mod;
		}
	}
	return res;
}
void Pow(int b)
{
	while(b)
	{
		if(b&1) ans=mul(ans,base);
		base=mul(base,base);
		b=b>>1;
	}
}
char s[N+1];
int main()
{
	scanf("%d%d",&n,&t);
	for(int i=1;i<=n;i++)
	{
		scanf("%s",s+1);
		for(int j=1;j<=n;j++)
		a[i][j]=s[j]-'0';
	}
	m=n*9;
	for(int i=1;i<=n;i++)
		for(int k=0;k<9;k++)
			pos[i][k]=(i-1)*9+k+1;
	for(int i=1;i<=n;i++)
		for(int k=0;k<8;k++)
			b[pos[i][k+1]][pos[i][k]]=1;
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
		{
			if(a[i][j]!=0)
			{
				b[pos[i][0]][pos[j][a[i][j]-1]]=1;
			}
		}
	}	
	for(int i=1;i<=m;i++)
	for(int j=1;j<=m;j++)
	ans.a[i][j]=base.a[i][j]=b[i][j];

	t-=1; 
	Pow(t);	
	printf("%d",ans.a[pos[1][0]][pos[n][0]]%mod);
	return 0;
} 

P4159 [SCOI2009] 迷路

题目描述

该有向图有 $n$ 个节点，节点从 $1$ 至 $n$ 编号，windy 从节点 $1$ 出发，他必须恰好在 $t$ 时刻到达节点 $n$。

现在给出该有向图，你能告诉 windy 总共有多少种不同的路径吗？

答案对 $2009$ 取模。

注意：windy 不能在某个节点逗留，且通过某有向边的时间严格为给定的时间

解题思路

和上一道题差不多，但是边权很难受，我们考虑拆点，把每个点拆成9个点，类似分层图思想
对于一个点$u$,将其拆成$u1$，$u2$，$u3$……$u9$，并连边$u2?>u1$,$u3>u2$,……$u9?>u8$，边权是1，那么连一条权值为w的边时，连向$uw$，权值为1就行了，再用上一题的小trick
对邻接矩阵矩阵乘法就行了

#include<bits/stdc++.h>
using namespace std;
const int N = 12;
const int mod = 2009;
int a[N][N],b[9*N][9*N];
int n,m,t;
struct node
{
	int a[9*N][9*N];
	node()
	{
		memset(a,0,sizeof(a));
	}
}ans,base;
int pos[N][N];
node mul(node x,node y)
{
	node res;
	for(int i=1;i<=m;i++)
	{
		for(int j=1;j<=m;j++)
		{
			for(int k=1;k<=m;k++)
			res.a[i][j]=(res.a[i][j]+x.a[i][k]*y.a[k][j])%mod;
		}
	}
	return res;
}
void Pow(int b)
{
	while(b)
	{
		if(b&1) ans=mul(ans,base);
		base=mul(base,base);
		b=b>>1;
	}
}
char s[N+1];
int main()
{
	scanf("%d%d",&n,&t);
	for(int i=1;i<=n;i++)
	{
		scanf("%s",s+1);
		for(int j=1;j<=n;j++)
		a[i][j]=s[j]-'0';
	}
	m=n*9;
	for(int i=1;i<=n;i++)
		for(int k=0;k<9;k++)
			pos[i][k]=(i-1)*9+k+1;
	for(int i=1;i<=n;i++)
		for(int k=0;k<8;k++)
			b[pos[i][k+1]][pos[i][k]]=1;
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
		{
			if(a[i][j]!=0)
			{
				b[pos[i][0]][pos[j][a[i][j]-1]]=1;
			}
		}
	}	
	for(int i=1;i<=m;i++)
	for(int j=1;j<=m;j++)
	ans.a[i][j]=base.a[i][j]=b[i][j];

	t-=1; 
	Pow(t);	
	printf("%d",ans.a[pos[1][0]][pos[n][0]]%mod);
	return 0;
} 

[SDOI2009]HH去散步

题目描述

HH有个一成不变的习惯，喜欢饭后百步走。所谓百步走，就是散步，就是在一定的时间
内，走过一定的距离。

但是同时HH又是个喜欢变化的人，所以他不会立刻沿着刚刚走来的路走回。

又因为HH是个喜欢变化的人，所以他每天走过的路径都不完全一样，他想知道他究竟有多
少种散步的方法。

现在给你学校的地图（假设每条路的长度都是一样的都是1），问长度为t，从给定地
点A走到给定地点B共有多少条符合条件的路径。

解题思路

题目的限制很难受，所以我们要用到另一个小trick，点边互换，把无向边看成两条有向边，然后把有递进关系的两条边看成点，并在这两点之间连边，然后我们发现只要将新图中同一条无向边拆成的两条有向边对应的点的边删掉，就可以满足题目限制，强烈建议自己举例手推
然后应用小trick，套矩阵乘法

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N =200;
const int mod = 45989;
int n,m,T,s,t;
struct mat
{
	int a[N][N];
	mat()
	{
		memset(a,0,sizeof(a));
	}
	void Make(int x,int y,int v)
	{
		a[x][y]=v;
	}
	void PREE(int R)
	{
		for(int i=1;i<=R;i++)
		a[i][i]=1;
	}
}A;
int From[N],To[N],cnt=0;
mat mul(mat a,mat b)
{
	mat res;
	for(int i=1;i<=cnt;i++)
	for(int j=1;j<=cnt;j++)
	for(int k=1;k<=cnt;k++)
	res.a[i][j]=(res.a[i][j]+a.a[i][k]*b.a[k][j]%mod)%mod;
	return res;
}
mat Pow(mat x,int b)
{
	mat res;
	res.PREE(cnt);
	while(b)
	{
		if(b&1) res=mul(res,x);
		x=mul(x,x);
		b>>=1;
	}
	return res;
}
int main()
{
	freopen("walker.in","r",stdin);
	freopen("walker.out","w",stdout);
	cin>>n>>m>>T>>s>>t;
	s++;
	t++;
	From[++cnt]=0;
	To[cnt]=s;
	for(int i=1;i<=m;i++)
	{
		int x,y;
		scanf("%d%d",&x,&y);
		x++;
		y++;
		From[++cnt]=x;
		To[cnt]=y;[添加链接描述](https://www.luogu.com.cn/problem/CF576D)
		From[++cnt]=y;
		To[cnt]=x;
	}
	for(int i=1;i<=cnt;i++)
	{
		for(int j=1;j<=cnt;j++)
		{
			if(i!=j&&i!=(j^1))
			{
				if(To[i]==From[j])
				{
					A.Make(i,j,1);
				}
			}
		}
	}
	mat Ans=Pow(A,T);
	int ans=0;
	for(int i=1;i<=cnt;i++)
	{
		if(To[i]==t)
		{
			ans=(ans+Ans.a[1][i])%mod;
		}
	}
	cout<<ans;
	return 0;
}

CF576D Flights for Regular Customers

题目描述

给定一张$n$ 个点 $m$ 条边的有向图。
一开始你在 $1$ 号节点，你要走到$n$ 号节点去。
只有当你已经走过了至少 $d_i$条边时，你才能走第 $i$ 条边。
问最少要走多少条边，或判断无法到达。
$m,n\le 150,d_i\le 10^9$

解题思路

我们按照$d_i$从小到大加入每一条边，可以用矩阵乘法维护加入这条边之前走过$d_i$条边能到达的位置，然后就变成了多源bfs，也可以写成floyd，然后把邻接矩阵加入这条边
只要用ans记录最小值就行了
但是矩阵乘法是$n^3$的，而我们只需要记录能不能到达，所以我们可以用bitset加速

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int INF = 1e18+7;
const int N = 160;
int n,m;
struct mat
{
	bitset<N> a[N];
	void clear()
	{
		for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
		a[i][j]=0;
	}
	void Make(int x,int y)
	{
		a[x][y]=1;
	}
	void E()
	{
		for(int i=1;i<=n;i++)
		a[i][i]=1;
	}
}Now,G;
struct node
{
	int x,y,d;
}e[2*N];
int t=0;
void add(int x,int y,int d)
{
	e[++t].x=x;
	e[t].y=y;
	e[t].d=d;
}
bool cmp(node a,node b)
{
	return a.d<b.d;
}
mat mul(mat a,mat b)
{
	mat res;
	res.clear();
	for(int k=1;k<=n;k++)
	for(int i=1;i<=n;i++)
	if(a.a[i][k]) res.a[i]|=b.a[k];
	return res;
}
LL ans;
mat Pow(mat res,mat a,int b)
{
	while(b)
	{
		if(b&1) res=mul(res,a);
		a=mul(a,a);
		b>>=1;
	}
	return res; 
}
LL dis[N][N];
int main()
{
	freopen("flights.in","r",stdin);
	freopen("flights.out","w",stdout);
	cin>>n>>m;
	for(int i=1;i<=m;i++)
	{
		int x,y,d;
		scanf("%d%d%d",&x,&y,&d);
		add(x,y,d);
	}
	sort(e+1,e+t+1,cmp);
	ans=INF;
	Now.E();
	for(int k=1;k<=m;k++)
	{
		Now=Pow(Now,G,e[k].d-e[k-1].d);
		G.Make(e[k].x,e[k].y);
		for(int i=1;i<=n;i++)
		{
			for(int j=1;j<=n;j++)
			{
				if(G.a[i][j]==0)
				dis[i][j]=(i==j?0:INF);
				else dis[i][j]=1; 
			}
		}
		for(int l=1;l<=n;l++)
		for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
		dis[i][j]=min(dis[i][j],dis[i][l]+dis[l][j]);
		for(int i=1;i<=n;i++)
		if(Now.a[1][i]) ans=min(ans,e[k].d+dis[i][n]);
	}
	if(ans==INF) printf("Impossible");
	else cout<<ans;
	return 0;
}

P2579 [ZJOI2005]沼泽鳄鱼

题目描述

潘塔纳尔沼泽地号称世界上最大的一块湿地，它地位于巴西中部马托格罗索州的南部地区。每当雨季来临，这里碧波荡漾、生机盎然，引来不少游客。

为了让游玩更有情趣，人们在池塘的中央建设了几座石墩和石桥，每座石桥连接着两座石墩，且每两座石墩之间至多只有一座石桥。这个景点造好之后一直没敢对外开放，原因是池塘里有不少危险的食人鱼。

豆豆先生酷爱冒险，他一听说这个消息，立马赶到了池塘，想做第一个在桥上旅游的人。虽说豆豆爱冒险，但也不敢拿自己的性命开玩笑，于是他开始了仔细的实地勘察，并得到了一些惊人的结论：食人鱼的行进路线有周期性，这个周期只可能是2，3或者4个单位时间。每个单位时间里，食人鱼可以从一个石墩游到另一个石墩。每到一个石墩，如果上面有人它就会实施攻击，否则继续它的周期运动。如果没有到石墩，它是不会攻击人的。

借助先进的仪器，豆豆很快就摸清了所有食人鱼的运动规律，他要开始设计自己的行动路线了。每个单位时间里，他只可以沿着石桥从一个石墩走到另一个石墩，而不可以停在某座石墩上不动，因为站着不动还会有其它危险。如果豆豆和某条食人鱼在同一时刻到达了某座石墩，就会遭到食人鱼的袭击，他当然不希望发生这样的事情。

现在豆豆已经选好了两座石墩Start和End，他想从Start出发，经过K个单位时间后恰好站在石墩End上。假设石墩可以重复经过（包括Start和End），他想请你帮忙算算，这样的路线共有多少种（当然不能遭到食人鱼的攻击）。

解题思路

发现路线长度只有2,3，4，那么最小公倍数是12，也就是说鳄鱼可能存在的位置的局面最多只有12种，因为大于12就会循环重复，所以我们可以预处理每个局面可以走的邻接矩阵，则转移矩阵等于12个邻接矩阵的乘积，然后对于k，由循环出现的和最后剩下的几步构成，循环出现的可以用转移矩阵做矩阵快速幂，剩下的乘上对应的转移矩阵就行了

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL INF = 1e18+7;
const int N =220;
int dct[2*N],tot=0;
int segx[N],segy[N];
LL segv[N];
bool cmp(int x,int y)
{
	return x<y; 
}
struct mat
{
	LL a[N][N];
	void clear(int x)
	{
		for(int i=1;i<=x;i++)
		for(int j=1;j<=x;j++)
		a[i][j]=INF;
	}
	void Make(int x,int y,LL v)
	{
		a[x][y]=v;
	}
}F,T;
int n,t,s,e;
int m;
void put(mat p,int x)
{
	for(int i=1;i<=x;i++)
	{
		for(int j=1;j<=x;j++)
		{
			if(p.a[i][j]==INF) cout<<"INF";
			else cout<<p.a[i][j];
			cout<<' ';
		}
		cout<<endl;
	}
	cout<<"-----------"<<endl;
}
mat mul(mat x,mat y)
{
	mat res;
	res.clear(m);	
	for(int k=1;k<=m;k++)
	for(int i=1;i<=m;i++)
	for(int j=1;j<=m;j++)
	res.a[i][j]=min(res.a[i][j],x.a[i][k]+y.a[k][j]);
	return res;
}
void Pow(int b)
{
	while(b)
	{

		if(b&1) F=mul(F,T);
		T=mul(T,T);
		b>>=1;
	}
}
int id(int x)
{
	return lower_bound(dct+1,dct+m+1,x)-dct;
}
int main()
{
	freopen("relays.in","r",stdin);
	freopen("relays.out","w",stdout);
	cin>>n>>t>>s>>e;
	for(int i=1;i<=t;i++)
	{
		scanf("%d%d%lld",&segv[i],&segx[i],&segy[i]);
		dct[++tot]=segx[i];
		dct[++tot]=segy[i];
	}
	sort(dct+1,dct+tot+1,cmp);
	m=unique(dct+1,dct+tot+1)-dct-1;
	T.clear(m); 
	for(int i=1;i<=t;i++)
	{
		int x=lower_bound(dct+1,dct+m+1,segx[i])-dct;
		int y=lower_bound(dct+1,dct+m+1,segy[i])-dct;
		T.Make(x,y,segv[i]);
		T.Make(y,x,segv[i]);  
	}
	F=T;
	Pow(n-1);
	cout<<F.a[id(s)][id(e)];
	return 0;
}

[USACO07NOV]Cow Relays G

题目描述

给出一张无向连通图，求S到E经过k条边的最短路。

解题思路

设$f[i][j][k]$表示从$i$到$j$经过$k$条边的最短路
则$f[i][j][k]=min(f[i][p][k?1]+dis[p][j])$
这里的dis如果不连通为INF
如果把$k$看成阶段，滚动数组优化的话
$F[i][j]=mi{n}_{p=1}^{n}G[i][p]+dis[p][j]$
换一种字母就是
$A[i][j]=mi{n}_{p=1}^{n}B[i][p]+C[p][j]$
观察我们矩阵乘法的式子
$A[i][j]={\sum }_{p=1}^{n}B[i][p]\times C[p][j]$
发现两个式子非常像，猜想上面一种运算也具有结合律，证明略
这就是广义矩阵乘法，也可以叫成矩阵乘法的重定义，我们对其进行矩阵快速幂就行了

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL INF = 1e18+7;
const int N =220;
int dct[2*N],tot=0;
int segx[N],segy[N];
LL segv[N];
bool cmp(int x,int y)
{
	return x<y; 
}
struct mat
{
	LL a[N][N];
	void clear(int x)
	{
		for(int i=1;i<=x;i++)
		for(int j=1;j<=x;j++)
		a[i][j]=INF;
	}
	void Make(int x,int y,LL v)
	{
		a[x][y]=v;
	}
}F,T;
int n,t,s,e;
int m;
void put(mat p,int x)
{
	for(int i=1;i<=x;i++)
	{
		for(int j=1;j<=x;j++)
		{
			if(p.a[i][j]==INF) cout<<"INF";
			else cout<<p.a[i][j];
			cout<<' ';
		}
		cout<<endl;
	}
	cout<<"-----------"<<endl;
}
mat mul(mat x,mat y)
{
	mat res;
	res.clear(m);	
	for(int k=1;k<=m;k++)
	for(int i=1;i<=m;i++)
	for(int j=1;j<=m;j++)
	res.a[i][j]=min(res.a[i][j],x.a[i][k]+y.a[k][j]);
	return res;
}
void Pow(int b)
{
	while(b)
	{

		if(b&1) F=mul(F,T);
		T=mul(T,T);
		b>>=1;
	}
}
int id(int x)
{
	return lower_bound(dct+1,dct+m+1,x)-dct;
}
int main()
{
	cin>>n>>t>>s>>e;
	for(int i=1;i<=t;i++)
	{
		scanf("%d%d%lld",&segv[i],&segx[i],&segy[i]);
		dct[++tot]=segx[i];
		dct[++tot]=segy[i];
	}
	sort(dct+1,dct+tot+1,cmp);
	m=unique(dct+1,dct+tot+1)-dct-1;
	T.clear(m); 
	for(int i=1;i<=t;i++)
	{
		int x=lower_bound(dct+1,dct+m+1,segx[i])-dct;
		int y=lower_bound(dct+1,dct+m+1,segy[i])-dct;
		T.Make(x,y,segv[i]);
		T.Make(y,x,segv[i]);  
	}
	F=T;
	Pow(n-1);
	cout<<F.a[id(s)][id(e)];
	return 0;
}

P6569 [NOI Online #3 提高组] 魔法值

题目描述

H 国的交通由 $n$ 座城市与 $m$ 条道路构成，城市与道路都从 $1$ 开始编号，其中 $1$ 号城市是 $H$ 国的首都。$H$ 国中一条道路将把两个不同城市直接相连，且任意两个城市间至多有一条道路。
$H$ 国是一个信奉魔法的国家，在第 $j$ 天，$i$ 号城市的魔法值为 $f_{i,j}$ 。H 国的魔法师已观测到第 0 天时所有城市的魔法值 $f_{i,0}$ ，且他们还发现，之后的每一天每个城市的魔法值，都将会变为所有与该城市直接相连的城市的前一天魔法值的异或值，即
$$f_{x,j}=f_{y_1,j?1}?\oplus f_{y_2,j?1}?\oplus f_{y_3,j?1}?\oplus \dots \dots$$
现在 H 国的国王问了你 $q$ 个问题，对于第$i$个问题你需要回答：第 $a_i$ 天时首都的魔法值是多
少。

解题思路

像上一道题目一样定义广义矩阵乘法
$A[i][j]=XO{R}_{k=1}^{n}B[i][k]\times C[k][j]$
易证该矩阵具有结合律，进行矩阵快速幂转移
但是我们发现A、B矩阵只有一行，但是我们枚举时枚举到n，所以复杂度n^3，肯定是不行的
但是我们矩阵快速幂时依旧需要处理C矩阵的矩阵乘法，依旧是n^3
所以我们应用到一个小技巧，预处理C矩阵的2的整次幂时的形态
这样我们询问时将$a_i$ 2进制拆分，然后用A乘上对应的C的2的整次幂，这样每次就不用对C矩阵矩阵乘法了，优化掉一个n，足以通过此题

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 120;
struct mat
{
	LL a[N][N];
	int len,wid;
	mat()
	{
		len=0;
		wid=0;
		memset(a,0,sizeof(mat));
	}
}trans[40];
mat mul(mat a,mat b)
{
	mat res;
	res.len=a.len;
	res.wid=b.wid;
	for(int i=1;i<=res.len;i++)
	for(int j=1;j<=res.wid;j++)
	for(int k=1;k<=a.wid;k++)
	res.a[i][j]^=(a.a[i][k]*b.a[k][j]);
	return res;
}
int n,m,q;
LL w[N];
int main()
{
	freopen("magic.in","r",stdin);
	freopen("magic.out","w",stdout);
	cin>>n>>m>>q;
	for(int i=1;i<=n;i++)
	scanf("%lld",&w[i]);
	trans[0].len=n;
	trans[0].wid=n;
	for(int i=1;i<=m;i++)
	{
		int x,y;
		scanf("%d%d",&x,&y);
		trans[0].a[x][y]=1;
		trans[0].a[y][x]=1;
	}
	for(int i=1;i<=32;i++)
	trans[i]=mul(trans[i-1],trans[i-1]);
	while(q--)
	{
		LL k;
		scanf("%lld",&k);
		mat ans;
		for(int i=1;i<=n;i++)
		ans.a[1][i]=w[i];
		ans.len=1;
		ans.wid=n;
		for(int i=0;i<32;i++)
		{
			if((k>>i)&1) ans=mul(ans,trans[i]);
		}
		printf("%lld\n",ans.a[1][1]);
	}
	return 0;
}

Part 3 进阶题目

P6772 [NOI2020] 美食家

坐落在 Bzeroth 大陆上的精灵王国击退地灾军团的入侵后，经过十余年的休养生息，重新成为了一片欣欣向荣的乐土，吸引着八方游客。小 W 是一位游历过世界各地的著名美食家，现在也慕名来到了精灵王国。
精灵王国共有 $n$ 座城市，城市从 $1$ 到 $n$ 编号，其中城市 $i$ 的美食能为小 W 提供 $c_i$的愉悦值。精灵王国的城市通过 $m$ 条单向道路连接，道路从 $1$ 到 $m$ 编号，其中道路 $1$ 的起点为城市 $u_i$，终点为城市 $v_i$ ，沿它通行需要花费 $w_i$天。也就是说，若小 W 在第 $d$ 天从城市 $u_i$ 沿道路 $i$ 通行，那么他会在第 $d+w_i$ 天到达城市 $v_i$ ?。
小 W 计划在精灵王国进行一场为期 $T$ 天的旅行，更具体地：他会在第 $0$ 天从城市 $1$ 出发，经过 $T$ 天的旅行，最终在恰好第 $T$ 天回到城市 $1$ 结束旅行。由于小 W 是一位美食家，每当他到达一座城市时（包括第 $0$ 天和第 $T$ 天的城市 $1$），他都会品尝该城市的美食并获得其所提供的愉悦值，若小 W 多次到达同一座城市，他将获得多次愉悦值。注意旅行途中小 W 不能在任何城市停留，即当他到达一座城市且还未结束旅行时，他当天必须立即从该城市出发前往其他城市。

对于上图，小 W 一种为期 $11$ 天的可行旅游方案为 $1\to 2\to 1\to 2\to 3\to 1$
第 $0$ 天，小 W 从城市 $1$ 开始旅行，获得愉悦值 $1$ 并向城市 $2$ 出发。
第 $1$ 天，小 W 到达城市 $2$，获得愉悦值 $3$ 并向城市 $1$ 出发。
第 $4$ 天，小 W 到达城市 $1$，获得愉悦值 $1$ 并向城市 $2$ 出发。
第 $5$ 天，小 W 到达城市 $2$，获得愉悦值 $3$ 并向城市 $2$ 出发。
第 $7$ 天，小 W 到达城市 $3$，获得愉悦值 $4$ 并向城市 $1$ 出发。
第 $11$ 天，小 W 到达城市 $1$，获得愉悦值 $1$ 并结束旅行。
小 W 在该旅行中获得的愉悦值之和为 $13$。

此外，精灵王国会在不同的时间举办 $k$ 次美食节。具体来说，第 $i$ 次美食节将于第 $t_i$
天在城市 $x_i$ 举办，若小 W 第 $t_i$天时恰好在城市 $x_i$ ，那么他在品尝城市 $x_i$的美食时会额外得到 $y_i$ 的愉悦值。现在小 W 想请作为精灵王国接待使者的你帮他算出，他在旅行中能获得的愉悦值之和的最大值。

解题思路

这道题是上面那些题的大杂烩
· 首先这个数据范围就很矩阵乘法
· 先考虑没有美食节的情况
· 接着我们发现要实现矩阵乘法，边权很难受，所以类似于SDOI迷路那道题，把边拆成$w_i$条边，不过我们又发现这个图的点数小于边数，所以我们选择拆点而不选择拆边
· 然后是点权的处理，我们为了矩阵乘法方便，把点权附在连向它的边上，得到邻接矩阵C
· 然后类似 [USACO07NOV]Cow Relays G]，进行广义矩阵乘法$A[i][j]=ma{x}_{p=1}^{n}B[i][p]+C[p][j]$，进行矩阵快速幂就行了
· 再考虑有美食节的情况，先按照t排序，接着从上一个时间点通过矩乘转移到当前时间点，再在$x_i$号点的值加$y_i$
· 但是这样依然复杂度很大，所以类似 [P6569 [NOI Online #3 提高组] 魔法值]，二进制拆分然后转移，可以优化掉一维n

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL N = 300;
struct mat
{
	LL a[N][N];
	LL H,W;
	mat()
	{
		memset(a,0xc0,sizeof(a));
	}
	void clear()
	{
		memset(a,0,sizeof(a));
	}
	void Make(LL x,LL y,LL v)
	{
		a[x][y]=v;
	}
}trans[N],A;
LL n,m,T,k,tot=0;
LL c[N];
LL Note[N][50];
struct node
{
	LL t,x,v;
}F[N];
bool cmp(node a,node b)
{
	return a.t<b.t;
}
mat mul(mat a,mat b)
{
	mat res;
	res.H=a.H;
	res.W=b.W;
	for(LL i=1;i<=res.H;i++)
	{
		for(LL k=1;k<=a.W;k++)
		{
			if(a.a[i][k]<0) continue;
			for(LL j=1;j<=res.W;j++)
			res.a[i][j]=max(res.a[i][j],a.a[i][k]+b.a[k][j]);
		}
	}
	return res;
}
void pre()
{
	for(LL i=1;i<=30;i++)
	trans[i]=mul(trans[i-1],trans[i-1]);
	for(LL i=1;i<=tot;i++)
	A.a[1][i]=-1e18;
	A.a[1][1]=c[1];
	A.H=1;
	A.W=tot;
}
int main()
{
	freopen("delicacy.in","r",stdin);
	freopen("delicacy.out","w",stdout);
	scanf("%lld%lld%lld%lld",&n,&m,&T,&k);
	for(LL i=1;i<=n;i++)
	scanf("%lld",&c[i]);
	for(LL i=1;i<=n;i++)
	Note[i][0]=i;
	tot=n;
	for(LL i=1;i<=m;i++)
	{
		LL x,y,w;
		scanf("%lld%lld%lld",&x,&y,&w);
		for(LL j=1;j<w;j++)
		{
			if(!Note[x][j]) Note[x][j]=++tot;
			trans[0].Make(Note[x][j-1],Note[x][j],0);
		}
		trans[0].Make(Note[x][w-1],y,c[y]);
	}
	trans[0].W=tot;
	trans[0].H=tot;
	for(LL i=1;i<=k;i++)
	scanf("%lld%lld%lld",&F[i].t,&F[i].x,&F[i].v);
	sort(F+1,F+k+1,cmp);
	F[0].t=0;
	F[0].x=0;
	F[0].v=0;
	F[k+1].t=T;
	F[k+1].x=0;
	F[k+1].v=0;
	pre();
	for(LL i=1;i<=k+1;i++)
	{
		LL D=F[i].t-F[i-1].t;
		for (int j=0;j<=30;j++) 
		if((D>>j)&1) A=mul(A,trans[j]);
		A.a[1][F[i].x]+=F[i].v;
	}
	if(A.a[1][1]<0) printf("-1");
	else printf("%lld\n",A.a[1][1]);
	return 0;
} 

P3502 [POI2010]CHO-Hamsters

题意描述

给出 $n$ 个互不包含的字符串，要求你求出一个最短的字符串 $S$ ，使得这 $n$ 个字符串在 $S$ 中总共至少出现 $m$ 次，问 $S$ 最短是多少。

解题思路

设$F[i][j]$表示出现的第$i$ 个字符串，是$j$的方案数，枚举上一个出现的串$k$
$F[i][j]=Mi{n}_{k=1}^{n}F[i?1][k]+D[k][j]$
其中$D[k][j]$表示在$k$后面接上$j$所需要的的最少字符数，可以用kmp简单求出
上面的式子可以写成这个形式
$A[i][j]=Mi{n}_{k=1}^{n}B[i][k]+C[k][j]$
广义矩阵乘法快速幂转移即可
ps:KM不开O2会T，建议写AC自动机

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N =207;
const int S = 1e5+7;
const LL INF = 1e18;
struct mat
{
	LL a[N][N];
	void clear()
	{
		memset(a,120,sizeof(a));
	}
	void Make(int x,int y,LL v)
	{
		a[x][y]=v;
	}
}F,T;
char s[N][S];
int len[N];
int Next[N][S];
int n,m;
void Get_Next(int x)
{
	Next[x][1]=0;
	for(int i=2,j=0;i<=len[x];i++)
	{
		while(j&&s[x][j+1]!=s[x][i]) j=Next[x][j];
		if(s[x][j+1]==s[x][i]) j++;
		Next[x][i]=j;
	}
}
LL KMP(int x,int y)
{
	for(int i=2,j=0;i<=len[x];i++)
	{
		while(j&&s[y][j+1]!=s[x][i]) j=Next[y][j];
		if(s[y][j+1]==s[x][i]) j++;
		if(i==len[x]) return len[y]-j;
	}
}
mat mul(mat a,mat b)
{
	mat res;
	res.clear();
	for(int i=0;i<=n;i++)
	for(int j=0;j<=n;j++)
	for(int k=0;k<=n;k++)
	res.a[i][j]=min(res.a[i][j],a.a[i][k]+b.a[k][j]);
	return res;
}
void Pow(int b)
{
	while(b)
	{
		if(b&1) F=mul(F,T);
		T=mul(T,T);
		b>>=1;
	}
}
int main()
{
	cin>>n>>m;
	for(int i=1;i<=n;i++)
	{
		scanf("%s",s[i]+1);
		len[i]=strlen(s[i]+1);
	}
	for(int i=1;i<=n;i++)
	Get_Next(i);
	T.Make(0,0,INF);
	for(int i=1;i<=n;i++)
	{
		T.Make(0,i,len[i]);
		T.Make(i,0,INF);  
		for(int j=1;j<=n;j++)
		{
			T.Make(i,j,KMP(i,j)); 
		} 
	}
	F=T;
	Pow(m-1);
	LL ans=INF;
	for(int i=1;i<=n;i++)
	ans=min(ans,F.a[0][i]);
	cout<<ans;
	return 0;
}

压轴题

CF1067D Computer Game

题目描述

有n个游戏，每个游戏有收益$a_i$，升级后的收益$b_i$，每次成功概率$p_i$。每秒可以玩一个游戏，如果成功则得到当前收益，并且可以升级任意某个游戏。求t秒后的期望收益的最大值。n≤1e5，t≤1e10，a<b

解题思路

根据一个显然的贪心策略，我们将某个游戏升级后肯定会一直玩这个游戏，那么期望收益是$b_i\times p_i$，所以我们选择这个值最大一个游戏升级就好了
设其为W
设$f_i$表示剩下$i$时刻，且没有升级的期望收益
则$f_i=max((1?p_j)f_{i?1}+p_{j}M(i?1)+p_j{a}_j$
改式的含义为，枚举打的游戏$j$,有$1?p_j$的概率失败，剩余$i?1$时刻，有$p_j$概率成功，这一时刻收益$p_j{a}_j$，之后的收益就用W，所以是$p_{j}W(i?1)$
设$S_i=W\times i?f_i$,方程变为$f_i=max(f_{i?1}+p_j{S}_{i?1}+p_j{a}_j)$,这个式子很显然可以斜率优化，不过$t=10^{1}0$，不能直接转移，洛谷上题解能证明这个转移点具有单调性，也就是决策单调性，所以某个转移点会控制好几个区间可以用矩阵乘法优化转移过程，用二分枚举下一个转移点，更快的算法是用倍增的思想矩阵乘法，复杂度$O(nlogt)$

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef double db;
const db eps = 1e-14;
const int N =1e5+7; 
const int M = 34;
LL T;
int n;
bool equal(db a,db b)
{
	if(fabs(a-b)<eps) return 1;
	return 0;
}
struct Game
{
	int a,b;
	db p;
	bool operator < (const Game &x) const
	{
		if(equal(p,x.p)) return a*p>x.a*x.p;
		return p<x.p;
	}
}G[N];
db B;
struct Matrix
{
	int n,m;
	db w[5][5];
	Matrix()
	{
		memset(w,0,sizeof(w));
	}
}trans[N],Ans;
Matrix mul(Matrix a,Matrix b)
{
	Matrix res;
	int n=a.n,m=b.m;
	res.n=n;
	res.m=m;
	for(int i=1;i<=n;i++)
	for(int j=1;j<=m;j++)
	for(int k=1;k<=a.m;k++)
	res.w[i][j]+=a.w[i][k]*b.w[k][j];
	return res;
}
Matrix Fst(int x)
{
	Matrix res;
	res.n=3;
	res.m=3;
	res.w[1][1]=1-G[x].p;
	res.w[2][1]=G[x].p*B;
	res.w[3][1]=G[x].p*G[x].a;
	res.w[2][2]=1;
	res.w[3][2]=1;
	res.w[3][3]=1;
	return res;
}
db X(int x)
{
	return G[x].p;

}
db Y(int x)
{
	return G[x].p*G[x].a;
}
db slope(int x,int y)
{
	return (Y(y)-Y(x))/(X(y)-X(x));
}
int q[N],l=1,r=1;
db Get(db S,int x)
{
	return G[x].p*(G[x].a-S);
}
int main()
{
	cin>>n>>T;
	for(int i=1;i<=n;i++)
	{
		scanf("%d%d%lf",&G[i].a,&G[i].b,&G[i].p);
		B=max(B,G[i].b*G[i].p);
	}
	sort(G+1,G+n+1);
	for(int i=1;i<=n;i++)
	{
		if(i>1&&equal(G[i].p,G[i-1].p)) continue;
		while(l<r&&slope(q[r-1],q[r])<=slope(q[r],i)) r--;
		q[++r]=i;
	}
	Ans.n=1;
	Ans.m=3;
	Ans.w[1][3]=1;
	LL x=0;
	db k=0;
	while (x != T) 
	{
		while (l < r && Get(k, q[l]) < Get(k, q[l + 1])) l++;
		trans[0] = Fst(q[l]);
		for (int i = 1; i < M; i++) trans[i] = mul(trans[i - 1],trans[i-1]);
		for (int i = M - 1; ~i; i--) {
			if (x + (1ll << i) < T) {
				Matrix C = mul(Ans,trans[i]);
				double nK = C.w[1][1] - B * C.w[1][2];
				if (l == r || Get(nK, q[l]) > Get(nK, q[l + 1]))
					Ans = C, x += 1ll << i;
			} 
		}
		Ans = mul(Ans,trans[0]), k = Ans.w[1][1] - B * Ans.w[1][2], ++x;
	}	
	printf("%.16lf\n",Ans.w[1][1]);
	return 0;
}