【学习小记】常系数齐次线性递推

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/hzj1054689699/article/details/85683342

问题引入

给出数列gg,满足当n>mn>m
gn=i=1mgni×aig_n=\sum\limits_{i=1}^{m}g_{n-i}\times a_i
n<=mn<=m时,gn=cng_n=c_n

m比较小,n特别大,快速计算gng_n

Newbie的解法

暴力递推计算

时间复杂度O(nm)O(nm)

Pupil的解法

可以将转移和数列都写成m×mm\times m的矩阵的形式,矩阵快速幂即可

时间复杂度O(m3logn)O(m^3\log n)

Master的解法

我们需要一些数学知识进行铺垫:

Part 1 矩阵的特征值与特征多项式

我们知道一个矩阵乘一个列向量仍然是一个列向量。

若对于m阶矩阵A,有常数λ\lambda,非零列向量v\vec v,满足λv=Av\lambda\vec v=A\vec v则称λ\lambda为矩阵A的特征值,v\vec v为矩阵的特征向量

上式也可以写作(λIA)v=0(\lambda I-A)\vec v=0其中II为单位矩阵
此式有解的充要条件是λIA=0|\lambda I-A|=0,即矩阵λIA\lambda I-A的行列式为0

λIA|\lambda I-A|可以看做是关于λ\lambda的一个m次多项式,记作f(λ)f(\lambda)f(λ)f(\lambda)称作矩阵A的特征多项式,对于矩阵A的任意一个特征值λ0\lambda_0,都有f(λ0)=0f(\lambda_0)=0

Part 2 Hamilton-Cayley theorem

对于矩阵,也一样的定义多项式运算(把多项式中的x换乘矩阵A),加法就是直接对应相加,常数乘法就按位相乘,乘法是矩阵乘法,0次方是单位矩阵,它的结果仍然是一个矩阵。

显然,矩阵多项式满足交换律,即f(A)g(A)=g(A)f(A)f(A)g(A)=g(A)f(A)成立。
简单证明:考虑某两项相乘的结果Ax×AyA^x\times A^y,由于前后都是A,矩阵乘法满足结合律,因此指数可以任意分配,换成Ay×AxA^y \times A^x也是可以的

哈密顿—凯莱定理:对于矩阵A的特征多项式f(x)f(x),满足f(A)=0f(A)=0

证明网上到处都有,此处就不赘述了。

Part 3 求解转移矩阵的特征多项式

回到原题,我们对于Pupil解法的转移矩阵A,求解它的特征多项式
考虑矩阵λIA\lambda I-A

它长这样:
(1)λIA=(λa1a2am1am1λ000100001λ) \lambda I-A= \left( { \begin{matrix} \lambda-a_1 & -a_2 & \cdots &-a_{m-1} & -a_m \\ -1 & \lambda & \cdots & 0 &0 \\ 0 & -1 &\cdots & 0 & 0\\ \vdots & \vdots & \ddots & \vdots &\vdots \\ 0 & 0 & \cdots & -1 & \lambda \end{matrix} \tag{1} } \right)

根据行列式的定义,将第一行展开
λIA=(λa1)A1,1+a2×A1,2++am×A1,m|\lambda I-A|=(\lambda-a_1)A_{1,1}+a_2\times A_{1,2}+\cdots+a_m\times A_{1,m}
其中Ai,jA{i,j}表示矩阵A的代数余子式,即挖掉第i行和第j列以后剩下的矩阵的行列式。

我们发现所有的余子矩阵都是下三角矩阵,行列式就是对角线乘积。

化简整理,可得f(λ)=λIA=λmi=0m1amiλif(\lambda)=|\lambda I-A|=\lambda^m-\sum\limits_{i=0}^{m-1}a_{m-i}\lambda ^i
负号都被行列式里面逆序对个数的负号消掉了。

Part 4 计算答案

我们设要求的数列gg的初始矩阵为GG,它是一个m行1列的矩阵(列向量),从第m行到第1行分别为g1mg_{1\dots m}(注意顺序是反的)
实际上我们想知道的gng_n就是矩阵An1GA^{n-1}G的第m行第一列的值。

此时的关键就是An1A^{n-1},因为n1n-1非常大,无法直接计算

