SkyWT / 博客 / 矩阵乘法在图论中的简单应用

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

2018 年 10 月 30 日 12:15

文章目录

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

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

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

代码模板

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

Github Link

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, m(0 < n <= 20, m <= 100) 表示校园内共有 nn 个点, 为了方便起见，点从 00 到 n−1n-1 编号。接着有 mm 行, 每行有两个整数 s,t(0<=s,t<n)s, t (0<=s,t<n) 表示从 ss 点能到 tt 点，注意图是有向的。接着的一行是两个整数 TT，表示有 TT 组询问 (1<=T<=100)(1<=T<=100)。

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

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

Output

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

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

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

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

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

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

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

F(i,j,t)=∑k=1nF(i,k,t−1)∗G(k,j,1)F(i,j,t) = \sum_{k=1}^{n} F(i,k,t-1) \ast G(k,j,1)

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

Code

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

POJ 3613 Cow Relays

POJ 3613 Link

Description

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

Each trail connects two different intersections (1≤I1i≤1,000;1≤I2i≤1,000)(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 lengthilength_i of each trail (1≤lengthi≤1,000)(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 NN 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 (SS) and the ending intersection (EE) and traverses exactly NN cow trails.

Input

  • Line 1: Four space-separated integers: N,T,SN, T, S , and EE
  • Lines 2..TT+1: Line ii+1 describes trail ii with three space-separated integers: lengthi,I1ilength_i , I_{1i} , and I2iI_{2i}

Output

  • Line 1: A single integer that is the shortest distance from intersection SS to intersection EE that traverses exactly NN 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

Translation

给你一张无向图，有 TT 条边，问你从起点 SS 到终点 EE，正好经过 NN 条边的最短路长度。

Analysis

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

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

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

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

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

考虑经过 tt 个点，则是

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

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

C[i,j]=∑k=1bA[i,k]∗B[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

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

Output

输出一个整数。

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

Hint

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

Analysis

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

考虑如果最后怪兽数量进入了不变的状态，那么必然经过了多少次繁殖，数量都不变了。那么“进入不变的状态”要多少步呢？我们不妨来一番暴力的分析：直接走 10910^9 步，与 2∗1092\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;
}

暂无评论

发表新的评论

所有评论都将经过博主审核。请勿填写无意义邮箱或发表无关评论、广告等，否则会被视为垃圾评论。

提交评论即表明你同意本网站使用 Cookie，并允许本站在后台记录你的邮箱、IP 地址等必要信息。这些信息不会被透露给其他用户。（提交一次评论后，本提示将不再展示）