矩阵乘法在图论中的简单应用

This blog post is only available in Simplified Chinse.

矩阵乘法的很多应用都是围绕矩阵乘法的定义式展开的：

$C[i,j]=\sum_{k=1}^{b} A[i,k]\ast B[k,j]$

本质上是一种动态规划吧。

代码模板

矩阵各种操作的代码模板……

struct matrix{
	static const int maxn=55;
	int n,m,a[maxn][maxn];

	matrix(){
		memset(a,0,sizeof(a));n=m=0;
	}

	void set_zero(int x,int y){
		n=x;m=y;
		memset(a,0,sizeof(a));
	}

	void set_unit(int x,int y){
		n=x;m=y;
		memset(a,0,sizeof(a));
		for (int i=0;i<min(n,m);i++) a[i][i]=1;
	}

	bool operator ==(matrix &b){
		if (n!=b.n || m!=b.m) return false;
		for (int i=0;i<n;i++){
			for (int j=0;j<n;j++){
				if (a[i][j]!=b.a[i][j]) return false;
			}
		}
		return true;
	}

	bool operator !=(matrix &b){
		return !(*this == b);
	}

	matrix operator +(matrix &b){
		matrix ret;ret.set_zero(n,m);
		for (int i=0;i<n;i++){
			for (int j=0;j<m;j++){
				ret.a[i][j]=(a[i][j]+b.a[i][j])%tt;
			}
		}
		return ret;
	}
	
	matrix operator -(matrix &b){
		matrix ret;ret.set_zero(n,m);
		for (int i=0;i<n;i++){
			for (int j=0;j<m;j++){
				ret.a[i][j]=(a[i][j]-b.a[i][j]+tt)%tt;
			}
		}
		return ret;
	}

	matrix operator *(matrix &b){
		matrix ret;ret.set_zero(n,b.m);
		for (int k=0;k<m;k++){
			for (int i=0;i<n;i++){
				for (int j=0;j<b.m;j++){
					ret.a[i][j]=(ret.a[i][j]+a[i][k]*b.a[k][j]%tt)%tt;
				}
			}
		}
		return ret;
	}

	matrix operator ^(int b){
		matrix ret;ret.set_unit(n,m);
		matrix w;w=*this;
		while (b){
			if (b&1) ret=ret*w;
			b>>=1;w=w*w;
		}
		return ret;
	}

	int count_sum(){
		int ret=0;
		for (int i=0;i<n;i++){
			for (int j=0;j<m;j++){
				ret=(ret+a[i][j])%tt;
			}
		}
		return ret;
	}
};

HDU 2157 How many ways??

HDU 2157 Link

Description

春天到了，HDU 校园里开满了花，姹紫嫣红，非常美丽。葱头是个爱花的人，看着校花校草竞相开放，漫步校园，心情也变得舒畅。为了多看看这迷人的校园，葱头决定，每次上课都走不同的路线去教室，但是由于时间问题，每次只能经过 k 个地方，比方说，这次葱头决定经过 2 个地方，那他可以先去问鼎广场看看喷泉，再去教室，也可以先到体育场跑几圈，再到教室。他非常想知道，从 A 点恰好经过 k 个点到达 B 点的方案数，当然这个数有可能非常大，所以你只要输出它模上 1000 的余数就可以了。你能帮帮他么？？你可决定了葱头一天能看多少校花哦。

Input

输入数据有多组，每组的第一行是 2 个整数 $n, m(0 < n <= 20, m <= 100)$ 表示校园内共有 $n$ 个点, 为了方便起见，点从 $0$ 到 $n-1$ 编号。接着有 $m$ 行, 每行有两个整数 $s, t (0<=s,t<n)$ 表示从 $s$ 点能到 $t$ 点，注意图是有向的。接着的一行是两个整数 $T$ ，表示有 $T$ 组询问 $(1<=T<=100)$ 。

接下来的 $T$ 行, 每行有三个整数 $A, B, k$ ，表示问你从 $A$ 点到 $B$ 点恰好经过 $k$ 个点的方案数 $(k < 20)$ ，可以走重复边。如果不存在这样的走法，则输出 0。

