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

2018 年 10 月 30 日 12:15

## 文章目录

$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??

### Input

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

### Sample Input

4 4
0 1
0 2
1 3
2 3
2
0 3 2
0 3 3
3 6
0 1
1 0
0 2
2 0
1 2
2 1
2
1 2 1
0 1 3
0 0


### Sample Output

2
0
1
3


### Analysis

$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) = \sum_{k=1}^{n} F(i,k,t-1) \ast G(k,j,1)$

## POJ 3613 Cow Relays

### 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

2 6 6 4
11 4 6
4 4 8
8 4 9
6 6 8
2 6 9
3 8 9


### Sample Output

10


### Analysis

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

$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]$

    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;
}
};

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(){
matrix now;now.set_infi(101,101);
for (int i=1;i<=T;i++){
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 怪兽繁殖

### Samples

#### Sample #1

Input:

1
1


Output

1


#### Sample #2

Input

1
1 1


Output

-1


#### Sample #3

Input

3
2
3
1


Output

1


#### Sample #4

Input

4
1
2 4
2
2


Output

1


#### Sample #5

Input

7
2 2
3
4 4 4
5
6
7 7 7 7
7


Output

24


#### Sample #6

Input

7
2 3
5 7
2 4
5
6
4
7


Output

5


### 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("-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("-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;
}