矩阵乘法略解

矩阵乘法

对于一个 $n\times m$ 的矩阵 $A$ 和一个 $m\times k$ 的矩阵 $B$，有且仅有一个 $n\times k$ 的矩阵 $C$ 使$$A\times B=C$$且对于 $\forall i\in [1,n]\&\forall j\in [1,k]$，有 $$C_{i,j}=\sum _{k=1}^m{A}_{ik}\times B_{kj}$$这就是矩阵乘法。比如：
$$\left[\begin{array}{ccc}1& 2& 3\\ 3& 2& 1\end{array}\right]\times \left[\begin{array}{cc}1& 1\\ 2& 2\\ 3& 3\end{array}\right]=\left[\begin{array}{cc}14& 14\\ 10& 10\end{array}\right]$$

题目描述 $\text{loj10220}$

大家都知道 Fibonacci 数列吧，$f_n=f_{n-1}+f_{n-2}$。

现在问题很简单，输入 $n$ 和 $m$，求 $f_n\mathrm{m}\mathrm{o}\mathrm{d}m$。

输入格式

输入 $n,m$。

输出格式

输出 $f_n\mathrm{m}\mathrm{o}\mathrm{d}m$。

样例输入

5 1000

样例输出

5

数据范围与提示

对于 $100\%$ 的数据，$1\le n\le 2\times 10^9,1\le m\le 10^9+10$ 。

$\text{Solution 10220}$

考虑使用矩阵快速幂优化时间复杂度。
若矩阵 $A$ 使$$\left[\begin{array}{cc}{f}_{i-2}& f_{i-1}\end{array}\right]\times A=\left[\begin{array}{cc}{f}_{i-1}& f_i\end{array}\right]$$则
$$\left[\begin{array}{cc}{f}_n& f_{n-1}\end{array}\right]=\left[\begin{array}{cc}{f}_1& f_2\end{array}\right]\times A^{n-1}$$而$A^{n-1}$可以使用快速幂求出。
下面推导 $A$。 由矩阵乘法知 $A$ 是一个 $2\times 2$ 的矩阵，且
$$f_{i-1}=A_{11}\times f_{i-2}+A_{21}\times f_{i-1}$$$\therefore A_{11}=0,A_{21}=1,$
同理得 $A_{12}=1,A_{22}=1,$
$\therefore A=\left[\begin{array}{cc}0& 1\\ 1& 1\end{array}\right].$

#include<cstdio>
#include<cstdlib>
#include<cstring>

#define reg register

typedef long long ll;

int n,m;

struct node{
	ll a[3][3];
	int x,y;
	node(){
		x=y=0;
		memset(a,0,sizeof(a));
	}
	void pt(){
		printf("%lld",a[1][2]);
	}
}s,t,ans;
node T(node a,node b){
	node c;c.x=a.x;c.y=b.y;
	for(reg int i=1;i<=c.x;++i)
		for(reg int j=1;j<=c.y;++j)
			for(reg int k=1;k<=a.y;++k)
				c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%m;
	return c;
}
node QT(node a,int b){
	if(b<=1) return a;
	node c=QT(a,b/2);
	if(b%2)
		return T(T(c,c),a);
	else
		return T(c,c);
}

int main(){
	scanf("%d%d",&n,&m);
	s.x=1;s.y=2;s.a[1][1]=s.a[1][2]=1;
	t.x=2;t.y=2;t.a[1][1]=0;t.a[1][2]=t.a[2][1]=t.a[2][2]=1;
	ans=T(s,QT(t,n-2));
	ans.pt();
}

题目描述 $\text{loj10221}$

大家都知道 Fibonacci 数列吧，$f_n=f_{n-1}+f_{n-2}$。

求$s_n={\sum }_{k=1}^n{f}_k\mathrm{m}\mathrm{o}\mathrm{d}m$的值。

样例输入

5 1000

样例输出

12

数据范围与提示

对于 $100\%$ 的数据，$1\le n\le 2\times 10^9,1\le m\le 10^9+10$ 。

$\text{Solution 10221}$

对于Fibonacci数列前 $n$ 项和 $s_n$，显然有 $s_n=s_{n-1}+f_n$。类似的，我们可以构造一矩阵 $B$，使 $\left[\begin{array}{ccc}{f}_{i-2}& f_{i-1}& s_{i-1}\end{array}\right]\times B=\left[\begin{array}{ccc}{f}_{i-1}& f_i& s_i\end{array}\right]$，就能完成此题。