当 $n, m$ 都为 0 的时候输入结束。

Output

计算每次询问的方案数。由于走法很多，输出其对 1000 取模的结果。

Sample Input

Sample Output

Analysis

这题表面看起来难以着手。不妨首先考虑 $k=2$ 如何处理：我们可以定义 $F(i,j)$ 表示 $i$ 到 $j$ 经过两个点的方案数。显然可以枚举一个中间点 $k$ ， $G$ 是题中给出的邻接矩阵。那么：

$F(i,j) = \sum_{k=1}^{n} G(i,k)\ast G(k,j)$

有没有感觉这个东西很像矩阵乘法的定义式？（~~明明就是一模一样~~）

$C[i,j]=\sum_{k=1}^{b} A[i,k]\ast B[k,j]$

也就是说，经过两个点的邻接矩阵就是原来邻接矩阵的平方。那么经过三个点的邻接矩阵就是三次方吗？拓展一下，假设现在 $F(i,j,t)$ 表示 $i$ 到 $j$ 经过 $t$ 个点的方案数，那么

$F(i,j,t) = \sum_{k=1}^{n} F(i,k,t-1) \ast G(k,j,1)$

通过数学归纳法可以得到， $F(i,j,k)$ 形成的矩阵 $\bold A$ 就是初始的邻接矩阵 $\bold S$ 的 $k$ 次幂。这个问题很容易用矩阵快速幂解决。

Code

直接套模板即可↑，这里不贴了……

POJ 3613 Cow Relays

POJ 3613 Link

Description

For their physical fitness program, $N (2 ≤ N ≤ 1,000,000)$ cows have decided to run a relay race using the $T (2 ≤ T ≤ 100)$ cow trails throughout the pasture.

Each trail connects two different intersections $(1 ≤ I_{1i} ≤ 1,000; 1 ≤ I_{2i} ≤ 1,000)$ , each of which is the termination for at least two trails. The cows know the $length_i$ of each trail $(1 ≤ length_i ≤ 1,000)$ , the two intersections the trail connects, and they know that no two intersections are directly connected by two different trails. The trails form a structure known mathematically as a graph.

To run the relay, the $N$ cows position themselves at various intersections (some intersections might have more than one cow). They must position themselves properly so that they can hand off the baton cow-by-cow and end up at the proper finishing place.

Write a program to help position the cows. Find the shortest path that connects the starting intersection ( $S$ ) and the ending intersection ( $E$ ) and traverses exactly $N$ cow trails.

Input

Line 1: Four space-separated integers: $N, T, S$ , and $E$
Lines 2.. $T$ +1: Line $i$ +1 describes trail $i$ with three space-separated integers: $length_i , I_{1i}$ , and $I_{2i}$

Output

Line 1: A single integer that is the shortest distance from intersection $S$ to intersection $E$ that traverses exactly $N$ cow trails.

Sample Input

Sample Output

Translation

给你一张无向图，有 $T$ 条边，问你从起点 $S$ 到终点 $E$ ，正好经过 $N$ 条边的最短路长度。

Analysis

观察到 $N ≤ 1000000$ ，以我们（我）现在的知识水平，正常的方法肯定做不了这么大的数据范围。基本上只能考虑矩阵乘法了。

题目里虽然没有告诉你点的数量（看到那个 $l_{i1},l_{i2} \leqslant 1000$ 是不是不知所措了？！），但是注意到边的数量 $T\leqslant 100$ ，而题目又保证了每个点的度数至少为 2，也就是说实际上最多只有 100 个点。我们可以先对点进行离散化。这点很重要。

这题其实和上面一题很类似，都是正好经过 $k$ 个点（题目中保证了每个点度数至少为 2，其实也告诉了你其实就是经过 $N+1$ 个点）。受到上一题的启发，能不能和上一题一样，从简单入手进行分析呢？

考虑只经过两个点， $F(i,j)$ 表示经过两个点， $i$ 到 $j$ 的最短路，

$F(i,j)=\min(G(i,k)+G(k,j))$

考虑经过 $t$ 个点，则是

