矩阵乘法
矩阵乘法顾名思义,就是将两个矩阵做乘法运算(相当于在矩阵意义下重载乘法运算符?)
矩阵乘法的定义
矩阵A $$\times $$ 矩阵B = 矩阵C
\[C_{i,j}=\sum_{k=1}^{m}A_{i,k}\times B_{k,j} \]矩阵C的第i行第j列元素是矩阵A第i行的所有元素与矩阵B的第j列的所有元素对应相乘的乘积的总和。
-
对应就要求矩阵A的列数要等于矩阵B的行数。
-
即A是\(M\times N\)的矩阵,B是\(N\times K\)的矩阵
同时矩阵乘法,满足结合律,分配律,但不满足交换律(一是最终形成的矩阵的行数和列数不一定一样,即便是一样,矩阵里的数据也不一定一样)。
矩阵乘法的特殊情形
矩阵A是一个\(N\times N\)的矩阵,矩阵F是一个\(N\times 1\)的矩阵,设\(F_{1}=A\times F\),发现\(F_{1}\)也是一个\(N\times 1\)的矩阵,是一个有只有一列元素的元素,我们不妨把这些元素看成是一个个变量,而矩阵\(F_{1}\) 的变量是通过矩阵\(F\)的变量的某种关系来实现更新,而这种关系是记录到矩阵A里面的(所以矩阵A也可以称作是系数矩阵?)
- 有将一维度的矩阵F定义为状态矩阵,并将不变且用来相乘的矩阵A称作转移矩阵。
比如斐波那契数列,\(F_{n+2}=F_{n+1}+F_{n}\),在不考虑对所有的\(F_{i}\)的值就进行记录时,即\(\begin{cases} F_{n+2}=F_{n+1}+F_{n}\\F_{n+1}=F_{n+1}\end{cases}\)在空间上进行优化可优化为,\(\begin{cases} F_{n+1}=F_{n+1}+F_{n}\\F_n=F_{n+1}\end{cases}\),转化为矩阵算法有\(\begin{bmatrix} F_{n+1}&F_{n} \end{bmatrix}=\begin{bmatrix} F_{n+1}&F_{n} \end{bmatrix}\times \begin{bmatrix} 1&1\\1&0 \end{bmatrix}\)
应用
上述的矩阵乘法的特殊情形常用于递推式,特别是当递推的层次达到一定的数量级,就可以用矩阵快速幂来减少计算次数,优化计算速度,从而提高效率,而一般的递推只能重复性的计算(和覆盖)。
实现
-
矩阵C的第i行第j列元素是矩阵A第i行的所有元素与矩阵B的第j列的所有元素对应相乘的乘积的总和。
for(int i=0;i<N;i++) for(int j=0;j<K;j++)//合成完后是一个N*K的矩阵 for(int k=0;k<M;k++) C[i][j]+=A[i][k]*B[k][j];
-
而如果要从\(F_{1}\)出发来得到\(F_{N}\),能不能有一种方法来减少计算量呢?
-
发现\(F_{N}=A\times F_{N-1}=A^{N}\times F_{0}\)(由递推得到)
-
A的n次方是否要一个一个A来连乘来获得?是否有什么方法来快速得到答案?
-
快速幂?!矩阵快速幂?!于是这块运算的级别就可以从\(O(N)\)级别降低到 \(O(logN)\)级别,最终再乘上内部单轮矩阵乘法共需要的时间(三重循环)来得到最终的复杂度。
-
代码框架
定义一个存放矩阵的结构体
struct Mat{
ll mat[110][110];
Mat() {}
//结构体内重载运算符
Mat operator*(Mat const &b)const
{
Mat res;
memset(res.mat,0,sizeof(res.mat));
for(int i=0;i<maxn;i++)
for(int j=0;j<maxn;j++)
for(int k=0;l<maxn;k++)
res.mat[i][j]=(res.mat[i][j]+this->mat[i][k]*b.mat[k][j]);
return res;
}
}
定义一个返回类型为mat的,作用为实现矩阵快速幂的函数
Mat pow_mod(Mat base,ll n)
{
Mat res;
memset(res,0,sizeof(res));
for(int i=0;i<maxn;i++) //定义一个单位矩阵
res.mat[i][i] = 1;
while(n>0)
{
if(n&1)res=res*base;
base=base*base;
n>>=1;
}
return res;
}
AcWing 205. 斐波那契
#include <bits/stdc++.h>
#define MEM(a,x) memset(a,x,sizeof(a))
#define W(a) while(a)
#define gcd(a,b) __gcd(a,b)
#define pi acos(-1.0)
#define PII pair<int,int>
#define pb push_back
#define mp make_pair
#define fi first
#define se second
#define ll long long
#define ull unsigned long long
#define rep(i,x,n) for(int i=x;i<n;i++)
#define repd(i,x,n) for(int i=x;i<=n;i++)
#define MAX 1000005
#define MOD 10000
#define INF 0x3f3f3f3f
#define lowbit(x) (x&-x)
using namespace std;
int n,m;
void mul(int f[2],int A[2][2])//f[1][2]
{
int C[2];//C[1][2]
memset(C,0,sizeof(C));
for(int j=0;j<2;j++)
for(int k=0;k<2;k++)
C[j]=(C[j]+(ll)f[k]*A[k][j])%MOD;//C[i]相当于是C的第一行的第i个元素
memcpy(f,C,sizeof(C));
}
void mulself(int A[2][2])
{
int C[2][2];
memset(C,0,sizeof(C));
for(int i=0;i<2;i++)
for(int j=0;j<2;j++)
for(int k=0;k<2;k++)
C[i][j]=(C[i][j]+(ll)A[i][k]*A[k][j])%MOD;
memcpy(A,C,sizeof(C));
}
int main()
{
while(cin>>n && n!=-1)
{
int f[2] = {0,1};
int A[2][2]={
{0,1},
{1,1}
};
while(n>0)
{
if(n&1)mul(f,A);
mulself(A);
n=n>>1;
}
cout<<f[0]<<endl;
}
return 0;
}
题意:
Reziba拥有许多神奇的宝石。每一颗魔法石都可以切分成M块普通石头。从占用的角度上来说,无论是神奇的宝石还是普通的宝石都会占据一个单位的空间,并且普通的宝石无法继续切割。
Rezilba想要选择一套魔法石并且切割其中的一部分,以至于最终这些宝石能够占用到N个单位空间。如果神奇的宝石有被切割成M份,那就当做是占据了M个空间,否则还是当做只占用了一个空间。
给定一个数字N和M,N是最终构成的宝石数目,M是普通宝石和神奇宝石的比值。请问在这种转换比的条件下,有多少种方案能到达N颗。
补题思路:
要达到刚刚好占用n个空间有两种办法,一种是在已经占用了n-1个空间的时候再添置一个普通宝石,另一种的方法是在已经占用n-m个空间的时候,释放一颗魔法石,来占用m个空间。
f[n]=f[n-1]+f[n-m];
仅当n>=m时,否则的话f[n]=f[n-1];
#include <bits/stdc++.h>
#define MEM(a,x) memset(a,x,sizeof(a))
#define W(a) while(a)
#define gcd(a,b) __gcd(a,b)
#define pi acos(-1.0)
#define PII pair<int,int>
#define pb push_back
#define mp make_pair
#define fi first
#define se second
#define ll long long
#define ull unsigned long long
#define rep(i,x,n) for(int i=x;i<n;i++)
#define repd(i,x,n) for(int i=x;i<=n;i++)
#define MAX 1000005
#define MOD 1000000007
#define INF 0x3f3f3f3f
#define lowbit(x) (x&-x)
using namespace std;
ll n;int m;
struct Mat{
ll mat[110][110];
Mat(){};
Mat operator*(Mat const &b)const
{
Mat res;
memset(res.mat,0,sizeof(res.mat));
rep(k,0,m)
rep(i,0,m)
rep(j,0,m)
res.mat[i][j]=(res.mat[i][j]+(this->mat[i][k]*b.mat[k][j]))%MOD;
return res;
}
};
Mat mat_pow(Mat A,ll k)
{
Mat res;
memset(res.mat,0,sizeof(res.mat));
rep(i,0,m)
res.mat[i][i]=1;
while(k>0)
{
if(k&1)res=res*A;
A=A*A;
k=k>>1;
}
return res;
}
int main()
{
Mat base,A;
scanf("%lld%d",&n,&m);//注意n是lld
if(n<m)
{
printf("1\n");
return 0;
}
memset(base.mat,0,sizeof(base.mat));
memset(A.mat,0,sizeof(A.mat));
A.mat[0][0]=1,A.mat[0][m-1]=1;
rep(i,1,m)
A.mat[i][i-1]=1;
rep(i,0,m)
base.mat[i][0]=1;
A=mat_pow(A,n-m+1);
base=A*base;
printf("%lld",base.mat[0][0]);//在乘法的过程中已经做过了取余
return 0;
}
其他
单位矩阵
一个矩阵乘上单位矩阵后保持不变(前提是能够相乘)
memcpy函数
- C 库函数 void *memcpy(void *str1, const void *str2, size_t n) 从存储区 str2 复制 n 个字节到存储区 str1。
矩阵乘法的优化
采取k,i,j的三重循环顺序相比于i,j,k能提高cache的命中率,从而提高运行效率。