#include<cstdio>
#include<cstdlib>
#include<cstring>

#define reg register

typedef long long ll;

int n,m;

struct node{
	ll a[4][4];
	int x,y;
	node(){
		x=y=0;
		memset(a,0,sizeof(a));
	}
	void pt(){
		printf("%lld",a[1][3]);
	}
}s,t,ans;
node T(node a,node b){
	node c;c.x=a.x;c.y=b.y;
	for(reg int i=1;i<=c.x;++i)
		for(reg int j=1;j<=c.y;++j)
			for(reg int k=1;k<=a.y;++k)
				c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%m;
	return c;
}
node QT(node a,int b){
	if(b<=1) return a;
	node c=QT(a,b/2);
	if(b%2)
		return T(T(c,c),a);
	else
		return T(c,c);
}

int main(){
	scanf("%d%d",&n,&m);
	s.x=1;s.y=3;s.a[1][1]=s.a[1][2]=1;s.a[1][3]=2;
	t.x=3;t.y=3;t.a[1][1]=t.a[3][1]=t.a[3][2]=0;t.a[1][2]=t.a[1][3]=t.a[2][1]=t.a[2][2]=t.a[2][3]=t.a[3][3]=1;;
	ans=T(s,QT(t,n-2));
	ans.pt();
}

题目描述 $\text{loj10222}$

大家都知道 Fibonacci 数列吧，$f_n=f_{n-1}+f_{n-2}$。

求 $T_n={\sum }_{k=1}^{n}k\cdot f_k\mathrm{m}\mathrm{o}\mathrm{d}m$ 的值。

样例输入

5 5

样例输出

4

数据范围与提示

对于 $100\%$ 的数据，$1\le n,m\le 2^{31}-1$ 。

$\text{Solution 10222}$

考虑构造一 $1\times 4$ 矩阵
$$A_i=\left[\begin{array}{cccc}{f}_{i-2}& f_{i-1}& T_{i-1}& i-1\end{array}\right]$$
使
$$A_i\times B=A_{i+1}$$
其中 $B$ 是一个 $4\times 4$ 的矩阵。同上文，我们有
$$f_{i-1}=B_{11}\cdot f_{i-2}+B_{21}\cdot f_{i-1}+B_{31}\cdot T_{i-1}+B_{41}\cdot (i-1)$$$$f_i=B_{12}\cdot f_{i-2}+B_{22}\cdot f_{i-1}+B_{32}\cdot T_{i-1}+B_{42}\cdot (i-1)$$$$T_i=\cdot \cdot \cdot$$
那么，$T_i$ 怎么用形如上式的表达式表示呢？
事实上，我们无法用形如上式的表达式表示 $T_i=T_{i-1}+i\cdot f_i$，因为我们无法用方阵 $B$ 的若干次幂的形式表示 $i$。

考虑如何将 $T_i$ 表示成一次递推式。
$$T_n=\sum _{k=1}^{n}k\cdot f_k\text{\hspace{1em}(1)}$$$$n{s}_n=\sum _{k=1}^{n}n\cdot f_k\text{\hspace{1em}(2)}$$
$(2)-(1)$ 得
$$n{s}_n-T_n=\sum _{k=1}^n(n-k)\cdot f_k\text{\hspace{1em}(3)}$$
$(3)+s_n$ 得
$$\begin{array}{rl}n{s}_n-T_n+s_n& =\sum _{k=1}^n(n-k)\cdot f_k+s_n\\ & =\sum _{k=1}^n(n-k)\cdot f_k+\sum _{k=1}^n{f}_k\\ & =\sum _{k=1}^n(n+1-k)\cdot f_k\\ & =\sum _{k=1}^{n+1}(n+1-k)\cdot f_k\\ & =(n+1)s_{n+1}-T_{n+1}\end{array}$$
设 $p_k=k{s}_k-T_k$，于是我们得到 $p$ 的递推关系式
$$p_k=p_{k-1}+s_{k-1}$$所以我们可以构造矩阵
$$A_i^{\prime }=\left[\begin{array}{cccc}{p}_i& s_i& f_i& f_{i-1}\end{array}\right]$$$$B^{\prime }=\left[\begin{array}{cccc}1& 1& 0& 0\\ 0& 1& 1& 0\\ 0& 0& 1& 1\\ 0& 0& 1& 0\end{array}\right]$$
使
$$A_i^{\prime }\times B^{\prime }=A_{i+1}^{\prime }$$

