特征多项式实现线性递推学习笔记

用途

求解形如

\[g_i=\sum_{i=1}^kg_{n-i}a_i(i>k)\\ g_i=s_i(1\le i\le k) \]

的常系数线性齐次递推的\(g_n\)，其中\(n\)很大，\(k\)较小。

暴力

老师，我会暴力！

直接\(O(nk)\)暴力递推，没什么分。

老师，我会矩阵快速幂！

\(O(k^3\log n)\)，比上面的做法优一些，但\(k\le 2000\)时就上天了。

优化

最暴力的做法没什么营养，不管它，考虑优化矩阵快速幂的做法。

矩阵快速幂慢在了哪里呢？

设转移矩阵为\(A\)，复杂度瓶颈就在于在于要求出\(A^{n-1}\)。

设初始列向量\(s={[g_k,g_{k-1},\cdots,g_1]}^T\),\(A\times s\)即为\({[g_{k+1},g_{k},\cdots,g_2]}^T\)。

接下来引入一个概念：

特征多项式

设\(A\)为一矩阵，那么\(A\)的特征多项式就为\(f(\lambda)=|\lambda I-A|\)，其中\(I\)为单位矩阵，\(|A|\)表示\(A\)的行列式。

特征多项式的\(\lambda\)可以代入实数、多项式、矩阵……几乎所有东西。

回到正题。考虑\(A\)的特征多项式，可以发现

\[\lambda I-A=\left\{ \begin{matrix} \lambda-a_1&-a_2&-a_3&-a_4&\cdots&-a_m\\ -1&\lambda&0&0&\cdots&0\\ 0&-1&\lambda&0&\cdots&0\\ 0&0&-1&\lambda&\cdots&0\\ \cdots&\cdots&\cdots&\cdots&\ddots&\vdots\\ 0&0&0&0&-1&\lambda \end{matrix} \right\} \]

把它在第一行展开，用余子式计算行列式。

可以发现每一个余子式的主对角线都是满的，而且长得像是一个台阶的样子，可以很方便地削成上三角/下三角矩阵。

于是它的行列式就是主对角线的乘积。

于是推一波就会发现

\[f(\lambda)=\lambda^k-\sum_{i=0}^{k-1} a_{k-i}\lambda^i \]

同时，由\(f(\lambda)​\)的定义可知，\(f(A)=0​\)。

现在要求\(A^{n-1}\)。

设\(Q(A)=A^{n-1}\ \text{mod}\ f(A)\)，由\(f(A)=0\)可知\(A^{n-1}=Q(A)\)。

（也许你会觉得不能对\(0\)取模，但\(f(A)\)是一个多项式，对其取模相当于减去若干个\(f(A)\)，所以没有问题）

又因为\(Q(A)\)是不超过\(k\)次的多项式，所以我们大大减小了计算量，现在只需要求\(Q(A)\times G\)。

设

\[Q(A)=\sum_{i=0}^{k-1} c_iA^i \]

那么答案就是

\[\sum_{i=0}^{k-1}c_iA^iG \]

的第\(k\)行。

仔细思考，\(A^iG\)的第\(k\)行的意义是什么？

就是\(s_{i+1}\)！

于是答案就是

\[\sum_{i=0}^{k-1}c_is_{i+1} \]

大功告成。

复杂度瓶颈在多项式快速幂+取模，暴力\(k^2\log n\)，NTT、FFT优化可以到\(k\log k\log n\)。（加上\(\ln,\exp\)甚至可以\(O(k(\log k+\log n))\)？不过不能存下这么大的多项式，应该还是要快速幂+取模才行）

代码

这里以bzoj4161 Shlw loves matrixI为例。

注意该题中给出的是\(s_{0,\cdots, k-1}\)，需要对上方式子进行一定修改。

