矩阵乘法

1. 矩阵的定义

一个 m × n m\times n m×n 的矩阵就是 m m m 行 n n n 列的数字阵列，如 2 × 3 2\times 3 2×3 的矩阵： [ 1 5 3 2 8 10 ] \left[ \begin{matrix} 1&5&3\\2&8&10\end{matrix} \right] [12​58​310​]

实际上，矩阵类似二维数组。

2. 矩阵的运算

（1）矩阵的加法和减法

矩阵的加法和减法就是将两个矩阵对应位置上的数相加减。因此，相加减的两个矩阵A，B的行列必须相同。

（2）矩阵乘法

A ， B ， C A，B，C A，B，C 是三个矩阵，若 A × B = C A\times B=C A×B=C，需要满足：

• A A A 的列数必须和 B B B 的行数相等；
• 设 A A A 是一个 n × r n\times r n×r 的矩阵， B B B 是一个 r × m r\times m r×m 的矩阵，则矩阵 A A A 乘 矩阵 B B B 的乘积 C C C 是一个 n × m n\times m n×m 的矩阵；
• C ( i , j ) = A ( i , 1 ) × B ( 1 , j ) + A ( i , 2 ) × B ( 2 , j ) + . . . A ( r , 1 ) × B ( r , j ) C_{(i,j)}=A_{(i,1)}\times B_{(1,j)}+A_{(i,2)}\times B_{(2,j)}+...A_{(r,1)}\times B_{(r,j)} C(i,j)​=A(i,1)​×B(1,j)​+A(i,2)​×B(2,j)​+...A(r,1)​×B(r,j)​

矩阵 C C C 的第 i i i 行第 j j j 列元素 = = = 矩阵 A A A 的第 i i i 行元素与矩阵 B B B的第 j j j 列对应元素乘积之和。

例如： A = [ a b c d e f ] A=\left[ \begin{matrix} a&b\\c&d\\e&f\end{matrix} \right] A=⎣⎡​ace​bdf​⎦⎤​， B = [ g h i j k l ] B=\left[ \begin{matrix} g&h&i\\j&k&l\end{matrix} \right] B=[gj​hk​il​]，求 C = A × B C=A\times B C=A×B 。

A × B = [ a × g + b × j a × h + b × k a × i + b × l c × g + d × j c × h + d × k c × i + d × l e × g + f × j e × h + f × k e × i + f × l ] A\times B=\left[ \begin{matrix} a\times g+b\times j&a\times h+b\times k&a\times i+b\times l\\c\times g+d\times j&c\times h+d\times k&c\times i+d\times l\\e\times g+f\times j&e\times h+f\times k&e\times i+f\times l\end{matrix} \right] A×B=⎣⎡​a×g+b×jc×g+d×je×g+f×j​a×h+b×kc×h+d×ke×h+f×k​a×i+b×lc×i+d×le×i+f×l​⎦⎤​

可以发现 A × B A\times B A×B 和 B × A B\times A B×A 将得到两种不同的结果，因此矩阵不满足乘法交换律。

（3）方阵次幂

如果矩阵的行和列相同，称矩阵为方阵。若 A A A 是一个方阵，方阵的幂是指，将 A A A 连乘 n n n 次， A n A^n An 。

若不是方阵则无法进行乘幂运算。

矩阵乘法不满足交换律，但是满足结合律，因此可以用快速幂的思想来求解方程次幂。

3. 矩阵乘法的应用

在信息学竞赛中考矩阵乘法不是为了考基本运算，而是利用矩阵乘法的特性，配合快速幂求解递推式的第 n n n 项，以及一些图论的构造转化。

【例1】 ( P O J   3070 ) (POJ \ 3070) (POJ 3070) F i b o n a c c i Fibonacci Fibonacci 第 n n n 项

【问题描述】

F i b o n a c c i Fibonacci Fibonacci 数列满足， f 0 = 0 f_0=0 f0​=0 ， f 1 = 1 f_1=1 f1​=1 ， f n = f n − 1 + f n − 2 （ n ≥ 2 ） f_n=f_{n-1}+f_{n-2}（n\geq 2） fn​=fn−1​+fn−2​（n≥2）。给定整数 n （ 0 ≤ n ≤ 1 0 9 ） n（0\leq n \leq 10^9） n（0≤n≤109） ，输出 f n f_n fn​ 。

【输入格式】

有多组测试数据，每组测试输出包含一个数 n n n ，当 n n n 为 − 1 -1 −1 时结束。

【输出格式】

