矩阵乘法
矩阵乘法(英文: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]\)
这样每乘上一个\(\binom{1,1}{1,0}\),斐波拉锲数列就会向前推一个数(由\(f[n-1]变成了f[n]\))。
同样,将\(\binom{f[n-1]}{f[n-2]}\)展开
\(\binom{f[n-1]}{f[n-2]}\)=\(\binom{1,1}{1,0}*\binom{f[n-2]}{f[n-3]}\)
这样不断展开,最终
\(\binom{f[n]}{f[n-1]}\)=\(\binom{1,1}{1,0}*\binom{1,1}{1,0}*\binom{1,1}{1,0}*\binom{1,1}{1,0}*\binom{1,1}{1,0}* ... *\binom{1,1}{1,0}*\binom{f[2]}{f[1]}\)
又因为矩阵乘法具有结合律
\(\binom{f[n]}{f[n-1]}=\binom{1,1}{1,0}^{n-2}*\binom{f[2]}{f[1]}\),
即\(\binom{f[n]}{f[n-1]}=\binom{1,1}{1,0}^{n-2}*\binom{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的矩阵。
如\(\binom{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~ \bmod 10^9 + 7\)的值。
输入格式
一行一个正整数 n
输出格式
输出一行一个整数表示答案。
输入 #1
5
输出 #1
5
输入 #2
10
输出 #2
55
说明/提示
【数据范围】
对于 60% 的数据,\(1≤n≤92\);
对于 100%的数据,\(1\le n < 2^{63}\);
思路
这是道裸的矩阵乘法。
构造矩阵\(\binom{1,1}{1,0}\)
那么 \(\binom{f[n]}{f[n-1]}=\binom{1,1}{1,0}^{n-2}*\binom{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\n",c.g[1][1]%MOD);
return 0;
}