zoukankan html css js c++ java

矩阵加速学习笔记

基本知识

首先关于矩阵的一些基本常识在这篇博客里已经介绍的很清楚了，在这里不在赘述。

但是矩阵乘的图片是一定放的，因为比较重要好像只有矩阵乘，~~其实看了图片好像根本不用看博文~~。

观察图片我们可以发现，对于一个新矩阵的 (F[i][j]) 其实就是第一个矩阵的第 (i) 行和对应第二个矩阵的 (j) 列的每一个数乘起来再相加。这也就解释了为什么只有大小为 (n*m) 和 (m*k) 的矩阵才可以相乘，且新矩阵为大小为 (n*k) ，如果不是这样，就会出现对不上项的情况。

和实数一样，在定义了矩阵乘之后矩阵幂也就自然出来了，只有 (n*n) 的矩阵才可以不断连乘。

就像上面的新矩阵 (n*k) ，如果 (n) (m) (k) 不相等，那么就会发现没有旧矩阵和新矩阵可以不断相乘。

还有就是单位矩阵，任何矩阵乘这个矩阵都是其本身，就像实数里的 (1) 一样，定义上面的博文里已经很详细了，这里只是口糊一下，就是只有对角线元素是 (1) 其余元素都为 (0) 的矩阵。

典型例题

P3390 【模板】矩阵快速幂

模板题，直接放代码

个人写的一个矩阵的模板。

#include<bits/stdc++.h>
#define ll long long
const int N=110;
const ll MOD=1e9+7;
ll n,k;
struct Matrix{
	int n;
	ll a[N][N];
	inline void clear(int size)
	{
		n=size;
		memset(a,0,sizeof(a));
		return ;
	}
	inline void build(int size)
	{
		n=size;
		memset(a,0,sizeof(a));
		for (int i=1;i<=n;i++)
		a[i][i]=1;
		return ;
	}
	inline void out()
	{
		for (int i=1;i<=n;i++)
		{
			for (int j=1;j<=n;j++)
			printf("%lld ",a[i][j]);
			printf("
");
		}
	}
};
inline Matrix operator * (Matrix a,Matrix b)
{
	Matrix c;
	c.clear(a.n);
	for (int i=1;i<=c.n;i++)
	for (int j=1;j<=c.n;j++)
	for (int k=1;k<=c.n;k++)
	c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j]%MOD)%MOD;
	return c;
}
inline Matrix Matrix_qpow(Matrix a,ll p)
{
	Matrix ans,base=a;
	ans.build(a.n);
	for (;p;p>>=1,base=base*base)
	if (p&1) ans=ans*base;
	return ans;
}
int main()
{
	scanf("%lld%lld",&n,&k);
	Matrix now;
	now.clear(n);
	for (int i=1;i<=n;i++)
	for (int j=1;j<=n;j++)
	scanf("%lld",&now.a[i][j]);
	Matrix_qpow(now,k).out();
	return 0;
}

需要注意的一个问题是如果存矩阵的结构体占用空间过大，当其传入自定义函数的时候会引起爆栈。

可以在编译命令里加入无限栈空间的命令。

-Wl,--stack=1234567890

那么说了那么多，矩阵乘和矩阵快速幂有什么用呢？？

矩阵快速幂可以优化一些递推式，使 (Theta(n)) 的时间复杂度变为 (Theta (log_2n)) 。

下面的题目就是一个比较简单的应用。

P1962 斐波那契数列

思路主要来自这篇博客，这里主要就是口述一下过程。

首先我们知道斐波那契数列 (F[n]=F[n-1]+F[n-2]) ，对于 (F[n]) 这和 (F[n-1]) 和 (F[n-2]) 有关，我们设 (Fib[n-1]) 为一个 (F[n-1]) 和 (F[n-2]) 的二元组，所以 (Fib[n]) 就是 (F[n]) 和 (F[n-1]) ，考虑 (Fib[n]) 怎样由 (Fib[n-1]) 得到。我们发现