对于每组数据，输出 f n   m o d   10000 f_n \ mod \ 10000 fn​ mod 10000 的值。

【样例输入】

0
9
999999999
1000000000
-1

【样例输出】

0
34
626
6875

看似一个简单的题目，但是当看到数据规模时，显然并不简单。这个时候，矩阵乘法就派上用场了。

可以知道：

f i = 1 × f i − 1 + 1 × f i − 2     （ 1 ） f_i=1\times f_{i-1}+1\times f_{i-2} \ \ \ （1） fi​=1×fi−1​+1×fi−2​   （1）

f i − 1 = 1 × f i − 1 + 0 × f i − 2     （ 2 ） f_{i-1}=1\times f_{i-1}+0\times f_{i-2} \ \ \ （2） fi−1​=1×fi−1​+0×fi−2​   （2）

将上式 （ 1 ） （ 2 ） （1）（2） （1）（2） 转化为矩阵运算：

[ f i f i − 1 ] = [ f i − 1 f i − 2 ] × [ 1 1 1 0 ] \left[ \begin{matrix} f_i&f_{i-1}\end{matrix} \right]=\left[ \begin{matrix} f_{i-1}&f_{i-2}\end{matrix} \right]\times\left[ \begin{matrix} 1&1\\1&0\end{matrix} \right] [fi​​fi−1​​]=[fi−1​​fi−2​​]×[11​10​]

同理，可以将 [ f i − 1 f i − 2 ] \left[ \begin{matrix} f_{i-1}&f_{i-2}\end{matrix} \right] [fi−1​​fi−2​​] 也进行转化：

[ f i f i − 1 ] = [ f i − 1 f i − 2 ] × [ 1 1 1 0 ] = [ f i − 2 f i − 3 ] × [ 1 1 1 0 ] × [ 1 1 1 0 ] \left[ \begin{matrix} f_i&f_{i-1}\end{matrix} \right]=\left[ \begin{matrix} f_{i-1}&f_{i-2}\end{matrix} \right]\times\left[ \begin{matrix} 1&1\\1&0\end{matrix} \right]=\left[ \begin{matrix} f_{i-2}&f_{i-3}\end{matrix} \right]\times\left[ \begin{matrix} 1&1\\1&0\end{matrix} \right]\times\left[ \begin{matrix} 1&1\\1&0\end{matrix} \right] [fi​​fi−1​​]=[fi−1​​fi−2​​]×[11​10​]=[fi−2​​fi−3​​]×[11​10​]×[11​10​]

所以，最后可得：

[ f n f n − 1 ] = [ f 1 f 0 ] × [ 1 1 1 0 ] n − 1 \left[ \begin{matrix} f_n&f_{n-1}\end{matrix} \right]=\left[ \begin{matrix} f_1&f_0\end{matrix} \right]\times\left[ \begin{matrix} 1&1\\1&0\end{matrix} \right]^{n-1} [fn​​fn−1​​]=[f1​​f0​​]×[11​10​]n−1

可以通过快速幂的方法求解矩阵 [ 1 1 1 0 ] n − 1 \left[ \begin{matrix} 1&1\\1&0\end{matrix} \right]^{n-1} [11​10​]n−1 ，因此求解 f n f_n fn​ 可以在 O ( 2 3 × log ⁡ N ) O(2^3\times \log N) O(23×logN) 的时间求得。

对于 2 × 2 2\times 2 2×2 的任意方阵，乘以方阵 [ 1 0 0 1 ] \left[ \begin{matrix} 1&0\\0&1\end{matrix} \right] [10​01​] 都不变，因此可以设置初始值。

可以推导出，对于 m × m m\times m m×m 的任意方阵，乘以 m × m m\times m m×m 大小的方阵 [ 1 0 . . . 0 0 1 . . . 0 . . . . . . . . . . . . 0 0 . . . 1 ] \left[ \begin{matrix} 1&0&...&0\\0&1&...&0\\...&...&...&...\\0&0&...&1\end{matrix} \right] ⎣⎢⎢⎡​10...0​01...0​............​00...1​⎦⎥⎥⎤​ 都不变，即对角线为 1 1 1 ，其余位置均为 0 0 0 。

【参考程序】

