矩阵加速递推浅谈


前言:

以下是一位矩阵蒟蒻学习矩阵加速的悲哀历史

(我才不会说是我不认真学)


前置芝士——矩阵快速幂:

矩阵加速,顾名思义,用到矩阵进行加速,而这个算法的加速核心就是矩阵快速幂 。(我只会背板子,逃~

在此对矩阵快速幂略作解释 :

快速幂大家应该都知道,矩阵快速幂就是把普通乘法换为矩阵乘法。

而矩阵乘法大家也应都有了解,大致方式如下:

[egin{bmatrix}a&c\b&dend{bmatrix} imes egin{bmatrix}e&g\f&hend{bmatrix}= egin{bmatrix}ae!+!cf&ag!+!ch\be!+!df&bg!+!dhend{bmatrix} quad ]

因此矩阵平方:

[egin{bmatrix}a_1&a_3\a_2&a_4end{bmatrix}^2= egin{bmatrix}a_1!a_1!+!a_3!a_2&a_1!a_3!+!a_3!a_4 \a_2!a_1!+!a_4!a_2&a_2!a_3!+!a_4!a_4end{bmatrix} ]

各位应该也能发现一个矩阵乘法的性质:如果两个矩阵能相乘,则第一个矩阵的行数与第二个矩阵的列数必须相等。

所以我们接下来构造矩阵,一般构造成 (n! imes!n)的形式,以便计算。


切入正题——矩阵加速:

DP or 递推是我们在经常用到的做题方法,一般时间复杂度为(O(n))(O(n^2)),但遇到一些比较恶心的出题人,给你来一个(nle10^{16}),你就可以和这个可爱的世界说再见了。而这个时候,便可能要用到矩阵加速。

例题1 :(P1962)斐波那契数列

这道题我们已经很熟了,但大家也应该清楚,随着数据量的增加,橙->绿->紫,紫题今天尚且不讲 (其实是我不会),我们来谈谈(n<2^{63})

本蒟蒻语文不好,没法子引入,所以也不拐弯抹角了,直接上矩阵。

众所周知,(f[n]=f[n-1]+f[n-2]),所以我们可以看成$ f[n+2]=f[n+1]* 1+f[n]* 1 $ 因此矩阵就构成了:

[egin{bmatrix}f[n]&f[n+1]end{bmatrix} imes egin{bmatrix}?&1\?&1end{bmatrix}= egin{bmatrix}?&f[n+2]end{bmatrix} ]

为了保证矩阵不断相乘的连续性,即因为第一个矩阵左边(f)的下标比右边小一,所以得出矩阵也应如此,所以:

[egin{bmatrix}f[n]&f[n+1]end{bmatrix} imes egin{bmatrix}?&1\?&1end{bmatrix}= egin{bmatrix}f[n+1]&f[n+2]end{bmatrix} ]

接下来中间矩阵大家就应该知道长啥样了,所以整个柿子为:

[egin{bmatrix}f[n]&f[n+1]end{bmatrix} imes egin{bmatrix}0&1\1&1end{bmatrix}= egin{bmatrix}f[n+1]&f[n+2]end{bmatrix} ]

又因为:

[egin{bmatrix}f[n+1]&f[n+2]end{bmatrix} imes egin{bmatrix}0&1\1&1end{bmatrix}= egin{bmatrix}f[n+2]&f[n+3]end{bmatrix} ]

所以以此类推:

[egin{bmatrix}f[1]&f[2]end{bmatrix} imes egin{bmatrix}0&1\1&1end{bmatrix}^{n-1}= egin{bmatrix}f[n]&f[n+1]end{bmatrix} ]

柿子推完了,接下来就是矩阵快速幂的事了,代码网上满天飞,我就不贴了。(就是我懒)

另外,

这题的柿子还可以写成如下形式:

[egin{bmatrix}f[1]\f[2]end{bmatrix} imes egin{bmatrix}0&1\1&1end{bmatrix}^{n-1}= egin{bmatrix}f[n]\f[n+1]end{bmatrix} ]

虽然说相乘方式变了,但为什么就不能说第一个矩阵是第二个,第二个矩阵是第一个?大家可以选择喜欢的方式推柿子。

(本人推荐第二种,原因在当你做了一道duliu矩阵加速题时便会知道)


例题2 :(P1707)刷题比赛

本题比较复杂 (恶心出题人),一大堆常数一下让本蒟蒻不知所措,但按照套路列两个矩阵,在慢慢算中间矩阵,还是可以的,就是比较烦。。。

但如此傲娇的我,怎能在这种事情上停住脚步?

矩阵如下:

[egin{bmatrix}a[n]\a[n+1]\b[n]\b[n+1]\c[n]\c[n+1]\k^2\k\1\w^{k}\z^{k}end{bmatrix} imes egin{bmatrix}0&1&0&0&0&0&0&0&0&0&0\q&p&0&1&0&1&r&t&1&0&0\0&0&0&1&0&0&0&0&0&0&0\0&1&v&u&0&1&0&0&0&1&0\0&0&0&0&0&1&0&0&0&0&0\0&1&0&1&y&x&0&1&2&0&1\0&0&0&0&0&0&1&2&1&0&0\0&0&0&0&0&0&0&1&1&0&0\0&0&0&0&0&0&0&0&1&0&0\0&0&0&0&0&0&0&0&0&w&0\0&0&0&0&0&0&0&0&0&0&zend{bmatrix}= egin{bmatrix}a[n+1]\a[n+2]\b[n+1]\b[n+2]\c[n+1]\c[n+2]\k^2+2k+1\k+1\1\w^{k+1}\z^{k+1}end{bmatrix} ]

终于打完了,累呀~~~

但我相信大家以后遇到这种duliu题应该都能有思路了。

终于到贴代码的时候了:

#include<bits/stdc++.h>
#define re register
using namespace std;
long long p,q,r,t,u,v,w,x,y,z;
long long n,mod,f[12],mp[12][12];
inline long long Mul(long long a,long long b)
{
	re long long ans=0;
	for(;b;b>>=1)
	{
		if(b&1)ans=(ans+a)%mod;
		a=(a<<1)%mod;
	}
	return ans;
}
inline void work()
{
	re long long mid[12];
	memset(mid,0,sizeof(mid));
	for(re int i=1;i<=11;i++)
		for(re int j=1;j<=11;j++)
			mid[i]=(mid[i]+Mul(f[j],mp[i][j]))%mod;
	for(re int i=1;i<=11;i++)f[i]=mid[i];
}
inline void Pow()
{
	re long long mid[12][12];
	memset(mid,0,sizeof(mid));
	for(re int i=1;i<=11;i++)
		for(re int j=1;j<=11;j++)
			for(re int k=1;k<=11;k++)
				mid[i][j]=(mid[i][j]+Mul(mp[i][k],mp[k][j]))%mod;
	for(re int i=1;i<=11;i++)
		for(re int j=1;j<=11;j++)
			mp[i][j]=mid[i][j];
}
int main()
{
	scanf("%lld%lld%lld%lld%lld%lld%lld%lld%lld%lld%lld%lld",&n,&mod,&p,&q,&r,&t,&u,&v,&w,&x,&y,&z);
	f[1]=f[3]=f[5]=f[7]=f[8]=f[9]=1;
	f[2]=f[4]=f[6]=3,f[10]=w,f[11]=z;
	mp[1][2]=mp[2][4]=mp[2][6]=mp[2][9]=mp[3][4]=mp[4][2]=mp[4][6]=mp[4][10]=mp[5][6]=mp[6][2]=mp[6][4]=mp[6][8]=mp[6][11]=mp[7][7]=mp[7][9]=mp[8][8]=mp[8][9]=mp[9][9]=1;
	mp[2][1]=q,mp[2][2]=p,mp[2][7]=r,mp[2][8]=t,mp[4][3]=v,mp[4][4]=u,mp[6][5]=y,mp[6][6]=x,mp[10][10]=w,mp[11][11]=z,mp[6][9]=mp[7][8]=2;
	re long long res=n-1ll;
	for(;res;Pow(),res>>=1)
		if(res&1)work();
	printf("nodgd %lld
Ciocio %lld
Nicole %lld
",f[1],f[3],f[5]);
	return 0;
}
/*
nodgd 74
Ciocio 80
Nicole 59
*/

提醒:本题需要用到龟速乘(会爆long long)~

代码如下:

const long long mod=998244353;
inline long long Mul(long long a,long long b)
{
	re long long ans=0;
	for(;b;b>>=1)
	{
		if(b&1)ans=(ans+a)%mod;
		a=(a<<1)%mod;
	}
	return ans;
}

当然你如果高兴也可以打高精。


都看到这里了,给个赞吧!

原文地址:https://www.cnblogs.com/jkzcr01-QAQ/p/13575194.html