[学习笔记]矩阵乘法学习笔记

1.矩阵乘法基本思想

如果 $A(n\times m)$ 数组乘上 $B(m\times q)$ 数组：

$C[i][j]=\sum _{k=1}^{m}a[i][k]\times b[k][j]$

就这个样子。。。绿色的等于两个紫色部分分别的乘积之和。。。

2.矩阵乘法的性质

不满足交换律。。。显然没法交换。。。
满足结合律。。。。。
显然满足分配律。矩阵加一下继续搞。。

3.构造矩阵解决题目

P1939 矩阵加速(数列)

考虑构造矩阵 $buf$ 解决问题，应该满足 $F(n)\times buf=F(n+1)$ 。
那么想一想 $F(n+1)$ 依赖于 $F(n)$ 的什么？因为 $F(n+1)$ 要从 $F(n-2)$ 递推而来，所以我们对于 $F(n)$ 记录三项：$F(n),F(n-1),F(n-2)$ 。
可以在这个矩阵中记录 $F(n+1)$ 依赖于 $F(n)$ 的系数。而对于 $f(n+1)$ ，我们记录 $f(n),f(n-1)$ 。即它前面两个。
第 $i$ 行第 $j$ 列代表了第 $j$ 个式子的组成，这里的系数是 $i$ 项。
上面的可能听不懂，但是看到表格一看就懂。

	$f(n+1)$	$f(n)$	$f(n-1)$
$f(n)$	1	1	0
$f(n-1)$	0	0	1
$f(n-2)$	1	0	0

矩阵 $buf$ ？就是上面的常数项。
初始值只用单位矩阵即可。

#include<bits/stdc++.h>
using namespace std;typedef int I;typedef long long LL;const int inf=0x3fffffff;int CH,FL;template<typename T>bool in(T&a){for(CH=getchar(),FL=1,a=0;!isdigit(CH)&&CH!=EOF;CH=getchar())FL=(CH=='-')?-1:1;for(;isdigit(CH)&&CH!=EOF;CH=getchar())a=a*10+CH-'0';return a*=FL,CH!=EOF;}template<typename T,typename...Args>int in(T&a,Args&...args){return in(a)+in(args...);}
const int maxn=3,maxm=3;
const int mod=1e9+7;
void pls(int&a,int b){
	a+=b;if(a>mod)a-=mod;
}
struct matrix{
	int n,m,a[maxn][maxm];
	void reset(){
		n=0;m=0;
		memset(a,0,sizeof(a));
	}
	void print(){
		for(I i=0;i<n;++i){
			for(I j=0;j<m;++j)printf("%d ",a[i][j]);
			putchar(10);
		}putchar(10);
	}
	void construct(int N,int M,int A[maxn][maxm]){
		n=N;m=M;
		for(I i=0;i<n;++i)for(I j=0;j<m;++j)a[i][j]=A[i][j];
	}
	friend matrix operator *(matrix x,matrix y){
		matrix z;z.reset();z.n=x.n;z.m=y.m;
		#define A(i,j) (x.a[i][j])
		#define B(i,j) (y.a[i][j])
		#define C(i,j) (z.a[i][j])
		for(I i=0;i<z.n;++i)
			for(I k=0,AIK=A(i,k);k<x.m;++k,AIK=A(i,k))
			for(I j=0;j<z.m;++j)
				pls(C(i,j),1ll*AIK*B(k,j)%mod);
		return z;
	}
	friend matrix operator ^(matrix a,LL b){
		int st[3][3]={{1,0,0},{0,1,0},{0,0,1}};
		matrix ans;
		ans.construct(3,3,st);
		for(;b;b>>=1){
			if(b&1)ans=ans*a;
			a=a*a;
		}
		return ans;
	}
}A;
int main(){
	int t;in(t);
	LL n;
	int gz[3][3]={{1,1,0},{0,0,1},{1,0,0}};
	while(t--){
		in(n);--n;
		A.construct(3,3,gz);
		A=A^n;
		printf("%d\n",A.a[0][0]);
	}
	return 0;
}

---

P1349 广义斐波那契数列

依然是构造矩阵。简明扼要的系数表：

	$f(n+1)$	$f(n)$
$f(n)$	$p$	$1$
$f(n-1)$	$q$	$0$

推导很简单。分析 $f(n+1)$ 的性质即可。
初始值也可以推导推导，显然得出初始值 $[b,a]$ 。。