#include<cstdio>
#include<cstdlib>
#include<cstring>

#define reg register

typedef long long ll;

int n,m;

struct node{
	ll a[5][5];
	int x,y;
	node(){
		x=y=0;
		memset(a,0,sizeof(a));
	}
	void pt(){
		printf("%lld",(a[1][2]*n%m-a[1][1]+m)%m);
	}
}s,t,ans;
node T(node a,node b){
	node c;c.x=a.x;c.y=b.y;
	for(reg int i=1;i<=c.x;++i)
		for(reg int j=1;j<=c.y;++j)
			for(reg int k=1;k<=a.y;++k)
				c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%m;
	return c;
}
node QT(node a,int b){
	if(b<=1) return a;
	node c=QT(a,b/2);
	if(b%2)
		return T(T(c,c),a);
	else
		return T(c,c);
}

int main(){
	scanf("%d%d",&n,&m);
	s.x=1;s.y=4;s.a[1][2]=s.a[1][3]=s.a[1][4]=1;
	t.x=4;t.y=4;t.a[1][1]=t.a[2][1]=t.a[2][2]=t.a[3][2]=t.a[3][3]=t.a[4][3]=t.a[3][4]=1;
	ans=T(s,QT(t,n-1));
	ans.pt();
}

题目描述 $\text{loj10225}$

原题来自：SCOI 2009

Windy 在有向图中迷路了。 该有向图有 $n$ 个节点，Windy 从节点 $0$ 出发，他必须恰好在 $T$ 时刻到达节点 。

现在给出该有向图，你能告诉 Windy 总共有多少种不同的路径吗？

注意：Windy 不能在某个节点逗留，且通过某有向边的时间严格为给定的时间。

输入格式

第一行包含两个整数 $n$，$T$；
接下来有 $n$ 行，每行一个长度为 $n$ 的字符串。第 $i$ 行第 $j$ 列为 $0$ 表示从节点 $i$ 到节点 $j$ 没有边，为 $1$ 到 $9$ 表示从节点 $i$ 到节点 $j$ 需要耗费的时间。

输出格式

包含一个整数，可能的路径数，这个数可能很大，只需输出这个数除以 $2009$ 的余数。

样例输入 1

2 2
11
00

样例输出 1

1

样例输入 2

5 30
12045
07105
47805
12024
12345

样例输出 2

852

数据范围与提示

对于 $100\%$ 的数据，满足 $2\le n\le 10,1\le T\le 10^9$。

$\text{Solution 10225}$

$n$ 的规模和边权都很小，考虑拆点，则时间复杂度为 $O(n^T)$。
如何把 $T$ 优化成 $\log T$ 呢？观察输入格式，我们很自然地想到矩阵，尝试使用矩阵快速幂优化。
设从 $i$ 到 $j$ 的边权为 $u_{ij}$。考虑当边权为 $0$ 或 $1$ 时，从 $i$ 到 $j$ 路径长为 $T$ 的路径数
$$s_{ijT}=\sum _{k=1}^n(\sum u=T)$$
如果构造矩阵 U = { u i j ∣ i , j ∈ [ 1 , n ] } U=\{u_{ij}|i,j\in [1,n]\} U={uij​∣i,j∈[1,n]}$$S=U^T$$
则$$S_{ij}=s_{ijT}(\forall i,j)$$

#include<cstdio>
#include<cstdlib>
#include<cstring>

#define reg register

typedef long long ll;

int n,t;
const int m=2009;
char str[20][20];

struct node{
	ll a[110][110];
	int x,y;
	node(){
		x=y=0;
		memset(a,0,sizeof(a));
	}
}s,ans;
node T(node a,node b){
	node c;c.x=a.x;c.y=b.y;
	for(reg int i=1;i<=c.x;++i)
		for(reg int j=1;j<=c.y;++j)
			for(reg int k=1;k<=a.y;++k)
				c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%m;
	return c;
}
node QT(node a,int b){
	node c;c.x=n;c.y=n;
	for(reg int i=1;i<=n;++i)
		c.a[i][i]=1;
	while(b){
		if(b&1) c=T(c,a);
		a=T(a,a);b>>=1;
	}
	return c;
}