#include<bits/stdc++.h>
#define M 10000
using namespace std;
struct node
{
	int mat[3][3];
};
int n;
node z;
node matrix_mul(node x,node y)		//矩阵乘法 
{
	memset(z.mat,0,sizeof(z.mat));
	for(int i=1;i<=2;i++)
		for(int j=1;j<=2;j++)
			for(int k=1;k<=2;k++)
				z.mat[i][j]=(z.mat[i][j]+(x.mat[i][k])%M*(y.mat[k][j])%M)%M;
	return z;
}
node matrix_power(int k)			//矩阵快速幂 
{
	node t,a;
	t.mat[1][1]=1;				//初始矩阵 
	t.mat[1][2]=0;
	t.mat[2][1]=0;
	t.mat[2][2]=1;
	a.mat[1][1]=1;				//底数矩阵 
	a.mat[1][2]=1;
	a.mat[2][1]=1;
	a.mat[2][2]=0;
	while(k)
	{
		if(k&1)	t=matrix_mul(t,a);
		k>>=1;
		a=matrix_mul(a,a);
	}
	
	return t;
}
int main()
{
	while(1)
	{
		scanf("%d",&n);
		if(n==-1) break;
		if(n==0) printf("0\n");
		else	
		{
			node temp;
			temp=matrix_power(n-1);
			int ans=temp.mat[1][1]*1+temp.mat[2][1]*0; 
			printf("%d\n",ans%M);
		}
	}	
	return 0;
} 

【例2】 F i b o n a c c i Fibonacci Fibonacci 前 n n n 项和

【问题描述】

F i b o n a c c i Fibonacci Fibonacci 数列满足， f 1 = 1 f_1=1 f1​=1 ， f n = f n − 1 + f n − 2 （ n ≥ 2 ） f_n=f_{n-1}+f_{n-2}（n\geq 2） fn​=fn−1​+fn−2​（n≥2）。给定整数 n （ 0 ≤ n , m ≤ 2 × 1 0 9 ） n（0\leq n,m \leq 2\times 10^9） n（0≤n,m≤2×109） ，求前 n n n 项和模 m m m 的结果 。

【输入格式】

输入 n , m n,m n,m 。

【输出格式】

输出前 n n n 项和模 m m m 。

【样例输入】

5 1000

【样例输出】

12

【算法分析】

设 s n s_n sn​ 表示前 n n n 项和，则有：

s n = 1 × s n − 1 + 1 × f n + 0 × f n − 1 s_n=1\times s_{n-1}+1\times f_n+0\times f_{n-1} sn​=1×sn−1​+1×fn​+0×fn−1​

f n + 1 = 0 × s n − 1 + 1 × f n + 1 × f n − 1 f_{n+1}=0\times s_{n-1}+1\times f_n+1\times f_{n-1} fn+1​=0×sn−1​+1×fn​+1×fn−1​

f n = 0 × s n − 1 + 1 × f n + 0 × f n − 1 f_n=0\times s_{n-1}+1\times f_n+0\times f_{n-1} fn​=0×sn−1​+1×fn​+0×fn−1​

因此，可以构造矩阵：

[ s n f n + 1 f n ] = [ s n − 1 f n f n − 1 ] × [ 1 0 0 1 1 1 0 1 0 ] \left[ \begin{matrix} s_n&f_{n+1}&f_n\end{matrix} \right]=\left[ \begin{matrix} s_{n-1}&f_n&f_{n-1}\end{matrix} \right]\times \left[ \begin{matrix} 1&0&0\\1&1&1\\0&1&0\end{matrix} \right] [sn​​fn+1​​fn​​]=[sn−1​​fn​​fn−1​​]×⎣⎡​110​011​010​⎦⎤​

因此，可以得到：

[ s n f n + 1 f n ] = [ s 1 f 2 f 1 ] × [ 1 0 0 1 1 1 0 1 0 ] n − 1 \left[ \begin{matrix} s_n&f_{n+1}&f_n\end{matrix} \right]=\left[ \begin{matrix} s_1&f_2&f_1\end{matrix} \right]\times \left[ \begin{matrix} 1&0&0\\1&1&1\\0&1&0\end{matrix} \right]^{n-1} [sn​​fn+1​​fn​​]=[s1​​f2​​f1​​]×⎣⎡​110​011​010​⎦⎤​n−1

已知 s 1 = 1 ， f 2 = 1 ， f 1 = 1 s_1=1，f_2=1，f_1=1 s1​=1，f2​=1，f1​=1，使用矩阵快速幂可以求解 s n s_n sn​ 。

矩阵可以加速递推，这一类问题具有以下特点：

• 可以抽象为长度为 n n n 的一维矩阵，矩阵在每个单位时间变化一次。
• 变化的形式是线性递推（只有若干“加法”或“乘以一个系数”的运算）。
• 递推的轮数很大，但是矩阵长度 n n n 不大。