#include<bits/stdc++.h>
using namespace std;typedef int I;typedef long long LL;const int inf=0x3fffffff;int CH,FL;template<typename T>bool in(T&a){for(CH=getchar(),FL=1,a=0;!isdigit(CH)&&CH!=EOF;CH=getchar())FL=(CH=='-')?-1:1;for(;isdigit(CH)&&CH!=EOF;CH=getchar())a=a*10+CH-'0';return a*=FL,CH!=EOF;}template<typename T,typename...Args>int in(T&a,Args&...args){return in(a)+in(args...);}
const int maxn=2,maxm=2;
int mod;
void pls(int&a,int b){
	unsigned int c=a+b;
	if(c>=mod)c-=mod;
	a=c;
}
int st,st2;
struct matrix{
	int n,m,a[maxn][maxm];
	void reset(){
		n=0;m=0;
		memset(a,0,sizeof(a));
	}
	void print(){
		for(I i=0;i<n;++i){
			for(I j=0;j<m;++j)printf("%d ",a[i][j]);
			putchar(10);
		}putchar(10);
	}
	void construct(int N,int M,int A[maxn][maxn]){
		n=N;m=M;
		for(I i=0;i<n;++i)for(I j=0;j<m;++j)a[i][j]=A[i][j];
	}
	friend matrix operator *(matrix x,matrix y){
		matrix z;z.reset();z.n=x.n;z.m=y.m;
		#define A(i,j) (x.a[i][j])
		#define B(i,j) (y.a[i][j])
		#define C(i,j) (z.a[i][j])
		for(I i=0;i<z.n;++i)
			for(I k=0,AIK=A(i,k);k<x.m;++k,AIK=A(i,k))
			for(I j=0;j<z.m;++j)
				pls(C(i,j),1ll*AIK*B(k,j)%mod);
		return z;
	}
	friend matrix operator ^(matrix a,int b){
		matrix ans;
		ans.a[0][0]=st2;
		ans.a[0][1]=st;
		ans.n=ans.m=2;
		for(;b;b>>=1){
			if(b&1)ans=ans*a;
			a=a*a;
		}
		return ans;
	}
}A;
int main(){
	int p,q,n;
	in(p,q,st,st2,n,mod);
	A.a[0][0]=p;A.a[1][0]=q;
	A.a[0][1]=1;A.a[1][1]=0;
	A.n=A.m=2;
	A=A^(n-2);
	printf("%d\n",A.a[0][0]);
	return 0;
}

---

#10222. 「一本通 6.5 例 4」佳佳的 Fibonacci

转化 interesting 。
首先 $T(n)$ 不能直接转化，但是 $T(n)=n\times S(n)-P(n)$ 。$P(n)=\sum S(i)$ 。
于是构造矩阵：

	$P[n+1]$	$S[n+1]$	$F[n+1]$	$F[n]$
$P[n]$	1	0	0	0
$S[n]$	1	1	0	0
$F[n]$	0	1	1	1
$F[n-1]$	0	0	1	0

初始值 1,1,1,0 即可。

#include<bits/stdc++.h>
using namespace std;typedef int I;typedef long long LL;const int inf=0x3fffffff;int CH,FL;template<typename T>bool in(T&a){for(CH=getchar(),FL=1,a=0;!isdigit(CH)&&CH!=EOF;CH=getchar())FL=(CH=='-')?-1:1;for(;isdigit(CH)&&CH!=EOF;CH=getchar())a=a*10+CH-'0';return a*=FL,CH!=EOF;}template<typename T,typename...Args>int in(T&a,Args&...args){return in(a)+in(args...);}
const int maxn=4,maxm=4;
int mod;
void pls(int&a,int b){
	a+=b;if(a>=mod)a-=mod;
}
struct matrix{
	int n,m,a[maxn][maxm];
	void reset(){
		n=0;m=0;
		memset(a,0,sizeof(a));
	}
	void print(){
		for(I i=0;i<n;++i){
			for(I j=0;j<m;++j)printf("%d ",a[i][j]);
			putchar(10);
		}putchar(10);
	}
	void construct(int N,int M,int A[maxn][maxn]){
		n=N;m=M;
		for(I i=0;i<n;++i)for(I j=0;j<m;++j)a[i][j]=A[i][j];
	}
	friend matrix operator *(matrix x,matrix y){
		matrix z;z.reset();z.n=x.n;z.m=y.m;
		#define A(i,j) (x.a[i][j])
		#define B(i,j) (y.a[i][j])
		#define C(i,j) (z.a[i][j])
		for(I i=0;i<z.n;++i)
			for(I k=0,AIK=A(i,k);k<x.m;++k,AIK=A(i,k))
			for(I j=0;j<z.m;++j)
				pls(C(i,j),1ll*AIK*B(k,j)%mod);
		return z;
	}
	friend matrix operator ^(matrix a,int b){
		matrix ans;
		ans.reset();
		for(I i=0;i<3;++i)ans.a[0][i]=1;
		
		ans.n=ans.m=4;
		for(;b;b>>=1){
			if(b&1)ans=ans*a;
			a=a*a;
		}
		return ans;
	}
}A;
int main(){
	int buf[4][4]={{1,0,0,0},{1,1,0,0},{0,1,1,1},{0,0,1,0}},n;
	in(n,mod);
	A.construct(4,4,buf);
	A=A^(n);
	printf("%d\n",(1ll*A.a[0][1]*n%mod-A.a[0][0]+mod+1)%mod);
	return 0;
}

