矩阵乘法略解

矩阵乘法

对于一个 \(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^{m}{k=1}{A\times B_{kj}}$$这就是矩阵乘法。比如：

\[\left[\begin{array}{ccc}1&2&3\\3&2&1\\\end{array}\right]\times\left[\begin{array}{ccc}1&1\\2&2\\3&3\\\end{array}\right]=\left[\begin{array}{ccc}14&14\\10&10\\\end{array}\right] \]

题目描述 \(\text{loj10220}\)

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

现在问题很简单，输入 \(n\) 和 \(m\)，求 \(f_n\mod m\)。

输入格式

输入 \(n,m\)。

输出格式

输出 \(f_n\mod m\)。

样例输入

5 1000

样例输出

5

数据范围与提示

对于 \(100\%\) 的数据，\(1\leq n\leq2\times10^9,1\leq m\leq 10^9+10\) 。

\(\text{Solution 10220}\)

考虑使用矩阵快速幂优化时间复杂度。
若矩阵 \(A\) 使$$\left[
\begin{array}{ccc}
f_{i-2}&f_{i-1}
\end{array}
\right] \times A=\left[
\begin{array}{ccc}
f_{i-1} & f_i
\end{array}
\right]$$则

\[\left[ \begin{array}{ccc} f_n & f_{n-1} \end{array} \right]=\left[ \begin{array}{ccc} 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}{ccc} 0 & 1\\ 1 & 1 \end{array} \right].$ ```cpp #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}$](https://loj.ac/problem/10221) 大家都知道 Fibonacci 数列吧，$f_n=f_{n-1}+f_{n-2}$。 求$s_n=\sum^{n}_{k=1}{f_k} \mod m$的值。 ## 样例输入 ``` 5 1000 ``` ## 样例输出 ``` 12 ``` ## 数据范围与提示 对于 $100\%$ 的数据，$1\leq n\leq2\times10^9,1\leq m\leq 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]$，就能完成此题。 ```cpp #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}$](https://loj.ac/problem/10222) 大家都知道 Fibonacci 数列吧，$f_n=f_{n-1}+f_{n-2}$。 求 $T_n=\sum^{n}_{k=1}{k·f_k} \mod m$ 的值。 ## 样例输入 ``` 5 5 ``` ## 样例输出 ``` 4 ``` ## 数据范围与提示 对于 $100\%$ 的数据，$1\leq n,m\leq 2^{31}-1$ 。 ## $\text{Solution 10222}$ 考虑构造一 $1\times 4$ 矩阵 $$A_i=\left[ \begin{array}{ccc} 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}·f_{i-2}+B_{21}·f_{i-1}+B_{31}·T_{i-1}+B_{41}·(i-1)$$$$f_i=B_{12}·f_{i-2}+B_{22}·f_{i-1}+B_{32}·T_{i-1}+B_{42}·(i-1)$$$$T_i=\ ··· \]

那么，\(T_i\) 怎么用形如上式的表达式表示呢？
事实上，我们无法用形如上式的表达式表示 \(T_i=T_{i-1}+i·f_i\)，因为我们无法用方阵 \(B\) 的若干次幂的形式表示 \(i\)。

考虑如何将 \(T_i\) 表示成一次递推式。

\[T_n=\sum^{n}_{k=1}{k·f_k}\tag 1$$$$ns_n=\sum^{n}_{k=1}{n·f_k}\tag 2 \]

\((2)-(1)\) 得

\[ns_n-T_n=\sum^{n}_{k=1}(n-k)·f_k\tag3 \]

\((3)+s_n\) 得

\[\begin{aligned} ns_n-T_n+s_n&=\sum^{n}_{k=1}{(n-k)·f_k}+s_n\\ &=\sum^{n}_{k=1}{(n-k)·f_k}+\sum^n_{k=1}{f_k}\\ &=\sum^{n}_{k=1}{(n+1-k)·f_k}\\ &=\sum^{n+1}_{k=1}{(n+1-k)·f_k}\\ &=(n+1)s_{n+1}-T_{n+1} \end{aligned}\]

设 \(p_k=ks_k-T_k\)，于是我们得到 \(p\) 的递推关系式

\[p_k=p_{k-1}+s_{k-1}$$所以我们可以构造矩阵 $$A_i'=\left[ \begin{array}{ccc} p_i&s_i&f_i&f_{i-1} \end{array} \right]$$$$B'=\left[ \begin{array}{ccc} 1&1&0&0\\ 0&1&1&0\\ 0&0&1&1\\ 0&0&1&0 \end{array} \right]\]

使

\[A'_i\times B'=A'_{i+1} \]

#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\leq n\leq 10,1\leq T\leq 10^9\)。

\(\text{Solution 10225}\)

\(n\) 的规模和边权都很小，考虑拆点，则时间复杂度为 \(O(n^T)\)。
如何把 \(T\) 优化成 \(\log T\) 呢？观察输入格式，我们很自然地想到矩阵，尝试使用矩阵快速幂优化。
设从 \(i\) 到 \(j\) 的边权为 \(u_{ij}\)。考虑当边权为 \(0\) 或 \(1\) 时，从 \(i\) 到 \(j\) 路径长为 \(T\) 的路径数

\[s_{ijT}=\sum^{n}_{k=1}{(\sum u=T)} \]

如果构造矩阵$$U={u_{ij}|i,j\in [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}\quad(\mod m) \]

的值，其中 \(\forall i>2\) 有

\[f_i=a·f_{i-1}+b·f_{i-2} \]

输入描述

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

输出描述

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

输入样例

\(\ \ \ \ \text{1 1 1 3 5 1 1}\)

输出样例

\(\ \ \ \ 2\)

数据规模与约定

对于 \(100\%\) 的数据，有 \(0\leq a,b,c,n,m,f_i\leq2^{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}\quad\text{...(4)} \]

则

\[cS=\sum_{k=1}^{n}{c^{k+1}}\quad\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);
}
posted @ 2019-10-12 16:29  TeacherDai  阅读(229)  评论(0)    收藏  举报