#include<bits/stdc++.h>
//clock_t t=clock();
namespace my_std{
	using namespace std;
	#define pii pair<int,int>
	#define fir first
	#define sec second
	#define MP make_pair
	#define rep(i,x,y) for (int i=(x);i<=(y);i++)
	#define drep(i,x,y) for (register int i=(x);i>=(y);--i)
	#define go(x) for (int i=head[x];i;i=edge[i].nxt)
	#define templ template<typename T>
	#define sz 4040
	#define mod 1000000007ll 
	typedef long long ll;
	typedef double db;
//	mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
//	templ inline T rnd(T l,T r) {return uniform_int_distribution<T>(l,r)(rng);}
	templ inline bool chkmax(T &x,T y){return x<y?x=y,1:0;}
	templ inline bool chkmin(T &x,T y){return x>y?x=y,1:0;}
	templ inline void read(T& t)
	{
		t=0;char f=0,ch=getchar();double d=0.1;
		while(ch>'9'||ch<'0') f|=(ch=='-'),ch=getchar();
		while(ch<='9'&&ch>='0') t=t*10+ch-48,ch=getchar();
		if(ch=='.'){ch=getchar();while(ch<='9'&&ch>='0') t+=d*(ch^48),d*=0.1,ch=getchar();}
		t=(f?-t:t);
	}
//	template<typename T,typename... Args>inline void read(T& t,Args&... args){read(t); read(args...);}
	char __sr[1<<21],__z[20];int __C=-1,__zz=0;
	inline void Ot(){fwrite(__sr,1,__C+1,stdout),__C=-1;}
	inline void print(register int x)
	{
		if(__C>1<<20)Ot();if(x<0)__sr[++__C]='-',x=-x;
		while(__z[++__zz]=x%10+48,x/=10);
		while(__sr[++__C]=__z[__zz],--__zz);__sr[++__C]='\n';
	}
	void file()
	{
		#ifndef ONLINE_JUDGE
		freopen("a.in","r",stdin);
		#endif
	}
//	inline void chktime()
//	{
//		#ifndef ONLINE_JUDGE
//		cout<<(clock()-t)/1000.0<<'\n';
//		#endif
//	}
	#ifdef mod
	ll ksm(ll x,int y){ll ret=1;for (;y;y>>=1,x=x*x%mod) if (y&1) ret=ret*x%mod;return ret;}
	ll inv(ll x){return ksm(x,mod-2);}
	#else
	ll ksm(ll x,int y){ll ret=1;for (;y;y>>=1,x=x*x) if (y&1) ret=ret*x;return ret;}
	#endif
//	inline ll mul(ll a,ll b){ll d=(ll)(a*(double)b/mod+0.5);ll ret=a*b-d*mod;if (ret<0) ret+=mod;return ret;}
}
using namespace my_std;

int n,K;
ll a[sz],s[sz];
ll f[sz];
ll c[sz],tmp[sz],tmp2[sz];

void mul(ll *a,ll *b,ll *ret)
{
	rep(i,0,K+K) tmp2[i]=0;
	rep(i,0,K-1) rep(j,0,K-1) (tmp2[i+j]+=a[i]*b[j])%=mod;
	drep(i,K+K-2,K)
		drep(j,K-1,0)	
			(tmp2[i-(K-j)]-=f[j]*tmp2[i])%=mod;
	rep(i,0,K-1) ret[i]=tmp2[i];
}

int main()
{
	file();
	read(n),read(K);
	rep(i,1,K) read(a[i]);
	rep(i,0,K-1) read(s[i]);
	f[K]=1;rep(i,0,K-1) f[i]=mod-a[K-i];
	tmp[1]=1;c[0]=1;
	for (int y=n;y;y>>=1,mul(tmp,tmp,tmp)) if (y&1) mul(tmp,c,c);
	ll ans=0;
	rep(i,0,K-1) (ans+=c[i]*s[i])%=mod;
	cout<<(ans+mod)%mod;
	return 0;
}

毒瘤洛谷\(k\le 32000\)不能暴力了……改天我再补个NTT吧。

update：常数巨大的代码走起~

（不开O2基本7000ms左右）