将长度为 n n n 的一维矩阵称为 “状态矩阵”，把用于与“状态矩阵”相乘的固定不变的矩阵称为“转移矩阵” 。

若状态矩阵中的第 x x x 个数对下一单位时间状态矩阵中的第 y y y 个数产生影响，则把转移矩阵的第 x x x 行第 y y y 列赋值为合适的系数。

矩阵乘法加速递推关键在于定义出“状态矩阵”，并根据递推式构造出正确的“转移矩阵”，之后就可以利用矩阵快速幂实现。时间复杂度为 O ( n 3 log ⁡ T ) O(n^3\log T) O(n3logT) ， T T T 为递归的总轮数。

【例3】石头游戏

【问题描述】

石头游戏的规则是这样的。石头游戏在一个 n n n 行 m m m 列（ 1 ≤ n ≤ 8 1\leq n\leq 8 1≤n≤8）的方格阵上进行。每个格子对应了一个编号在 0 ∼ 9 0\sim 9 0∼9 之间的操作序列。

操作序列是一个长度不超过 6 6 6 且循环执行、每秒执行一个字符的字符串。它包括：

• 数字 0 ∼ 9 0\sim 9 0∼9 ：拿 0 ∼ 9 0\sim 9 0∼9 个石头到该格子。
• N W S E NWSE NWSE：把这个格子内所有的石头推到相邻的格子， N N N 表示上方， W W W 表示左方， S S S 表示下方， E E E 表示右方
• D D D ：拿走这个格子的石头。

石头游戏的问题是：当这个石头游戏进行了 t t t （ 1 ≤ t ≤ 1 0 9 1\leq t\leq 10^9 1≤t≤109）秒之后，所有方格中最多的格子有多少个石头。

注意：所有格子的操作同时执行。

【输入格式】

第一行三个整数 n , m , t , a c t n, m, t, act n,m,t,act 。

接下来 n n n 行，每行 m m m 个字符，表示每个格子对应的操作序列。

最后 a c t act act （ 1 ≤ a c t ≤ 10 1\leq act \leq 10 1≤act≤10）行，每行一个字符串，表示从 0 0 0 开始的每个操作序列。

【输出格式】

一个整数：游戏进行了 t t t 秒之后，所有方格中最多的格子有多少个石头。

【样例输入】

1 6 10 3
011112
1E
E
0

【样例输出】

3

【算法分析】

以样例为例，设定一维矩阵 F t = [ a 1 a 2 a 3 a 4 a 5 a 6 ] F_t=\left[ \begin{matrix} a_1&a_2&a_3&a_4&a_5&a_6\end{matrix} \right] Ft​=[a1​​a2​​a3​​a4​​a5​​a6​​] ，表示 t t t 秒时当前每个格子的石子数量。特殊的，添加一个 a 0 a_0 a0​ ，使 a 0 a_0 a0​ 始终等于 1 1 1。

初始状态矩阵就是： F 0 = [ a 0 = 1 a 1 = 0 a 2 = 0 a 3 = 0 a 4 = 0 a 5 = 0 a 6 = 0 ] F_0=\left[ \begin{matrix} a_0=1&a_1=0&a_2=0&a_3=0&a_4=0&a_5=0&a_6=0\end{matrix} \right] F0​=[a0​=1​a1​=0​a2​=0​a3​=0​a4​=0​a5​=0​a6​=0​] ，第 1 1 1 秒， 6 6 6 个格子的操作分别是 1 , E , E , E , E , 0 1,E,E,E,E,0 1,E,E,E,E,0 ，则可以构造 7 × 7 7\times7 7×7 的转移矩阵 T 1 T_1 T1​：

保证 a 0 a_0 a0​ 始终等于1，转移矩阵 T 1 T_1 T1​ 第 0 0 0 列设定为 [ 1   0   0   0   0   0   0 ] [1\ 0\ 0\ 0\ 0\ 0\ 0] [1 0 0 0 0 0 0] ,新的状态矩阵中 a 0 = 1 × a 0 + 0 × a 1 + 0 × a 2 + 0 × a 3 + 0 × a 4 + 0 × a 5 + 0 × a 6 a_0=1\times a_0+0\times a_1+0\times a_2+0\times a_3+0\times a_4+0\times a_5+0\times a_6 a0​=1×a0​+0×a1​+0×a2​+0×a3​+0×a4​+0×a5​+0×a6​ ；