[F[n]=F[n-1]*1+F[n-2]*1 ]

[F[n-1]=F[n-1]*1+F[n-2]*0 ]

如果放在矩阵上，它就是这样的。

我们可以发现，如果用这种方法构造矩阵，第 (i) 行矩阵的第 (j) 列的数就是 (Fib[n]) 中的第 (i) 项和 (Fib[n-1]) 中第 (j) 项的系数关系。

然后我们发现求第 (n) 项的过程其实就是不断进行上面的过程的过程，最终的答案就是求下面这个矩阵的第一行第一个元素

[left[egin{array}{cc}{F_{n-1}} & {F_{n-2}}end{array} ight] imesleft[egin{array}{cc}{1} & {1} \ {1} & {0}end{array} ight]^{n-2} ]

因为矩阵满足结合律，所以我们就先用矩阵快速幂求出来(left[egin{array}{ll}{1} & {1} \ {1} & {0}end{array} ight]^{n-2}) ，然后再求 (left.[egin{array}{ll}{1} & {1}end{array} ight] imesleft[egin{array}{ll}{1} & {1} \ {1} & {0}end{array} ight]^{n-2}) ，需要注意的是矩阵并不满足交换律，所以不能换过来乘。

关于单位矩阵的构造的问题。

其实就乘起来让目标矩阵和原矩阵一样的矩阵，可以先按原矩阵的第一行推矩阵，然后就会推出一个单位矩阵，然后你会发现这个矩阵适用于原矩阵每一行的情况。