$F(i,j,t)=\min(F(i,k,t-1) + G(k,j))$

再回来看这个矩阵乘法定义式，我们似乎发现了一点小问题……

$C[i,j]=\sum_{k=1}^{b} A[i,k]\ast B[k,j]$

一个是取乘积之和，一个是取和的 min，似乎没有什么直接关联？

我们可以直接修改矩阵乘法的定义！因为矩阵只是优化本题解决方法的一种工具，想怎么定义由我们自己决定。不妨直接修改矩乘的代码。

	inline matrix operator *(matrix &b){
		matrix ret;ret.set_infi(n,b.m);
		for (int k=0;k<m;k++)
			for (int i=0;i<n;i++)
				for (int j=0;j<b.m;j++)
					ret.a[i][j]=min(ret.a[i][j],a[i][k]+b.a[k][j]);
		return ret;
	}

（需要注意一个微小的问题：用这样新的定义，在快速幂里就不能用原来定义之下的单位矩阵了，所以需要特别处理一下……）

Code

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<vector>

using namespace std;

const int INF=0x3f3f3f3f;

struct matrix{
	static const int maxn=105;
	int n,m,a[maxn][maxn];

	inline void set_zero(int x,int y){
		n=x;m=y;
		memset(a,0,sizeof(a));
	}

	inline void set_unit(int x,int y){
		n=x;m=y;
		memset(a,0,sizeof(a));
		for (int i=0;i<min(n,m);i++) a[i][i]=1;
	}

	inline void set_infi(int x,int y){
		n=x;m=y;
		memset(a,0x3f,sizeof(a));
	}

	inline matrix operator *(matrix &b){
		matrix ret;ret.set_infi(n,b.m);
		for (int k=0;k<m;k++)
			for (int i=0;i<n;i++)
				for (int j=0;j<b.m;j++)
					ret.a[i][j]=min(ret.a[i][j],a[i][k]+b.a[k][j]);
		return ret;
	}

	inline matrix operator ^(int b){
		matrix ret;ret=*this;b--;
		matrix w;w=*this;
		while (b){
			if (b&1) ret=ret*w;
			b>>=1;w=w*w;
		}
		return ret;
	}
};

inline int read(){
	int ret=0,f=1;char ch=getchar();
	while (ch<'0'||ch>'9') {if (ch=='-') f=-1;ch=getchar();}
	while (ch>='0'&&ch<='9') ret=ret*10+ch-'0',ch=getchar();
	return ret*f;
}

int N,T,S,E;
int id[1005];

signed main(){
	N=read();T=read();S=read();E=read();
	matrix now;now.set_infi(101,101);
	for (int i=1;i<=T;i++){
		int len=read(),x=read(),y=read();
		if (!id[x]) id[x]=++id[0];
		if (!id[y]) id[y]=++id[0];
		now.a[id[x]][id[y]]=now.a[id[y]][id[x]]=min(now.a[id[x]][id[y]],len);
	}
	now=now^N;
	printf("%d\n",now.a[id[S]][id[E]]);
	return 0;
}

XJOI 3879 怪兽繁殖

Description

有 N 种不同的怪兽，标号为 1 到 N。
每天晚上，所有的怪兽都会死去，当一只怪兽死去的时候，它会生出一只或者多只怪兽，因此怪兽的数量从来不会减少。

给你一个字符数组表示每种怪兽死去的时候能生出哪些种类的怪兽。假设第 i 个字符串为 2 3 3，表示第 i 只怪兽死的时候，1 只 2 类型的怪兽和 2 只 3 类型的怪兽会生出来。显然，怪兽的数量可能会是无穷大的，也有些时候会变成一个定值。

一开始你只有一只 1 类型的怪兽，你想要知道最终会有多少只怪兽。对 1e9+7 取模。如果有无穷多只，输出 -1。

Input

第一行输入一个整数 $n$ 表示怪兽的种类数 $(1≤n≤50)$ ；接下来 $n$ 行每行一个字符串由若干个数字构成，意义见题面描述。

Output

输出一个整数。

Samples

Sample #1

Input:

1
1

Output