要保证放 1 1 1 个石子到第 1 1 1 个格子，同时，第 1 1 1 个格子不移动，还需要加上原本存在的石头数，则转移矩阵 T 1 T_1 T1​ 第 1 1 1 列设定为 [ 1   1   0   0   0   0   0 ] [1\ 1\ 0\ 0\ 0\ 0\ 0] [1 1 0 0 0 0 0]，新的状态矩阵中 a 1 = 1 × a 0 + 1 × a 1 + 0 × a 2 + 0 × a 3 + 0 × a 4 + 0 × a 5 + 0 × a 6 a_1=1\times a_0+1\times a_1+0\times a_2+0\times a_3+0\times a_4+0\times a_5+0\times a_6 a1​=1×a0​+1×a1​+0×a2​+0×a3​+0×a4​+0×a5​+0×a6​，实际上添加的 a 0 a_0 a0​ 就是“石头来源”，如果要放 x x x 个石子，则第 1 1 1 列设定为 [ x   1   0   0   0   0   0 ] [x \ 1\ 0\ 0\ 0\ 0 \ 0] [x 1 0 0 0 0 0] ；

将第 2 2 2 个格子的石子推向右方，则转移矩阵 T 1 T_1 T1​ 第 3 3 3 列设定为 [ 0   0   1   0   0   0   0 ] [0\ 0\ 1\ 0\ 0\ 0\ 0] [0 0 1 0 0 0 0]，新的状态矩阵中 a 3 = 0 × a 0 + 0 × a 1 + 1 × a 2 + 0 × a 3 + 0 × a 4 + 0 × a 5 + 0 × a 6 a_3=0\times a_0+0\times a_1+1\times a_2+0\times a_3+0\times a_4+0\times a_5+0\times a_6 a3​=0×a0​+0×a1​+1×a2​+0×a3​+0×a4​+0×a5​+0×a6​ ；

同理，转移矩阵 T 1 T_1 T1​ 第 4 4 4 列设定为 [ 0   0   0   1   0   0   0 ] [0\ 0\ 0\ 1\ 0\ 0\ 0] [0 0 0 1 0 0 0] ；第 5 5 5 列设定为 [ 0   0   0   0   1   0   0 ] [0\ 0\ 0\ 0\ 1\ 0\ 0] [0 0 0 0 1 0 0] ；第 6 6 6 列设定为 [ 0   0   0   0   0   1   1 ] [0\ 0\ 0\ 0\ 0\ 1\ 1] [0 0 0 0 0 1 1] ，因为第 6 6 6 个格子不移动，还要加上本身存在的。

因此： F 1 = F 0 × [ 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 ] F_1=F_0\times \left[ \begin{matrix} 1&1&0&0&0&0&0\\0&1&0&0&0&0&0\\0&0&0&1&0&0&0\\ 0&0&0&0&1&0&0\\0&0&0&0&0&1&0\\0&0&0&0&0&0&1\\0&0&0&0&0&0&1\end{matrix} \right] F1​=F0​×⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡​1000000​1100000​0000000​0010000​0001000​0000100​0000011​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤​

按照同样的方法，还可以找到 F 1 F_1 F1​ 转移到 F 2 F_2 F2​ 的转移矩阵 T 2 T_2 T2​ ， F 2 F_2 F2​ 转移到 F 3 F_3 F3​ 的转移矩阵 T 3 T_3 T3​ …

对于任意的 n n n 行 m m m 列的网格，由于数值范围较小，可以直接转化为长度为 n × m n\times m n×m 的一维矩阵 F t = [ a ( 1 , 1 ) a ( 1 , 2 ) . . . a ( 1 , m ) a ( 2 , 1 ) . . . a ( n , m ) ] F_t=\left[ \begin{matrix} a_{(1,1)}&a_{(1,2)}&...&a_{(1,m)}&a_{(2,1)}&...&a_{(n,m)}\end{matrix} \right] Ft​=[a(1,1)​​a(1,2)​​...​a(1,m)​​a(2,1)​​...​a(n,m)​​] ，其中 a ( i , j ) a_{(i,j)} a(i,j)​ 在矩阵第 ( i − 1 ) × m + j (i-1)\times m+j (i−1)×m+j 个位置，设 S ( i , j ) = ( i − 1 ) × m + j S(i,j)=(i-1)\times m+j S(i,j)=(i−1)×m+j ，同时，添加一个 a 0 a_0 a0​ ，使 a 0 a_0 a0​ 始终等于 1 1 1。