然而根据前面的铺垫,我们有f(A)=0f(A)=0An1A^{n-1}我们可以看做只有一项的一个关于A的多项式

那么根据多项式除法相关知识,可以得到An1=P(A)f(A)+Q(A)A^{n-1}=P(A)f(A)+Q(A),其中Q(A)Q(A)的次数是小于f(A)f(A)的次数也就是小于m的,Q(A)Q(A)相当于多项式An1A^{n-1}对多项式f(A)f(A)取模

可能会有这样的疑问,f(A)=0f(A)=0怎么能作除数呢?
其实并不要紧,我们并不需要知道f(A)f(A)的实际值,我们相当于将An1A^{n-1}减去了若干个f(A)f(A),将次数降低了,而结果不变。

实现上来说,由于ff的系数已知,我们可以先将式子里的矩阵A换成变量xx,代入,利用多项式取模算出Q的系数,然后再将x换回A,这样得出来的Q的系数是相同的。并且计算Q(A)×GQ(A)\times GAn1×GA^{n-1}\times G的结果是一样的。

为了求出Q(x)Q(x)的系数,我们可以采用快速幂的做法,初始Q0(x)=x1Q_0(x)=x^1,然后不断的自己与自己相乘,乘完对多项式f(x)f(x)取模
这一部分如果暴力取模,时间复杂度为O(m2logn)O(m^2\log n)
如果采用NTT优化多项式取模,时间复杂度为O(mlogmlogn)O(m\log m\log n)
这样求出了Q(A)Q(A)的系数,不妨设Q(A)=i=0m1diAiQ(A)=\sum\limits_{i=0}^{m-1}d_iA^i
要求矩阵Q(A)×GQ(A)\times G的第m行第一列的值

也就是i=0m1diAiG\sum\limits_{i=0}^{m-1}d_iA^iG的第m行第一列
然而AiGA^iG的第m行第一列的值就是gi+1g_{i+1}

所以gn=i=0m1digi+1=i=0m1dici+1g_n=\sum\limits_{i=0}^{m-1}d_ig_{i+1}=\sum\limits_{i=0}^{m-1}d_ic_{i+1}

还有一种情况,前m项并没有直接给出,也是通过递推得出的,暴力递推求前m项的复杂度是O(m2)O(m^2)
考虑优化

考虑数列gg的一般生成函数G(x)G(x)(与矩阵G不同)
转移序列aa的一般生成函数A(x)A(x)

由于G(x)G(x)是无限长的一个序列,我们可以得到G(x)=G(x)A(x)+rG(x)=G(x)A(x)+r
其中rr是一个常数,相当于第0项

移项,可以得到G(x)=r1A(x)G(x)={r\over 1-A(x)}
在模xm+1x^{m+1}意义下多项式求逆即可
时间复杂度是O(mlogm)O(m\log m)

模板题([BZOJ4161] Shlw loves matrixI)

Code

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <iostream>
#include <algorithm>
#define fo(i,a,b) for(int i=a;i<=b;++i)
#define fod(i,a,b) for(int i=a;i>=b;--i)
#define N 4005
#define mo 1000000007
#define LL long long
using namespace std;
LL f[N],g[N],h[N],s1[N],a[N],u1[N];
int n,m;
void mul(LL *x,LL *y,LL *z)
{
	fo(i,0,2*m-2) u1[i]=0;
	fo(i,0,m-1) fo(j,0,m-1) u1[i+j]=(u1[i+j]+x[i]*y[j])%mo;
	fod(i,2*m-2,m)
	{
		fo(j,0,m) u1[i-m+j]=(u1[i-m+j]-f[j]*u1[i])%mo; 
	}
	fo(i,0,m-1) z[i]=u1[i];
}
int main()
{
	cin>>n>>m;
	fo(i,1,m) scanf("%lld",&a[i]),f[m-i]=-a[i];
	f[m]=1;
	g[1]=1;
	s1[0]=1;
	for(int t=n;t;t>>=1)
	{
		if(t&1) mul(s1,g,s1);
		mul(g,g,g);
	}
	fo(i,0,m-1) scanf("%lld",&h[i]);
	LL ans=0;
	fo(i,0,m-1) ans=(ans+s1[i]*h[i]%mo+mo)%mo;
	printf("%lld\n",ans);
}
展开阅读全文

没有更多推荐了,返回首页