1351. 阅读程序写结果

不能够直接求出，那么不妨考虑递推。

设 $f[n]$ 为答案，则 $f[n]=f[n-1]+g[n]$ 其中 $g[n]$ 为 $n\cdot p^n$
然后维护 $g[n]=g[n-1]\times p+P[n]$ 则 $P[n]=P[n-1]\ast p$

那么，构造初始矩阵：

	$f[n]$	$g[n]$	$P[n]$
$f[n-1]$	1	0	0
$g[n-1]$	p	p	0
$P[n-1]$	p	p	p

承接这个初始矩阵的初值：

$f[1]$	p
$g[1]$	p
$P[1]$	p

发现高精度，于是转十进制快速幂做法。

#include<bits/stdc++.h>
using namespace std;typedef int I;typedef long long LL;const int inf=0x3fffffff;int CH,FL;template<typename T>bool in(T&a){for(CH=getchar(),FL=1,a=0;!isdigit(CH)&&CH!=EOF;CH=getchar())FL=(CH=='-')?-1:1;for(;isdigit(CH)&&CH!=EOF;CH=getchar())a=a*10+CH-'0';return a*=FL,CH!=EOF;}template<typename T,typename...Args>int in(T&a,Args&...args){return in(a)+in(args...);}
I p;
const I maxn=3,maxm=3,mod=2009;
I NI[1010],len;
void pls(int&a,int b){
	a+=b;if(a>=mod)a-=mod;
}
struct matrix{
	int n,m,a[maxn][maxm];
	void reset(){
		n=0;m=0;
		memset(a,0,sizeof(a));
	}
	void print(){
		for(I i=0;i<n;++i){
			for(I j=0;j<m;++j)printf("%d ",a[i][j]);
			putchar(10);
		}putchar(10);
	}
	void construct(int N,int M,int A[maxn][maxm]){
		n=N;m=M;
		for(I i=0;i<n;++i)for(I j=0;j<m;++j)a[i][j]=A[i][j];
	}
	friend matrix operator *(matrix x,matrix y){
		matrix z;z.reset();z.n=x.n;z.m=y.m;
		#define A(i,j) (x.a[i][j])
		#define B(i,j) (y.a[i][j])
		#define C(i,j) (z.a[i][j])
		for(I i=0;i<z.n;++i)
			for(I k=0,AIK=A(i,k);k<x.m;++k,AIK=A(i,k))
			for(I j=0;j<z.m;++j)
				pls(C(i,j),AIK*B(k,j)%mod);
		return z;
	}friend matrix operator ^(matrix a,LL b){
		int st[3][3]={{1,0,0},{0,1,0},{0,0,1}};
		matrix ans;
		ans.construct(3,3,st);
		for(;b;b>>=1){
			if(b&1)ans=ans*a;
			a=a*a;
		}
		return ans;
	}
}A;
void pwr(matrix a){//十进制快速幂做法
	int st[3][3]={{1,0,0},{0,1,0},{0,0,1}};
	matrix ans;
	ans.construct(3,3,st);
	for(I i=1;i<=len;++i){
		ans=ans*(A^NI[i]);
		A=A^10;
	}
	A=ans;
}
int main(){
	CH=getchar();
	while(!isdigit(CH))CH=getchar();
	while(isdigit(CH))NI[++len]=CH-'0',CH=getchar();
	in(p);
	reverse(NI+1,NI+len+1);
	I a[3][3]={{1,0,0},{p,p,0},{p,p,p}};
	A.construct(3,3,a);
	pwr(A);
	printf("%d\n",A.a[2][0]);
	return 0;
}