每个操作序列的长度不超过 6 6 6 ，则 1 ∼ 6 1\sim6 1∼6 的最小公倍数为 60 60 60 ，所以每经过 60 60 60 秒，操作序列又会从最开始的字符开始。因此需要构造 60 60 60 个 ( n × m + 1 ) × ( n × m + 1 ) (n\times m+1)\times (n\times m+1) (n×m+1)×(n×m+1)转移矩阵，包含第 0 ∼ ( n × m ) 0\sim (n\times m) 0∼(n×m) 行，第 0 ∼ ( n × m ) 0\sim (n\times m) 0∼(n×m) 列。

转移矩阵 T i （ 1 ≤ i ≤ 60 ） T_i（1\leq i\leq 60） Ti​（1≤i≤60） 的构造方法：

构造原理，状态矩阵 F i F_i Fi​ 每一行所有元素与转移矩阵 T i + 1 T_{i+1} Ti+1​ 第 i i i 列所有元素分别相乘的和，得到状态矩阵 F i + 1 F_{i+1} Fi+1​ 第 i i i 个元素的值。

若 F i [ S ( i , j ) ] F_i[S(i,j)] Fi​[S(i,j)] 的操作为数字 " 0 ∼ 9 " "0\sim9" "0∼9"，设数值为 x x x。则转移矩阵 T i + 1 T_{i+1} Ti+1​ 第 0 0 0 行，第 S ( i , j ) S(i,j) S(i,j) 列为 x x x ，这样保证能从石头来源获得 a 0 × x a_0\times x a0​×x 个石头。同时第 S ( i , j ) S(i,j) S(i,j) 行，第 S ( i , j ) S(i,j) S(i,j) 列为 1 1 1 ，因为没有移动，需要加上原本存在的石头，即加上 a ( i , j ) × 1 a_{(i,j)}\times1 a(i,j)​×1。

若 F i [ S ( i , j ) ] F_i[S(i,j)] Fi​[S(i,j)] 的操作为字符 " N " "N" "N" 。则转移矩阵 第 S ( i , j ) S(i,j) S(i,j) 行，第 S ( i − 1 , j ) S(i-1,j) S(i−1,j) 列为 1 1 1 ，这样计算 F i F_i Fi​ 中第 S ( i − 1 , j ) S(i-1,j) S(i−1,j) 个 元素时，会加上 a ( i , j ) × 1 a_{(i,j)}\times 1 a(i,j)​×1，相当于得到了下方的石子。字符 " W " , " S " , " E " "W","S","E" "W","S","E" 类似。

为了保证 F i [ a 0 ] F_i[a_0] Fi​[a0​] 始终为 1 1 1，则所有转移矩阵第 0 0 0 行，第 0 0 0 列一直为 1 1 1 。

需要将转移矩阵 T 1 ∼ T 60 T_1\sim T_{60} T1​∼T60​ 全部求解出来，假设 T T = T 1 × T 2 × T 3 . . . × T 60 TT=T_1\times T_2\times T_3...\times T_{60} TT=T1​×T2​×T3​...×T60​ ，则 t t t 秒后，状态矩阵 F t F_t Ft​ ：

F t = F 0 × T T t 60 × ( T 1 × T 2 × . . . × T r ) , r = t   m o d   60 F_t=F_0\times TT^{\frac t {60}}\times (T_1\times T_2\times...\times T_r),r=t \ mod \ 60 Ft​=F0​×TT60t​×(T1​×T2​×...×Tr​),r=t mod 60

其中， T T t 60 TT^{\frac t {60}} TT60t​ 可以使用矩阵加速求解，最后找出 F t F_t Ft​ 中的第 1 ∼ ( n × m ) 1\sim (n\times m) 1∼(n×m) 列的最大值即可。

【参考程序】

