常系数齐次线性递推优化矩阵快速幂-bzoj4161-4944

版权声明:欢迎转载(请附带原链接)ヾ(๑╹◡╹)ノ https://blog.csdn.net/corsica6/article/details/82696715


常系数齐次线性递推式

fk=i=1naifkif_k=\sum _{i=1}^{n} a_if_{k-i}

形如上式的dpdp转移式(ff表示dpdp状态,aa表示转移系数)即为常系数齐次线性递推式。对于这样的dpdp式,给定f1,2,..,n,a1,2,...,nf_{1,2,..,n},a_{1,2,...,n},即可递推出fk(k>n)f_k(k>n)

常系数:每一项均由前nn项乘以固定系数得到。齐次:都为一次式。线性递推:显然。

矩阵快速幂优化上式达到O(n3logk)O(n^3\log k)的复杂度(但这不是本文的重点)。

而常系数齐次线性递推式的转移矩阵的特征多项式可以O()O(读入)求出,所以通过特征多项式可以进一步优化到O(n2logk)O(n^2 \log k),再进一步优化多项式取模后可以做到O(nlognlogk)O(n\log n\log k )(本文讲解O(n2logk)O(n^2\log k)的暴力取模方法)。


特征值&特征向量&特征多项式

AAnn阶矩阵,如果存在常数λ\lambda非零向量v\vec v,使得Av=λvA\vec v=\lambda\vec v,则称λ\lambdaAA特征值(也称缩放系数),非零向量v\vec v为**AA的对应于特征值λ\lambda的特征向量**。

则:
Av=λv(λIA)v=0A\vec v=\lambda \vec v \Longleftrightarrow (\lambda I-A)\vec v=0