#include<bits/stdc++.h>
clock_t t=clock();
namespace my_std{
	using namespace std;
	#define pii pair<int,int>
	#define fir first
	#define sec second
	#define MP make_pair
	#define rep(i,x,y) for (int i=(x);i<=(y);i++)
	#define drep(i,x,y) for (int i=(x);i>=(y);i--)
	#define go(x) for (int i=head[x];i;i=edge[i].nxt)
	#define templ template<typename T>
	#define sz 300000
	#define mod 998244353ll 
	typedef long long ll;
	typedef double db;
	mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
	templ inline T rnd(T l,T r) {return uniform_int_distribution<T>(l,r)(rng);}
	templ inline bool chkmax(T &x,T y){return x<y?x=y,1:0;}
	templ inline bool chkmin(T &x,T y){return x>y?x=y,1:0;}
	templ inline void read(T& t)
	{
		t=0;char f=0,ch=getchar();double d=0.1;
		while(ch>'9'||ch<'0') f|=(ch=='-'),ch=getchar();
		while(ch<='9'&&ch>='0') t=t*10+ch-48,ch=getchar();
		if(ch=='.'){ch=getchar();while(ch<='9'&&ch>='0') t+=d*(ch^48),d*=0.1,ch=getchar();}
		t=(f?-t:t);
	}
	template<typename T,typename... Args>inline void read(T& t,Args&... args){read(t); read(args...);}
	char __sr[1<<21],__z[20];int __C=-1,__zz=0;
	inline void Ot(){fwrite(__sr,1,__C+1,stdout),__C=-1;}
	inline void print(register int x)
	{
		if(__C>1<<20)Ot();if(x<0)__sr[++__C]='-',x=-x;
		while(__z[++__zz]=x%10+48,x/=10);
		while(__sr[++__C]=__z[__zz],--__zz);__sr[++__C]='\n';
	}
	void file()
	{
		#ifndef ONLINE_JUDGE
		freopen("a.in","r",stdin);
		#endif
	}
	inline void chktime()
	{
		#ifndef ONLINE_JUDGE
		cout<<(clock()-t)/1000.0<<'\n';
		#endif
	}
	#ifdef mod
	ll ksm(ll x,int y){ll ret=1;for (;y;y>>=1,x=x*x%mod) if (y&1) ret=ret*x%mod;return ret;}
	ll inv(ll x){return ksm(x,mod-2);}
	#else
	ll ksm(ll x,int y){ll ret=1;for (;y;y>>=1,x=x*x) if (y&1) ret=ret*x;return ret;}
	#endif
//	inline ll mul(ll a,ll b){ll d=(ll)(a*(double)b/mod+0.5);ll ret=a*b-d*mod;if (ret<0) ret+=mod;return ret;}
}
using namespace my_std;