#include<bits/stdc++.h> 
#define M 10000
#define ll long long 
using namespace std;
struct node
{
    ll mat[70][70];
}a[65],b[65],T;
char ch[10][12];
char s[12][8];
int n,m,t,act;
int S(int x,int y)			//坐标(x,y)转化为一维矩阵的位置 
{
	return (x-1)*m+y;
}
node matrix_mul(node x,node y)      		//矩阵乘法 
{
	int nn=n*m;
	node z;
    memset(z.mat,0,sizeof(z.mat));
    for(int i=0;i<=nn;i++)
        for(int j=0;j<=nn;j++)
            for(int k=0;k<=nn;k++)
                z.mat[i][j]+=x.mat[i][k]*y.mat[k][j];
    return z;
}
void pre()		//构造前60s的转移矩阵
{
	for(int tt=1;tt<=60;tt++)		 
	{
	 	a[tt].mat[0][0]=1; 								//石头来源，保证[0][0]为1 
		for(int i=1;i<=n;i++)
			for(int j=1;j<=m;j++)		 
			{
				int u=ch[i][j]-'0';				//操作序列	 
				int v=(tt-1)%strlen(s[u]);		//t秒时，执行的操作字符为s[u][v] 
				if(s[u][v]>='0'&&s[u][v]<='9')	
				{
					a[tt].mat[0][S(i,j)]=s[u][v]-'0';	//从石头来源获取石头 
					a[tt].mat[S(i,j)][S(i,j)]=1;		//加上本身原有的个数 
				}
				if(s[u][v]=='N'&&i>1)				a[tt].mat[S(i,j)][S(i-1,j)]=1;			 
				if(s[u][v]=='W'&&j>1)				a[tt].mat[S(i,j)][S(i,j-1)]=1;
				if(s[u][v]=='S'&&i<n)				a[tt].mat[S(i,j)][S(i+1,j)]=1;
				if(s[u][v]=='E'&&j<m)				a[tt].mat[S(i,j)][S(i,j+1)]=1;
			} 		
	} 
	//b[i]表示前i秒前i个转移矩阵的乘积
	b[1]=a[1];
	for(int tt=2;tt<=60;tt++)
		b[tt]=matrix_mul(b[tt-1],a[tt]);	
} 
node matrix_power(ll k)         //方阵次幂 
{
    int nn=n*m;
 	for(int i=0;i<=nn;i++)		//初始矩阵
 		for(int j=0;j<=nn;j++)
 			if(i==j)	T.mat[i][j]=1;
 	node aa;
	aa=b[60];
    while(k)
    {
        if(k&1) T=matrix_mul(T,aa);
        k>>=1;
        aa=matrix_mul(aa,aa);
    }
    return T;
}
int main()
{
    scanf("%d%d%d%d",&n,&m,&t,&act);
	for(int i=1;i<=n;i++)
		scanf("%s",ch[i]+1);
	for(int i=0;i<act;i++)
		scanf("%s",s[i]);
	pre();		
	int k=t/60,r=t%60;
	node ans=matrix_power(k);	
	ans=matrix_mul(ans,b[r]);
	ll ans_max=0;
	for(int i=1;i<=n*m;i++)		//只需要找转移矩阵第0行最大值即可 
		ans_max=max(ans_max,ans.mat[0][i]);
	printf("%lld",ans_max);		  
    return 0;
} 

【例4】（TJOI 2017）可乐

【题目描述】

加里敦星球的人们特别喜欢喝可乐。因而，他们的敌对星球研发出了一个可乐机器人，并且放在了加里敦星球的 1 1 1 号城市上。这个可乐机器人有三种行为： 停在原地，去下一个相邻的城市，自爆。它每一秒都会随机触发一种行为。现在给加里敦星球城市图，在第 0 0 0 秒时可乐机器人在 1 1 1 号城市，问经过了 t t t 秒，可乐机器人的行为方案数是多少？

【输入格式】

第一行输入两个正整数 N N N， M M M。 N N N表示城市个数， M M M 表示道路个数。

接下来 M M M 行每行两个整数 u u u， v v v，表示 u u u， v v v 之间有一条道路。保证两座城市之间只有一条路相连，且没有任何一条道路连接两个相同的城市。

最后一行是一个整数 t t t，表示经过的时间。

【输出格式】

输出可乐机器人的行为方案数，答案可能很大，请输出对 2017 2017 2017 取模后的结果。

【样例输入】

输入 #1复制

3 2
1 2
2 3
2

【样例输出】

8

【样例解释】

• 1 1 1 ->爆炸。
• 1 1 1 -> 1 1 1 ->爆炸。
• 1 1 1 -> 2 2 2 ->爆炸。
• 1 1 1 -> 1 1 1 -> 1 1 1。
• 1 1 1 -> 1 1 1 -> 2 2 2。
• 1 1 1 -> 2 2 2 -> 1 1 1。
• 1 1 1 -> 2 2 2 -> 2 2 2。
• 1 1 1 -> 2 2 2 -> 3 3 3。

【数据范围与约定】

对于 20 % 20\% 20% 的数据，保证 t ≤ 1000 t \leq 1000 t≤1000。