Sample #2

Input

1
1 1

Output

-1

Sample #3

Input

Output

Sample #4

Input

Output

Sample #5

Input

Output

Sample #6

Input

Output

Hint

这题有非常多的做法，这里只介绍一种非常暴力的做法 =_= By dalao♂ → YYK 遇到这样的题目，矩阵是非常有用的武器。

Analysis

首先转化成图论问题，如果一个 $a$ 怪兽可以转化成 $k$ 个 $b$ 怪兽，那么就可以从 $a$ 向 $b$ 建一条长度为 $k$ 的边。形成的一个邻接矩阵，就是转移矩阵。那么可以可以把每个时刻怪兽的状态看成一个 1 行 n 列的矩阵，每个格子代表这种怪兽的数量。显然这样的矩阵乘以一次转移矩阵，就是进行了一天的繁殖产生的怪兽状态。这样就把问题模型进行了简化。

考虑如果最后怪兽数量进入了不变的状态，那么必然经过了多少次繁殖，数量都不变了。那么“进入不变的状态”要多少步呢？我们不妨来一番暴力的分析：直接走 $10^9$ 步，与 $2\ast 10^9$ 进行比较。如果数量一样基本就说明有环了。但是直接这样做由于最后真正的结果会非常大，有很大哈希碰撞的可能。我们干脆多模几个质数多做几次（选定一些~~国家Lingdao人~~的生日轮番作为模数）。这样就可以卡过去了……

Code

~~（当然还要加上 47 行优化……）~~

#pragma GCC target("avx")
#pragma GCC optimize(3)
#pragma GCC optimize("Ofast")
#pragma GCC optimize("inline")
#pragma GCC optimize("-fgcse")
#pragma GCC optimize("-fgcse-lm")
#pragma GCC optimize("-fipa-sra")
#pragma GCC optimize("-ftree-pre")
#pragma GCC optimize("-ftree-vrp")
#pragma GCC optimize("-fpeephole2")
#pragma GCC optimize("-ffast-math")
#pragma GCC optimize("-fsched-spec")
#pragma GCC optimize("unroll-loops")
#pragma GCC optimize("-falign-jumps")
#pragma GCC optimize("-falign-loops")
#pragma GCC optimize("-falign-labels")
#pragma GCC optimize("-fdevirtualize")
#pragma GCC optimize("-fcaller-saves")
#pragma GCC optimize("-fcrossjumping")
#pragma GCC optimize("-fthread-jumps")
#pragma GCC optimize("-funroll-loops")
#pragma GCC optimize("-fwhole-program")
#pragma GCC optimize("-freorder-blocks")
#pragma GCC optimize("-fschedule-insns")
#pragma GCC optimize("inline-functions")
#pragma GCC optimize("-ftree-tail-merge")
#pragma GCC optimize("-fschedule-insns2")
#pragma GCC optimize("-fstrict-aliasing")
#pragma GCC optimize("-fstrict-overflow")
#pragma GCC optimize("-falign-functions")
#pragma GCC optimize("-fcse-skip-blocks")
#pragma GCC optimize("-fcse-follow-jumps")
#pragma GCC optimize("-fsched-interblock")
#pragma GCC optimize("-fpartial-inlining")
#pragma GCC optimize("no-stack-protector")
#pragma GCC optimize("-freorder-functions")
#pragma GCC optimize("-findirect-inlining")
#pragma GCC optimize("-fhoist-adjacent-loads")
#pragma GCC optimize("-frerun-cse-after-loop")
#pragma GCC optimize("inline-small-functions")
#pragma GCC optimize("-finline-small-functions")
#pragma GCC optimize("-ftree-switch-conversion")
#pragma GCC optimize("-foptimize-sibling-calls")
#pragma GCC optimize("-fexpensive-optimizations")
#pragma GCC optimize("-funsafe-loop-optimizations")
#pragma GCC optimize("inline-functions-called-once")
#pragma GCC optimize("-fdelete-null-pointer-checks")

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<vector>
#include<queue>

#define int long long
using namespace std;

const int prime[4]={2333, 19260817, 19550701, 1e9+7};
int tt=1e9+7;