int limit,r[sz];
void NTT_init(int n)
{
    limit=1;int l=-1;
    while (limit<=n+n) limit<<=1,++l;
    rep(i,0,limit-1) r[i]=(r[i>>1]>>1)|((i&1)<<l);
}
void NTT(ll *a,int type)
{
    rep(i,0,limit-1) if (i<r[i]) swap(a[i],a[r[i]]);
    rep(i,0,limit-1) a[i]%=mod;
    for (int mid=1;mid<limit;mid<<=1)
    {
        ll Wn=ksm(3,(mod-1)/mid>>1);if (type==-1) Wn=inv(Wn);
        for (int j=0,len=mid<<1;j<limit;j+=len)
        {
            ll w=1;
            for (int k=0;k<mid;k++,w=w*Wn%mod)
            {
                ll x=a[j+k],y=a[j+k+mid]*w%mod;
                a[j+k]=(x+y)%mod;a[j+k+mid]=(x-y+mod)%mod;
            }
        }
    }
    if (type==1) return;
    ll I=inv(limit);
    rep(i,0,limit-1) a[i]=a[i]*I%mod;
}
ll tmp1[sz],tmp2[sz],tmp3[sz],tmp4[sz];
void PolyInv(ll *a,ll *f,int n) // f=a^{-1} mod x^n
{
    if (n==1) return (void)(f[0]=inv(a[0]));
    int mid=(n+1)>>1;
    PolyInv(a,f,mid);
    NTT_init(n);
    rep(i,0,mid-1) tmp1[i]=f[i];
    rep(i,0,n-1) tmp2[i]=a[i];
    NTT(tmp1,1);NTT(tmp2,1);
    rep(i,0,limit-1) tmp1[i]=tmp1[i]*(mod+2-tmp1[i]*tmp2[i]%mod)%mod;
    NTT(tmp1,-1);
    rep(i,0,n-1) f[i]=tmp1[i];
    rep(i,0,limit-1) tmp1[i]=tmp2[i]=0;
}
void Derivative(ll *a,ll *b,int n){rep(i,0,n-2) b[i]=a[i+1]*(i+1)%mod;b[n-1]=0;}
void Integrate(ll *a,ll *b,int n){drep(i,n-1,1) b[i]=a[i-1]*inv(i)%mod;b[0]=0;}
void PolyLn(ll *a,ll *f,int n) // f=ln a
{
    NTT_init(n);
    PolyInv(a,tmp3,n);Derivative(a,tmp4,n);
    NTT(tmp3,1);NTT(tmp4,1);
    rep(i,0,limit-1) tmp1[i]=tmp3[i]*tmp4[i]%mod;
    NTT(tmp1,-1);
    Integrate(tmp1,f,n);
    rep(i,0,limit-1) tmp1[i]=tmp3[i]=tmp4[i]=0;
}
void PolyExp(ll *a,ll *f,int n)
{
    if (n==1) return (void)(f[0]=1);
    int mid=(n+1)>>1;
    PolyExp(a,f,mid);
    rep(i,mid,n-1) f[i]=0;
    PolyLn(f,tmp2,n);
    rep(i,0,n-1) tmp1[i]=f[i];
    rep(i,0,n-1) tmp2[i]=(a[i]-tmp2[i]+mod)%mod;
    ++tmp2[0];
    NTT_init(n);
    NTT(tmp1,1);NTT(tmp2,1);
    rep(i,0,limit-1) tmp1[i]=tmp1[i]*tmp2[i]%mod;
    NTT(tmp1,-1);
    rep(i,0,n-1) f[i]=tmp1[i];
    rep(i,0,limit-1) tmp1[i]=tmp2[i]=0;
}
void PolyDiv(ll *a,ll *b,ll *c,ll *d,int n,int m) // a(n)=b(m)*c+d , c&d are the return values
{
    rep(i,0,n-m) tmp4[i]=(i>m?0:b[m-i]);
    PolyInv(tmp4,tmp3,n-m+1);
    rep(i,0,n-m) tmp4[i]=a[n-i];
    NTT_init(n-m+1);
    NTT(tmp3,1);NTT(tmp4,1);
    rep(i,0,limit-1) tmp3[i]=tmp3[i]*tmp4[i]%mod;
    NTT(tmp3,-1);
    reverse(tmp3,tmp3+n-m+1);
    rep(i,0,n-m) c[i]=tmp3[i];
    rep(i,0,m) tmp1[i]=b[i];
    rep(i,0,n) tmp2[i]=a[i];
    NTT_init(n);
    rep(i,n-m+1,limit-1) tmp3[i]=0;
    NTT(tmp1,1);NTT(tmp2,1);NTT(tmp3,1);
    rep(i,0,limit-1) tmp4[i]=(tmp2[i]-tmp1[i]*tmp3[i]%mod+mod)%mod;
    NTT(tmp4,-1);
    rep(i,0,m-1) d[i]=tmp4[i];
    rep(i,0,limit-1) tmp1[i]=tmp2[i]=tmp3[i]=tmp4[i]=0;
}

int n,K;
ll a[sz],s[sz],c[sz];
ll tmp[sz],f[sz];
ll t1[sz],t2[sz],t3[sz],t4[sz];
void mul(ll *a,ll *b,ll *ret)
{
	NTT_init(K+K+2);
	rep(i,0,K) t1[i]=a[i],t2[i]=b[i];
	NTT(t1,1);NTT(t2,1);
	rep(i,0,limit-1) t1[i]=t1[i]*t2[i]%mod;
	NTT(t1,-1);
	PolyDiv(t1,f,t4,t3,K+K,K);
	rep(i,0,K-1) ret[i]=t3[i];
	rep(i,0,limit-1) t1[i]=t2[i]=t3[i]=t4[i]=0;
}

int main()
{
    file();
    read(n,K);
    rep(i,1,K) read(a[i]);
    rep(i,0,K-1) read(s[i]);
    f[K]=1;rep(i,0,K-1) f[i]=mod-a[K-i];
    tmp[1]=1;c[0]=1;
    for (int y=n;y;y>>=1,mul(tmp,tmp,tmp)) if (y&1) mul(tmp,c,c);
    ll ans=0;
    rep(i,0,K-1) (ans+=c[i]*s[i])%=mod;
    cout<<(ans+mod)%mod;
    return 0;
}