对于 100 % 100\% 100%的数据，保证 1 < t ≤ 1 0 6 1 < t \leq 10^6 1<t≤106， 1 ≤ N ≤ 30 1 \leq N \leq30 1≤N≤30， 0 < M < 100 0 < M < 100 0<M<100， 1 ≤ u , v ≤ N 1 \leq u, v \leq N 1≤u,v≤N。

【算法分析】

使用邻接矩阵存储这个图，先不考虑自爆和原地停留的情况。如果点 u ， v u，v u，v 存在连边，则令 a[u][v]=a[v][u]=1。

假设存在下图的情况，右图有邻接矩阵 A A A ：

当只有一个邻接矩阵时，a[1][2]=1 可以看作 1 1 1 到 2 2 2 有 1 1 1 种方案，a[1][3] 可以看做 1 1 1 到 3 3 3 有 1 1 1 种方案，a[1][4] 可以看做 1 1 1 到 4 4 4 有 0 0 0 种方案。

两个邻接矩阵矩阵 A A A 相乘时，乘积矩阵中 a[1][1]=a[1][1]*a[1][1]+a[1][2]*a[2][1]+a[1][3]*a[3][1]+a[1][4]*a[4][1]=2 ,可以看做起点为 1 1 1 ，终点为 1 1 1 ，经过 2 2 2 步，有 2 2 2 种方案。其中 a[1][2]*a[2][1]=1 表示路径 1 − > 2 − > 1 1->2->1 1−>2−>1 ，a[1][3]*a[3][1]=1 表示路径 1 − > 3 − > 1 1->3->1 1−>3−>1 ，a[1][4]+a[4][1]=0 表示的路径 1 − > 4 − > 1 1->4->1 1−>4−>1 不存在。

同理 a[1][2]=a[1][1]*a[1][2]+a[1][2]*a[2][2]+a[1][3]*a[3][2]+a[1][4]*a[4][2]=1 ，起点为 1 1 1 ，终点为 2 2 2 ，经过 2 2 2 步，共有 1 1 1 种方案。其中 a[1][3]*a[3][2]=1 表示路径 1 − > 3 − > 2 1->3->2 1−>3−>2 。

同理可以得到a[1][3]=1,a[1][4]=1 。

因此，可以发现，邻接矩阵 A t A^t At 第 i i i 行 第 j j j 列的数字 a[i][j] 表示的就是起点为 i i i 终点为 j j j ，经过 t t t 步的方案数。

对于自爆的情况，可以添加一个结点 0 0 0 ，所有点都指向 0 0 0 ，到达 0 0 0 就表示自爆。

对于停留在原地，给每个结点添加一个自环，即 a[i][i]=1 。

因此最后答案 a n s = ∑ i = 1 n a [ 1 ] [ i ] ans=\sum_{i=1}^na[1][i] ans=∑i=1n​a[1][i] 。

【参考程序】

#include<bits/stdc++.h> 
using namespace std;
const int N=50;
int n,m,t;
struct node				//矩阵 
{
   int mat[N][N];
}s,a; 
node matrix_mul(node x,node y)      		//矩阵乘法 
{
   node z;
   memset(z.mat,0,sizeof(z.mat));
   for(int i=0;i<=n;i++)
       for(int j=0;j<=n;j++)
           for(int k=0;k<=n;k++)
            z.mat[i][j]=(z.mat[i][j]+(x.mat[i][k]*y.mat[k][j])%2017)%2017;   
   return z;
}
void matrix_power(int k)         //方阵次幂 
{
	for(int i=0;i<=n;i++)		//初始矩阵
		for(int j=0;j<=n;j++)
			if(i==j)	s.mat[i][j]=1;
   while(k)
   {
       if(k&1) s=matrix_mul(s,a);
       k>>=1;
       a=matrix_mul(a,a);
   }
}
int main()
{
   scanf("%d%d",&n,&m);
   int x,y;
   for(int i=1;i<=m;i++)
   {
   	scanf("%d%d",&x,&y);
   	a.mat[x][y]=a.mat[y][x]=1;		
   }
   for(int i=1;i<=n;i++)		
   {
   	a.mat[i][0]=1;		//设置0为自爆结点
   	a.mat[i][i]=1;		//原地停留视为自环
   }
   a.mat[0][0]=1;			//第一秒直接自爆 
   scanf("%d",&t);
   matrix_power(t); 		//求解矩阵a的t次方
   int ans=0;
   for(int i=0;i<=n;i++)
   	ans=(ans+s.mat[1][i])%2017;
   printf("%d\n",ans);
   return 0;
}