矩阵乘法

定义

设 (A) 是一个 (n imes m) 的矩阵，(B) 是一个 (m imes p) 的矩阵，那么 (A imes B) 的结果 (C) 中的元素为：

[C_{i,j}=sum_{i=1}^m A_{i,k}B_{k,j} ]

运算过程就是下面这样：

矩阵快速幂

由于矩阵满足结合律和分配律，所以可以使用快速幂。

斐波那契数列

一个经典的例子是 (mathcal O(log n)) 的时间复杂度求斐波那契数列的第 (n) 项。假设第 (i) 项为 (f_i)，那么有：

[left[ egin{matrix} f_{i-1} &f_{i} end{matrix} ight] imes left[ egin{matrix} 0 & 1\ 1 & 1 end{matrix} ight] =left[egin{matrix}f_i & f_{i+1}end{matrix} ight] ]

答案就是求 (left[egin{matrix}f_{0} & f_{1}end{matrix} ight] imesleft[egin{matrix}0 & 1\1 & 1end{matrix} ight]^{n-1}=left[egin{matrix}0 & 1end{matrix} ight] imesleft[egin{matrix}0 & 1\1 & 1end{matrix} ight]^{n-1}) 的第 (1) 行第 (2) 列的数，对 (left[egin{matrix}0 & 1\1 & 1end{matrix} ight]^{n-1}) 快速幂求解即可做到 (mathcal O(log n))。

求路径数量

给定一个邻接矩阵，每经过一条边的时间都是 (1)，求 (T) 秒之后 (s) 到 (t) 的路径数量。

考虑矩阵乘法的过程 (C_{i,j}=sum_{k=1}^nA_{i,k}B_{k,j})，实际上就是在将 (isim k) 和 (ksim j) 的路径排列组合，所以直接求邻接矩阵的 (T) 次幂即可。

LOJ#10225. 迷路

边权变成了 ([1,9])，无法直接处理，那么将它拆为若干个边权为 (1) 的边。具体地，若有一条 (u) 到 (v) 边权为 (w) 的边，那么连 ((k-1)n+u o kn+u;(1le k<w)) 的边，再连一条 (u+(w-1)n o j) 的边，边权都是 (1)，就可以直接跑矩阵快速幂了。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
using namespace std;
inline int read()
{
	int x=0,f=1;char c=getchar();
	while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
	while(c>='0'&&c<='9'){x=(x<<1)+(x<<3)+(c^48);c=getchar();}
	return x*f;
}
const int N=100;
typedef long long ll;
int n;
struct mat
{
	int a[N][N];
	void init(){for(int i=1;i<=n*9;i++)a[i][i]=1;}
	mat(){memset(a,0,sizeof(a));}
	mat operator *(const mat &x)const
	{
		mat ans;
		for(int k=1;k<=n*9;k++)
			for(int i=1;i<=n*9;i++)
				for(int j=1;j<=n*9;j++)
					ans.a[i][j]+=(ll)a[i][k]*x.a[k][j]%2009,ans.a[i][j]%=2009;
		return ans;
	}
};
mat qpow(mat a,int m)
{
	mat ans;ans.init();
	while(m)
	{
		if(m&1)ans=ans*a;
		a=a*a;
		m>>=1;
	}
	return ans;
}
char s[N][N];
int main()
{
	n=read();int t=read();
	mat x;
	for(int i=1;i<=n;i++)scanf("%s",s[i]+1);
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
		{
			int cost=s[i][j]^48;
			if(!cost)continue;
			for(int k=1;k<cost;k++)x.a[i+(k-1)*n][i+k*n]=1;
			x.a[i+(cost-1)*n][j]=1;
		}
	}
	x=qpow(x,t);
	printf("%d",x.a[1][n]);
	return 0;
}