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

$$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$ 个点, 为了方便起见，点从 $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

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)$ 表示 $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

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)$ 表示经过两个点，$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

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$ 怪兽可以转化成 $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;
}