struct matrix{
	static const int maxn=55;
	int n,m,a[maxn][maxn];

	matrix(){
		memset(a,0,sizeof(a));n=m=0;
	}

	void set_zero(int x,int y){
		n=x;m=y;
		memset(a,0,sizeof(a));
	}

	void set_unit(int x,int y){
		n=x;m=y;
		memset(a,0,sizeof(a));
		for (int i=0;i<min(n,m);i++) a[i][i]=1;
	}

	bool operator ==(matrix &b){
		if (n!=b.n || m!=b.m) return false;
		for (int i=0;i<n;i++){
			for (int j=0;j<n;j++){
				if (a[i][j]!=b.a[i][j]) return false;
			}
		}
		return true;
	}

	bool operator !=(matrix &b){
		return !(*this == b);
	}

	matrix operator +(matrix &b){
		matrix ret;ret.set_zero(n,m);
		for (int i=0;i<n;i++){
			for (int j=0;j<m;j++){
				ret.a[i][j]=(a[i][j]+b.a[i][j])%tt;
			}
		}
		return ret;
	}
	
	matrix operator -(matrix &b){
		matrix ret;ret.set_zero(n,m);
		for (int i=0;i<n;i++){
			for (int j=0;j<m;j++){
				ret.a[i][j]=(a[i][j]-b.a[i][j]+tt)%tt;
			}
		}
		return ret;
	}

	matrix operator *(matrix &b){
		matrix ret;ret.set_zero(n,b.m);
		for (int k=0;k<m;k++){
			for (int i=0;i<n;i++){
				for (int j=0;j<b.m;j++){
					ret.a[i][j]=(ret.a[i][j]+a[i][k]*b.a[k][j]%tt)%tt;
				}
			}
		}
		return ret;
	}

	matrix operator ^(int b){
		matrix ret;ret.set_unit(n,m);
		matrix w;w=*this;
		while (b){
			if (b&1) ret=ret*w;
			b>>=1;w=w*w;
		}
		return ret;
	}

	int count_sum(){
		int ret=0;
		for (int i=0;i<n;i++){
			for (int j=0;j<m;j++){
				ret=(ret+a[i][j])%tt;
			}
		}
		return ret;
	}
};

int n;

signed main(){
	scanf("%lld",&n);
	matrix chg;chg.set_zero(n,n);
	char ch=getchar();while (ch<'0'||ch>'9') ch=getchar();
	for (int i=1;i<=n;i++){
		for (;;){
			bool flg=true;
			int now=0;
			while (ch>='0'&&ch<='9') now=now*10+ch-'0',ch=getchar();
			while (ch<'0'||ch>'9'){
				if (ch==10||ch==13||ch==-1) {flg=false;break;}
				ch=getchar();
			}
			chg.a[i-1][now-1]++;
			// printf("%lld -> %lld ++\n",i,now);
			if (flg==false && i==n) break;
			while (ch<'0'||ch>'9') ch=getchar();
			if (flg==false) break;
		}
	}

	matrix now;now.set_zero(1,n);
	now.a[0][0]=1;

	matrix result1;matrix result2;
	matrix tmp1,tmp2;
	for (int i=0;i<4;i++){
		tt=prime[i];
		tmp1=chg^(1e9);
		tmp2=chg^(2e9);
		result1=now*tmp1;
		result2=now*tmp2;
		if (result1.count_sum()!=result2.count_sum()){
			printf("-1\n");
			return 0;
		}
	}
	int ans=result1.count_sum();
	printf("%lld\n",ans);
	return 0;
}

10/30/2018

矩阵乘法在图论中的简单应用

代码模板

HDU 2157 How many ways??

Description

Input

Output

Sample Input

Sample Output

Analysis

Code

POJ 3613 Cow Relays

Description

Input

Output

Sample Input

Sample Output

Translation

Analysis

Code

XJOI 3879 怪兽繁殖

Description

Input

Output

Samples

Sample #1

Sample #2

Sample #3

Sample #4

Sample #5

Sample #6

Hint

Analysis

Code

Post a New Comment