矩阵乘法在图论中的应用

/ 0评 /

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

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

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

n, 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=2 如何处理:我们可以定义 F(i,j) 表示 ij 经过两个点的方案数。显然可以枚举一个中间点 kG 是题中给出的邻接矩阵。那么:

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) 表示 ij 经过 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 Sk 次幂。这个问题很容易用矩阵快速幂解决。

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

Output

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

给你一张无向图,有 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) 表示经过两个点,ij 的最短路,

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

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

首先转化成图论问题,如果一个 a 怪兽可以转化成 kb 怪兽,那么就可以从 ab 建一条长度为 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;
}

知识共享许可协议 本文章采用 知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议 进行许可。
欢迎转载,如有错误欢迎指出。
本文链接:https://skywt.cn/posts/matrix-multiply-in-graph/


发表评论

电子邮件地址不会被公开。 必填项已用*标注