其有非零解的充要条件:det(λIA)=λIA=0det(\lambda I -A)=|\lambda I-A|=0
(det(A),Adet(A),|A|是矩阵AA的行列式的两种等价的表示方式)
(并不会证明,似乎det()=0det()=0的一个矩阵可以想象成拍扁成了向量?而一个可逆的(det()0det()\neq 0的矩阵乘上一个向量并不能得到0)。

方程λIA=0|\lambda I-A|=0称为AA特征方程

A=a(ij)A=a(ij)的秩(列秩/行秩)为nn(令AA是满秩的),则AA的有nn组线性不相关特征值和特征向量(唔,证明没有找到)。
则有:
(1)λ1+λ2+...+λn=a11+a22+...+ann(1)\lambda_1 +\lambda_2+...+\lambda_n=a_{11}+a_{22}+...+a_{nn}
(2)λ1λ2...λn=A(2)\lambda_1\lambda_2...\lambda_n=|A|
结合λIA=0|\lambda I-A|=0可以简单证明(理解)。

λIA=λn+k1λn1+...+kn1λ+kn|\lambda I-A|=\lambda^n+k_1\lambda^{n-1}+...+k_{n-1}\lambda +k^n这个λ\lambdann次多项式称为AA特征多项式ϕ(λ)\phi(\lambda)
也即ϕ(x)=xIA\phi(x)=|xI-A|
ϕ(λ)=0\phi(\lambda)=0nn个解即为AAnn个特征值。


特征矩阵&求解特征多项式

λIA=[λa11a12a1na21λa22a2nan1an2λann]\lambda I-A=\left [ \begin {matrix} \lambda -a_{11} & -a_{12} & \cdots & -a_{1n} \\ -a_{21} & \lambda -a_{22} &\cdots & -a_{2n} \\ \cdots & \cdots & \cdots & \cdots &\\ -a_{n1} & -a_{n2} &\cdots & \lambda -a_{nn}\\ \end {matrix} \right ]
该矩阵称为AA特征矩阵

回归到常系数齐次线性递推式fk=i=1naifkif_k=\sum _{i=1}^{n} a_if_{k-i}
[a1a2a3an1an1000001000001000000000010][fk1fk2fk3fk4fk(n1)fkn]=[fkfk1fk2fk3fk(n2)fk(n1)]\left [ \begin {matrix} a_1 &a_2 & a_3& \cdots &a_{n-1} &a_n\\ 1 & 0&0& \cdots & 0 & 0\\ 0& 1 & 0 &\cdots &0&0\\ 0 & 0&1& \cdots &0&0\\ \cdots &\cdots &\cdots &\cdots &\cdots &\cdots &\\ 0&0&0&\cdots &0 &0\\ 0&0&0&\cdots &1&0 \\ \end {matrix} \right] \left [ \begin{matrix} f_{k-1}\\f_{k-2}\\ f_{k-3}\\f_{k-4}\\ \vdots\\f_{k-(n-1)} \\f_{k-n} \end{matrix} \right ] = \left [ \begin{matrix} f_{k}\\f_{k-1}\\ f_{k-2}\\f_{k-3}\\ \vdots\\f_{k-(n-2)} \\f_{k-(n-1)} \end{matrix} \right ]
构造转移矩阵的特征矩阵:
ϕ(λ)=λIA=[λa1a2a3an1an1λ00001λ0000100000λ00001λ] \phi(\lambda)=|\lambda I-A|= \left [ \begin {matrix} \lambda -a_1 & -a_2 & -a_3& \cdots &-a_{n-1} &-a_n\\ -1 & \lambda&0& \cdots & 0 & 0\\ 0& -1 & \lambda &\cdots &0&0\\ 0 & 0&-1& \cdots &0&0\\ \cdots &\cdots &\cdots &\cdots &\cdots &\cdots &\\ 0&0&0&\cdots &\lambda &0\\ 0&0&0&\cdots &-1&\lambda \\ \end {matrix} \right]
这里可以强行手算一下行列式(全排列逆序对方法),因为这个矩阵每列最多只有三行不为0,所以可以模拟/枚举一下:

  • 选中第一列第一排的λa1\lambda -a_1,由于每排每列只能选一个,且不能选0,则只有一种方案:沿对角线选择。贡献为:λn1(λa1)\lambda^{n-1}(\lambda -a_1)
  • 选中第一列第二排的1-1,可以想到在有一列选中第一排之前,每一列ii都只能选第i+1i+1排的数1(第ii排被上一列占了),选了第一排后,为保证有解,剩下的列都只能选对角线上的数,所以结合逆序对,总贡献为:i=2n(ai)(1)i1λni(1)i1=i=2naiλni\sum_{i=2}^n(-a_i)(-1)^{i-1}\lambda^{n-i}(-1)^{i-1}=-\sum_{i=2}^n a_i\lambda ^{n-i}

综上:
(1)ϕ(λ)=λIA=λni=1naiλni\phi(\lambda)=|\lambda I-A|=\lambda ^{n}-\sum_{i=1}^n a_i\lambda ^{n-i}\tag 1
这样就得到了特征多项式的每项系数。


常系数齐次线性递推优化

ϕ(A)=AIA=0\phi(A)=|AI-A|=0得到ϕ(A)\phi(A)为0矩阵。

但似乎上面这个证明是错的?

前文已经说明ϕ(x)\phi(x)是一个nn次多项式,同样由ϕ(λ)=0\phi(\lambda)=0可以用插值法表达为:ϕ(x)=i=1n(λix)\phi(x)=\prod \limits_{i=1}^n (\lambda_i -x)

ϕ(A)=i=1n(λiIA)\phi(A)=\prod \limits_{i=1}^n (\lambda_iI-A)。显然对于所有1in1\leq i\leq n都有ϕ(A)vi=0\phi(A)\vec v_i=0,又因为vi\vec v_i互相线性不相关,所以AA变换得到的ϕ(A)\phi(A)矩阵为0矩阵。得证。

假设原本需要矩阵快速幂求得AkA^{k},现在去掉无用的00矩阵,答案是相同的,实际上就是多项式取模的过程,最终得到一个最高次幂为n1n-1次的多项式g(A)g(A)
(2)g(A)=Ak%ϕ(A)=i=0n1piAig(A)=A^k\%\phi(A)=\sum _{i=0}^{n-1} p_i A^i\tag 2

设初始递推矩阵为:f[1]=fn,fn1,...,f1f[1]=f_n,f_{n-1},...,f_1,则最终得到:f[1]Ak=f[1]g(A)=f[1+k]f[1]A^k=f[1]g(A)=f[1+k],进一步转化为f[1+k]=i=0n1pif[1+i]f[1+k]=\sum\limits_{i=0}^{n-1}p_if[1+i],但实际上我们只需要f[1+k]f[1+k]中的最后一项f1+kf_{1+k},可以将矩阵加法展开,故:
(3)f1+k=i=0n1pif1+if_{1+k}=\sum\limits_{i=0}^{n-1}p_if_{1+i}\tag 3

常系数齐次线性递推优化的过程即为:

(1)(1)求出转移矩阵AA的特征多项式ϕ(A)\phi(A)

(2)(2)ϕ(A)\phi(A)为模做AA的快速幂得到AkA^k

(3)(3)解出n1n-1次多项式g(A)g(A)的每项系数转移得到答案


暴力多项式取模

这里补充一点知识:

关于多项式取模法则:
如:5x5+3x4+x3+x2+6x+1mod  (x3+x2+x+1)5x^5+3x^4+x^3+x^2+6x+1 \mod (x^3+x^2+x+1)
我们先把 (x3+x2+x+1)(x^3+x^2+x+1)乘以 x2x^2,把被除的多项式中的 5x55x^5 消掉(做减法)
然后以此把次高位消掉,直到消到x3x^3为止。
来源

这是一个O(n2)O(n^2)的暴力取模。

O(nlogn)O(n\log n)的取模看这里


模板

bzoj4161

一道真·模板

#include<bits/stdc++.h>
#define gc getchar()
#define si isdigit(ch)
#define RI register
using namespace std;

const int N=4002,mod=1e9+7;
int ans,n,K,a[N],b[N],c[N],e[N],re[N];

char ch; 
inline void rd(int &x)
{
	ch=gc;x=0;int f=0;
	for(;!si;ch=gc) if(ch=='-') f=1;
	for(;si;ch=gc) x=x*10+(ch^48);
	if(f) x=mod-x;
} 

inline int ad(int x,int y){x+=y;return x>=mod? x-mod:x;}
inline int dc(int x,int y){x-=y;return x<0? x+mod:x;}
inline int mul(int x,int y){return 1ll*x*y%mod;}

inline void Mul(int *f,int *g,int n)
{
	RI int i,j;
	memset(re,0,sizeof(int)*(n+1));
	for(i=0;i<=K;++i)
	 for(j=0;j<=K;++j)
	  re[i+j]=ad(re[i+j],mul(f[i],g[j]));
	for(i=n;i>=K;--i){
		if(!re[i]) continue;
		for(j=0;j<K;++j)
		 re[i-K+j]=ad(re[i-K+j],mul(re[i],a[K-j]));
		re[i]=0;
	}
	for(i=0;i<=K;++i) g[i]=re[i];
}

int main(){
	RI int i,j,M;
	rd(n);rd(K);M=K<<1;
	for(i=1;i<=K;++i) rd(a[i]);
	for(i=1;i<=K;++i) rd(b[i]);
	if(n<K){printf("%d\n",b[n+1]);return 0;}
	e[0]=c[1]=1;
	for(;n;n>>=1,Mul(c,c,M))
	 if(n&1) Mul(c,e,M);
	for(i=0;i<K;++i) ans=ad(ans,mul(e[i],b[i+1]));
	printf("%d\n",ans);
}

bzoj4944 noi2017泳池

#include<bits/stdc++.h>
#define RI register
using namespace std;
const int mod=998244353;
const int N=2005;
int a[N],e[N],c[N],b[N],f[2][N],ss[N];
int pw[N],rvq,q,n,lim,xx,yy,ans;

inline int ad(int x,int y){x+=y;return x>=mod? x-mod:x;}
inline int dc(int x,int y){x-=y;return x<0? x+mod:x;}
inline int mul(int x,int y){return 1ll*x*y%mod;}

inline int fp(int x,int y)
{
	int re=1;
	for(;y;y>>=1,x=mul(x,x))
	 if(y&1) re=mul(re,x);
	return re;
}

inline void Mul(int *f,int *g,int n,int K)
{
	RI int i,j,k,t;
	memset(ss,0,sizeof(int)*(n+1));
	for(i=0;i<K;++i)
	 for(j=0;j<K;++j)
	  ss[i+j]=ad(ss[i+j],mul(f[i],g[j]));
	for(i=n;i>=K;--i){
		if(!ss[i]) continue;
		for(j=1;j<=K;++j)
		 ss[i-j]=ad(ss[i-j],mul(ss[i],b[j]));
		ss[i]=0;
	}
	for(i=0;i<=K;++i) g[i]=ss[i];
}

inline int cal(int n,int lim)
{
	RI int i,j,k,t,ori;
	memset(f[1],0,sizeof(f[1]));
	int *F=f[0],*G=f[1];G[0]=1;
	for(i=lim-1;i;--i,swap(F,G)){
		memset(F,0,sizeof(int)*(lim+1));
		for(j=0;i*j<lim;++j){
			F[j]=mul(pw[j],G[j]);
			for(k=0;k<j;++k)
			  F[j]=ad(F[j],mul(pw[k],mul(G[k],mul(rvq,F[j-k-1]))));
		}
	}
	for(b[1]=rvq,i=1;i<lim;++i)
	 b[i+1]=mul(mul(G[i],pw[i]),rvq);
    for(a[0]=1,i=1;i<lim;++i){
    	a[i]=mul(rvq,a[i-1]);
    	for(j=1;j<min(i,lim);++j)
    	    a[i]=ad(a[i],mul(b[j+1],a[i-1-j]));
    	if(i<lim) a[i]=ad(a[i],mul(pw[i],G[i]));
    }
    t=lim<<1;ori=lim;
    memset(e,0,sizeof(int)*(t+1));
    memset(c,0,sizeof(int)*(t+1));
    e[0]=c[1]=1;
    for(;n;n>>=1,Mul(c,c,t,ori))
     if(n&1) Mul(c,e,t,ori);
	for(ans=0,i=0;i<lim;++i)
	 ans=ad(ans,mul(e[i],a[i]));
	return ans;
} 
 

int main(){
    scanf("%d%d%d%d",&n,&lim,&xx,&yy);
	q=mul(xx,fp(yy,mod-2));rvq=dc(1,q);
	pw[0]=1;
	for(RI int i=1;i<=lim+1;++i) pw[i]=mul(pw[i-1],q); 
	printf("%d\n",dc(cal(n,lim+1),cal(n,lim)));
}

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