#include<bits/stdc++.h>
#define ll long long
struct Mat{
	int size;
	ll **M=NULL;
	inline ll Start()
	{
		if (M!=NULL) return M[1][1];
		else return LLONG_MIN;
	}
	inline void Clear(int sz)
	{
		if (M!=NULL)
		for (int i=0;i<=sz;i++)
		for (int j=0;j<=sz;j++)
		M[i][j]=0;
		return ;
	}
	inline void New(int sz)
	{
		size=sz;
		M=new ll*[sz+10];
		for (int i=0;i<sz+10;i++)
		M[i]=new ll[sz+10];
		Clear(sz);
		return ;
	}
	inline void Build(int sz)
	{
		size=sz;
		if (M==NULL) New(sz);
		for (int i=1;i<=sz;i++) M[i][i]=1;
		return ;
	}
	inline void Init(ll now[],int sz)
	{
		if (M==NULL) New(sz);
		int num=0;
		for (int i=1;i<=sz;i++)
		for (int j=1;j<=sz;j++)
		M[i][j]=now[++num];
		return ;
	}
	inline void Out()
	{
		if (M!=NULL)
		for (int i=1;i<=size;i++)
		{
			for (int j=1;j<=size;j++)
			printf("%lld ",M[i][j]);
			printf("
");
		}
		return ;
	}
	inline void Delete()
	{
		for (int i=0;i<size+10;i++)
		delete []M[i] ;
		delete []M;
		M=NULL;
		size=0;
		return ;
	}
};
ll mod;
inline Mat operator * (Mat a,Mat b)
{
	Mat c;
	c.New(a.size);
	c.Clear(a.size);
	for (int i=1;i<=c.size;i++)
	for (int j=1;j<=c.size;j++)
	for (int k=1;k<=c.size;k++)
	c.M[i][j]=(c.M[i][j]+a.M[i][k]*b.M[k][j]%mod)%mod;
	return c;
}
inline Mat Mat_qpow(Mat a,ll p)
{
	Mat ans,base;
	base.New(a.size);
	ans.Build(a.size);
	for (base=a;p;p>>=1,base=base*base)
	if (p&1) ans=ans*base;
	return ans;
}
ll n;
int main()
{
	scanf("%lld",&n);
	if (n==1||n==2)
	{
		printf("1
");
		return 0;
	}
	mod=(ll)1e9+7;
	ll a[]={0,1,1,1,0};
	ll b[]={0,1,1,0,0};
	Mat ans,use;
	use.Init(a,2);
	ans.New(2);
	ans=Mat_qpow(use,n-2);
	use.Init(b,2);
	ans=use*ans;
	printf("%lld
",ans.Start());
	ans.Delete();use.Delete();
	return 0;
}

P1939 【模板】矩阵加速（数列）

和上面的那一道题目推法差不多，一般来说要求的项数前与 (n) 项相关，矩阵就是 (n*n) 的。

#include<bits/stdc++.h>
#define ll long long
const int N=12;
const ll MOD=1e9+7;
struct Matrix{
	int n;
	ll a[N][N];
	inline void clear(int sz)
	{
		n=sz;
		memset(a,0,sizeof(a));
		return ;
	}
	inline void build(int sz)
	{
		n=sz;
		memset(a,0,sizeof(a));
		for (int i=1;i<=n;i++) a[i][i]=1;
		return ;
	}
	inline void out()
	{
		for (int i=1;i<=n;i++)
		{
			for (int j=1;j<=n;j++)
			printf("%lld ",a[i][j]);
			printf("
");
		}
		return ; 
	}
};
inline Matrix operator * (Matrix a,Matrix b)
{
	Matrix c;
	c.clear(a.n);
	for (int i=1;i<=c.n;i++)
	for (int j=1;j<=c.n;j++)
	for (int k=1;k<=c.n;k++)
	c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j]%MOD)%MOD;
	return c;
}
inline Matrix Matrix_qpow(Matrix a,ll p)
{
	Matrix ans,base=a;
	ans.build(a.n);
	for (;p;p>>=1,base=base*base)
	if (p&1) ans=ans*base;
	return ans;
}
int t;
int main()
{
	scanf("%d",&t);
	while (t--)
	{
		ll k;
		scanf("%lld",&k);
		if (k<=3)
		{
			printf("1
");
			continue;
		}
		Matrix now,ans;
		now.clear(3);
		now.a[1][1]=1;now.a[1][2]=1;now.a[1][3]=0;
		now.a[2][1]=0;now.a[2][2]=0;now.a[2][3]=1;
		now.a[3][1]=1;now.a[3][2]=0;now.a[3][3]=0;
		ans=Matrix_qpow(now,k-3);
		printf("%lld
",(ans.a[1][1]+ans.a[2][1]+ans.a[3][1])%MOD);
	}
	return 0;
}

不过在最后乘的时候一定要严格按照矩阵乘法的定义来，因为你永远不知道经过矩阵快速幂的矩阵是什么。

P1349 广义斐波那契数列

只是加了系数，和上面的推法一样，需要注意的是定义的递推矩阵和给的顺序正好相反。

#include<bits/stdc++.h>
#define ll long long
const int N=5;
struct Matrix{
	int n;
	ll a[N][N];
	inline void clear(int sz)
	{
		n=sz;
		memset(a,0,sizeof(a));
		return ;
	}
	inline void build(int sz)
	{
		n=sz;
		memset(a,0,sizeof(a));
		for (int i=1;i<=n;i++)
		a[i][i]=1;
		return ;
	}
};
ll q,p,a1,a2,n,mod;
inline Matrix operator * (Matrix a,Matrix b)
{
	Matrix c;
	c.clear(a.n);
	for (int i=1;i<=c.n;i++)
	for (int j=1;j<=c.n;j++)
	for (int k=1;k<=c.n;k++)
	c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j]%mod)%mod;
	return c;
}
inline Matrix Matrix_qpow(Matrix a,ll p)
{
	Matrix ans,base=a;
	ans.build(a.n);
	for (;p;p>>=1,base=base*base)
	if (p&1) ans=ans*base;
	return ans;
}
int main()
{
	scanf("%lld%lld%lld%lld%lld%lld",&p,&q,&a1,&a2,&n,&mod);
	Matrix now,ans;
	now.clear(2);
	now.a[1][1]=p;now.a[1][2]=1;
	now.a[2][1]=q;now.a[2][2]=0;
	ans=Matrix_qpow(now,n-2);
	printf("%lld
",(ans.a[1][1]*a2%mod+ans.a[2][1]*a1%mod)%mod);
	return 0;
}

P5343 【XR-1】分块

一看 (60) 分的暴力，开心 (DP) 。

然后我们发现对于每一个 (F[n]) 好像只与 (F[n-1]) ， (F[n-2]) …… 最多只与 (100) 个前面的状态有关，所以我们考虑矩阵加速。

我们设初始矩阵全部为 (1) ，就是当 (n) 等于 (0) 的方案数，转移矩阵也比较好推。

！！！矩阵不要乘反，矩阵不要乘反，矩阵不要乘反，重要的事情说三遍。

一定要看清数据范围，一定要分清矩阵的第 (i) 项代表什么。

#include<bits/stdc++.h>
#define ll long long
const ll MOD=1e9+7;
const int N=110;
struct Mat{
	int n;
	ll a[N][N];
	inline ll answer() {return a[1][1];}
	inline void clear(int sz)
	{
		n=sz;
		memset(a,0,sizeof(a));
		return ;
	}
	inline void build(int sz)
	{
		n=sz;
		memset(a,0,sizeof(a));
		for (int i=1;i<=n;i++)
		a[i][i]=1;
		return ;
	}
};
inline Mat operator * (Mat a,Mat b)
{
	Mat c;
	c.clear(a.n);
	for (int i=1;i<=c.n;i++)
	for (int j=1;j<=c.n;j++)
	for (int k=1;k<=c.n;k++)
	c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j]%MOD)%MOD;
	return c;
}
inline Mat Mat_qpow(Mat a,ll p)
{
	Mat ans,base=a;
	ans.build(a.n);
	for (;p;p>>=1,base=base*base)
	if (p&1) ans=ans*base;
	return ans;
}
ll l,PR,NF;
bool a[N],b[N];
int main()
{
	scanf("%lld",&l);
	scanf("%lld",&PR);
	for (int i=1;i<=PR;i++)
	{
		int x;
		scanf("%d",&x);
		a[x]=1;
	}
	scanf("%lld",&NF);
	for (int i=1;i<=NF;i++)
	{
		int x;
		scanf("%d",&x);
		b[x]=1;
	}
	for (int i=1;i<=100;i++)
	b[i]=b[i]&a[i];
	Mat now,ans;
	now.clear(100);
	for (int i=1;i<=100;i++)
	if (b[i]) now.a[i][1]=1;
	for (int i=2;i<=100;i++)
	now.a[i-1][i]=1;
	ans=Mat_qpow(now,l);
	now.clear(100);
	for (int i=1;i<=100;i++) now.a[i][1]=1;
	ans=now*ans;
	printf("%lld
",ans.answer());
	return 0;
}

需要注意的一些问题

矩阵不满足交换律，所以一定要分清是哪一个矩阵乘哪一个矩阵
一定要对应好初始矩阵中的每一项和递推项之间的关系
如果在最后想直接展开递推项乘矩阵得到答案一定要注意第二个问题和严格按照矩阵乘法的定义来计算

查看全文

相关阅读:
【NOIP2016提高组】 Day1 T3 换教室
 【NOIP2016 Day1 T2】天天爱跑步
 web@前端--html,css,javascript简介、第一个页面(常用标签简介)
MySQL--视图view、触发器trigger、事务(start transaction)、存储过程(特殊的数据逻辑处理函数)、流程控制(if，case....)
MySQL--pymysql模块
 MySQL--（了解）可能会用到的内置函数
 MySQL--详细查询操作(单表记录查询、多表记录查询(连表查询)、子查询)
MySQL--表操作(约束条件foreign key关联表多对1，多对多，1对1)
MySQL--表操作(innodb表字段数据类型、约束条件)、sql_mode操作
 MySQL基本语句、存储引擎

原文地址：https://www.cnblogs.com/last-diary/p/11724424.html

矩阵加速 学习笔记