int main(){
	scanf("%d%d",&n,&t);
	for(reg int i=1;i<=n;++i){
		scanf("%s",str[i]+1);
		for(reg int j=1;j<=n;++j)
			str[i][j]-='0';
	}
	for(reg int i=1;i<=n;++i){
		for(reg int j=1;j<=8;++j)
			s.a[(i-1)*9+j][(i-1)*9+j+1]=1;
		for(reg int j=1;j<=n;++j)
			if(str[i][j])
				s.a[(i-1)*9+str[i][j]][(j-1)*9+1]=1;
	}
	n*=9;s.x=s.y=n;
	ans=QT(s,t);
	printf("%d",ans.a[1][n-8]);
}

题目描述 情书

求
$$\sum _{k=1}^n{f}_k+\sum _{k=1}^n{c}^k(\mathrm{m}\mathrm{o}\mathrm{d}m)$$
的值，其中 $\forall i>2$ 有
$$f_i=a\cdot f_{i-1}+b\cdot f_{i-2}$$

输入描述

一行 $7$ 个整数，分别为 $a,b,c,n,m,f_1,f_2$，相邻两数以一个空格隔开。

输出描述

一行，一个整数，表示答案。

输入样例

$\text{1 1 1 3 5 1 1}$

输出样例

$2$

数据规模与约定

对于 $100\%$ 的数据，有 $0\le a,b,c,n,m,f_i\le 2^{31}$，其中 $i\in [1,n]$。

$\text{Solution 5}$

Pt.1 求 $s_k$ 的值

仿照例2知，$a,b$ 是矩阵 $B$ 的两个元素，请读者自己尝试推导。

Pt.2 求 ${\sum }_{k=1}^n{c}^k$

简单等比数列求和。
设
$$S=\sum _{k=1}^n{c}^k\text{...(4)}$$
则
$$cS=\sum _{k=1}^n{c}^{k+1}\text{...(5)}$$
$(5)-(4)$ 得
$$(c-1)S=c^{k+1}-c$$
$$S=\frac{c^{k+1}-c}{c-1}$$
使用乘法逆元和快速幂即可完成求解。

#include<cstdio>
#include<cstdlib>
#include<cstring>

#define reg register
typedef long long ll;

struct node{
	int x,y;
	ll a[5][5];
	node(){
		memset(a,0,sizeof(a));
		x=y=0;
	}
}s,t;

ll a,b,c;
ll n,m,x,y,son;
ll sum=0;

node T(node a,node b){
	node c;c.x=a.x;c.y=b.y;
	for(reg int i=1;i<=c.x;++i)
		for(reg int j=1;j<=c.y;++j)
			for(reg int k=1;k<=a.y;++k)
				c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%m;
	return c;
}
node QT(node a,int b){
	if(b<=1) return a;
	node c=QT(a,b/2);
	if(b%2)
		return T(T(c,c),a);
	else
		return T(c,c);
}
void exgcd(ll A,ll B,ll&x,ll&y){
	if(!B){
		x=1;y=0;
		return;
	}
	exgcd(B,A%B,x,y);
	ll t=x;x=y;y=t-A/B*y;
}
ll qt(ll a,int b){
	if(b==1) return a;
	ll e=qt(a,b/2)%m;
	if(b%2)
		return (e*e)%m*a%m;
	else
		return (e*e)%m;
}
int main(){
	scanf("%lld%lld%lld%lld%lld%lld%lld",&a,&b,&c,&n,&m,&s.a[1][1],&s.a[1][2]);c%=m;
	s.a[1][3]=s.a[1][1]+s.a[1][2];s.a[1][4]=s.a[1][1];
	s.x=1;s.y=4;t.x=t.y=4;
	t.a[2][1]=t.a[3][3]=t.a[3][4]=t.a[4][4]=1;
	t.a[1][2]=t.a[1][3]=b;t.a[2][2]=t.a[2][3]=a;
	s=T(s,QT(t,n-2));
	exgcd(c-1,m,x,y);son=(qt(c,n+1)-c+m)%m;x=(x%m+m)%m;
	printf("%lld",(((s.a[1][3]*n)%m-s.a[1][4]+m)%m+(x*son)%m+m)%m);
}