矩阵乘法

矩阵乘法（英文：matrix multiplication）是一种根据两个矩阵得到第三个矩阵的二元运算，第三个矩阵即前两者的乘积，称为矩阵积。
它只有在第一个矩阵的列数（column）和第二个矩阵的行数（row）相同时才有意义。
矩阵乘法满足结合律，不满足交换律；
矩阵乘法的定义为 (c[i][j]=sum_{k=1}^{n}a[i][k]*b[k][j])
如图，(C矩阵=A矩阵*B矩阵)
那么(c的行数c.n=a的行数a.n，c的列数c.m=b的列数b.m).，(a.m=b.n)
也只有这样时，(a[i][k]和b[k][i])才有意义

矩阵乘法的用处是配合快速幂在(log)的时间复杂度里计算有递推式的数列第n项。
就比如斐波拉锲数列的递推式是(f[n]=f[n-1]+f[n-2])
那就可以用如下矩阵优化

(f[n] =1*f[n-1]+1*f[n-2])
(f[n-1]=1*f[n-1]+0*f[n-2])
这样每乘上一个(inom{1,1}{1,0})，斐波拉锲数列就会向前推一个数（由(f[n-1]变成了f[n])）。
同样，将(inom{f[n-1]}{f[n-2]})展开
(inom{f[n-1]}{f[n-2]})=(inom{1，1}{1，0}*inom{f[n-2]}{f[n-3]})
这样不断展开，最终
(inom{f[n]}{f[n-1]})=(inom{1，1}{1，0}*inom{1，1}{1，0}*inom{1，1}{1，0}*inom{1，1}{1，0}*inom{1，1}{1，0}* ... *inom{1，1}{1，0}*inom{f[2]}{f[1]})

又因为矩阵乘法具有结合律
(inom{f[n]}{f[n-1]}=inom{1,1}{1,0}^{n-2}*inom{f[2]}{f[1]})，
即(inom{f[n]}{f[n-1]}=inom{1,1}{1,0}^{n-2}*inom{1}{1})，
恩~，这是快速幂的味道。。（逃）
以下是代码，
一般采用结构体存储矩阵

const int N=4;

struct matrix
{
	int n,m; //行数，列数
	LL g[N][N]; //矩阵
	void set(int x) //建立x*x的单位矩阵 
	{
		n=m=x;
		for(int i=1;i<=x;i++) 
			for(int j=1;j<=x;j++) 
				g[i][j]=(i==j);
		return;
	}
	void clear() //矩阵初始化 ，全部设成0
	{
		n=m=0;
		memset(g,0,sizeof(g));
		return;	
	} 
};

单位矩阵是指一个除对角线为1以外，其他数为0的矩阵。
如(inom{1,0}{0,1}),其性质是无论什么矩阵与之相乘，都得到原矩阵。
相当于整数中的1.

矩阵乘法

const LL MOD=1e9+7;

matrix operator*(matrix x,matrix y) //重载运算符
{
	matrix z;
	z.clear();
	z.n=x.n; z.m=y.m;
	for(int i=1;i<=z.n;i++) {
		for(int j=1;j<=z.m;j++) {
			LL ans=0;
			for(int k=1;k<=x.m;k++) 
				ans+=(LL)x.g[i][k]*y.g[k][j]%MOD;
			z.g[i][j]=ans%MOD;
		}
	}
	return z;
}

矩阵快速幂

就跟普通的快速幂一样啦。
不知道快速幂的同学移步这里。

matrix operator^(matrix x,LL k) //需要矩阵乘法
{
	matrix res;
	res.set(3); //注意根据题目修改
	while(k) {
		if(k&1) res=res*x;
		x=x*x; k>>=1;
	}
	return res;
}

输出矩阵

。。。再水一段代码。。

void matrix_put(matrix a)
{
	for(int i=1;i<=a.n;i++) {
		for(int j=1;j<=a.m;j++) {
			printf("%lld ",a.g[i][j]%MOD);
		}
		puts("");
	}
	return;
}

例题乱打

洛谷 P1962 斐波那契数列
题目背景
大家都知道，斐波那契数列是满足如下性质的一个数列：
(f[n]=1; (n<=2))
(f[n]=f[n-1]+f[n-2] (n>=3))
题目描述
请你求出 (F~n~ mod 10^9 + 7)的值。
输入格式
一行一个正整数 n
输出格式
输出一行一个整数表示答案。
输入 #1

5

输出 #1

5

输入 #2

10

输出 #2

55

说明/提示
【数据范围】
对于 60% 的数据，(1≤n≤92)；
对于 100%的数据，(1le n < 2^{63});
思路
这是道裸的矩阵乘法。
构造矩阵(inom{1,1}{1,0})
那么 (inom{f[n]}{f[n-1]}=inom{1,1}{1,0}^{n-2}*inom{1}{1})
剩下的交给模板解决。

#include<bits/stdc++.h>
#define rint register int
#define INF 2147483646l
#define LL long long
using namespace std;
const LL MOD=1e9+7;
LL n;
const int N=10;
LL ago[6][6]= //即斐波拉锲数列的矩阵
{
	{0},
	{0,1,1},
	{0,1,0}
};

struct matrix
{
	int n,m;
	LL g[N][N];
	void set(int x) 
	{
		n=m=x;
		for(int i=1;i<=x;i++) 
			for(int j=1;j<=x;j++) 
				g[i][j]=(i==j);
		return;
	}
};
matrix ss;
matrix operator*(matrix x,matrix y) //矩阵乘法
{
	matrix z;
	memset(z.g,0,sizeof(z.g));
	z.n=x.n; z.m=y.m;
	for(int i=1;i<=z.n;i++) {
		for(int j=1;j<=z.m;j++) {
			LL ans=0;
			for(int k=1;k<=x.m;k++) 
				ans+=(x.g[i][k]*y.g[k][j])%MOD;
			z.g[i][j]=ans%MOD;
		}
	}
	return z;
}
matrix operator^(matrix x,LL k) //矩阵快速幂
{
	matrix res;
	memset(res.g,0,sizeof(res.g));
	res.set(2);
	while(k)
	{
		if(k&1) res=res*x;
		x=x*x; k>>=1;
	}
	return res;
}
matrix c;
int main()
{
//	freopen("1.in","r",stdin);
	scanf("%lld",&n);
	for(int i=1;i<=2;i++) 
		for(int j=1;j<=2;j++)
			ss.g[i][j]=ago[i][j];
	ss.n=ss.m=2;
	c.n=2; c.m=1;
	c.g[1][1]=1; //原始矩阵
	c.g[2][1]=1;
	if(n==1||n==2) {
		puts("1"); return 0;
	}
	if(n==3) {
		puts("2"); return 0;
	}
	c=(ss^(n-2))*c;
	printf("%lld
",c.g[1][1]%MOD);
	return 0;
}