【瞎讲】 Cayley-Hamilton 常系数齐次线性递推式第n项的快速计算 (m=1e5,n=1e18)

【背诵瞎讲】 Cayley-Hamilton 常系数齐次线性递推式第n项的快速计算 (m=1e5,n=1e18)

看CSP看到一题“线性递推式”，不会做，去问了问zsy怎么做，他并不想理我并丢给我以下方法：

[ ext{Cayley-Hamilton} ]

下文会根据CH定理证明的思路证明，没有形式上使用特征系统，因为我也不会...

复读鸽子讲的话

一句话就是求：

[f_n=sum_{i=1}^m c_if_{n-i} mod 998244353 ]

但这个算法卡常,zsy说1e5估计要跑10s

先求前(m)项

前(m)项的话，递推关系就变成这样：

[f_i=sum_{j=1}^i c_jf_{i-j} ]

钦定(c_0=0)，则写成生成函数(F(x)=sum f_{i+1}x^i,C(x)=sum c_ix^i)，则

[F(x)-f_1=C(x)F(x) ]

为什么有个(f_1)？因为你会发现按照定义(F(x))少了常数项((c_0=0) )

解出来

[G(x)=dfrac {f_1} {1-C(x)} ]

就是：

[G(x)=f_1(1-C(x))^{-1} ]

所以直接就求出来了前(m)项

增广矩阵的特征多项式

我们规范一下平常用的矩阵快速幂的格式：

设(m)阶增广矩阵(A)，我们所谓的增广是这样的：

对于行向量(e=(f_1,f_2,dots,f_{m-1},f_{m}))，则：

[ecdot A^n=(f_{1+n},f_{2+n},f_{3+n},f_{4+n},dots ,f_{m+n}) ]

我们就是要求(f_{1+n})。归纳地理解上面这个式子，这是我们对(A)的定义，根据递推关系很容易构造出一个满足条件的(A)

[a_{j,m}=c_j,jin[1,m] \ a_{i,i+1}=1,iin[1,m-1] ]

大概是这样的

[egin{pmatrix} 0&0&0&0&dots&a_1 \ 1&0&0&0&dots&a_2 \ 0&1&0&0&dots&a_3 \ 0&0&1&0&dots&a_4 \ vdots&vdots&vdots&vdots&ddots&vdots \ 0&0&0&0&1&a_m end{pmatrix} ]

我们新建一个生成函数(F(x))(和前面那个没有关系，上面那个是一个“局部变量”hhh)

[F(x)=x^m-sum_{i=1}^{m}c_{i}x^{m-i} ]

这个多项式叫做A的特征多项式，我们带入(A)进去有：

[F(A)=A^m-sum_{i=1}^{m}c_{i}A^{m-i}= O ]

如何证明？先移项再左乘一个行向量(e),定义在上面

[eA^m=sum_{i=1}^{m}c_{i}eA^{m-i} ]

根据(e)和(A)之间的定义，直接变成

[(f_{1+m},f_{2+m},dots,f_{m+m})=sum_{i=1}^{m} c_{i}(f_{1+m-i},f_{2+m-i},dots,f_{m+m-i}) ]

考虑一下第(j)维的值，两边分别为

[f_{j+m}=sum_{i=1}^{m} c_{i} f_{j+m-i} ]

显然成立。故对于每一维相等，故左右两个向量相等。

故原命题成立。

好现在我们就证明了

[A^m=sum_{i=1}^m c_iA^{m-i} ]

这其实是我们的"初始条件"

假设会求(A^n)了，我们怎么得到(f_n)?

请仔细阅读标题

我们之前得到了:

[A^m=sum_{i=1}^m c_iA^{m-i} ]

考虑我们实质上干了什么？实际上我们就把下标换成了指数。也就是说，我们可以把(A)的任意次方化为若干(A^i,ile m)的和表示出来。 这意味着我们可能可以利用指数有结合律的优美性质。

形式化的，我们就得到了这样一个(b[])数组，使得：

[A^n=sum_{i=1}^m b_iA^{m-i} ]

然而，(A)是一个(m imes m)的方阵，(mle 1e5)。我们是不可能对于(A)在程序中进行任何矩阵(O(n^2),O(n^3))的操作的，否则超时(理解我的意思就好)。

但是我们只是要求(f_n)，所以两边直接同右乘以(e)就好了。

就是

[eA^n=sum_{i=1}^m b_ieA^{m-i} ]

然后根据定义，我们就直接得到了这个：

[(f_{1+n},dots,f_{m+n})=sum_{i=1}^m b_i (f_{1+m-i},dots,f_{m+m-i}) ]

然而我们只需要第一个(f_{1+n})

所以把第一维单独拿出来：

[f_{1+n}=sum_{i=1}^m b_if_{1+m-i} ]

矩阵已经消失了，剩下的只有一个(O(m))的式子。

现在问题就变成了要求(b[])数组，然后就直接(O(m))算(f_{1+n})了

那么到底如何球(b[])数组呢？

关键的系数(b_i)

现在的问题变成了如何得到(b_i)数组

递归地思考，我们考虑任何时刻用他的系数(b[])数组，不及(m)位的补(0)。例如：(A^1->b[]=underbrace {{0,0,0,underbrace {1}_{m-1\_th},dots}}_{m})来表示我们当前的状态。

两个(A^1)相乘就构成(A^2)，我们就可以倍增求我们要的(A^n) 。

然而每次相乘会出现(2m-2)次方，不符合要求啊，怎么办？但是我们有个关系式

[x^m=sum_{i=1}^mc_ix^{m-i} ]

代入我们之前得到的(2m-2)次方多项式里，所有次数大于(m)的项都会被化为小于(m)的形式，现在就是要快速带入这个关系式

最神奇的事情发生了！！

带入

[x^m=sum_{i=1}^mc_ix^{m-i} ]

就相当于认为

[x^m-sum_{i=1}^mc_ix^{m-i}=0 ]

也就是对(x^m-sum_{i=1}^mc_ix^{m-i})多项式取膜!!1

很拍脑袋的解释是吧，其实是有科学说明的:

考虑取模那个定义的式子，可以发现是F=MOD*Q+R，R是余数，由于我们带入前式也是带入MOD=0进去，现在就很显然啦

于是就一边从(A^1)对应的那个系数倍增，一边倍增一边对(x^m-sum_{i=1}^mc_ix^{m-i})多项式取膜就好了！！1。

多项式取膜一个log，加上倍增，总复杂度(O(mlog mlog n))，常数真 的 好 大 啊

orz-orz--orz--orz--orz--orz--orz--orz--orz--orz--orz--orz-分鸽线-orz--orz--orz--orz--orz--orz--orz--orz--orz--orz--orz--orz--orz--orz--orz--orz-

代码

在路上了在路上了

觚不觚，觚哉觚哉！

//@winlere
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<assert.h>
#define getchar() (__c==__ed?(__ed=__buf+fread(__c=__buf,1,1<<18,stdin),*__c++):*__c++)

using namespace std;  typedef long long ll;   char __buf[1<<18],*__c=__buf,*__ed=__buf;
inline int qr(){
	int ret=0,f=0,c=getchar();
	while(!isdigit(c))f|=c==45,c=getchar();
	while(isdigit(c)) ret=ret*10+c-48,c=getchar();
	return f?-ret:ret;
}
const int mod=998244353;
const int maxn=1<<16;
inline int MOD(const int&x){return x>=mod?x-mod:x;}
inline int MOD(const int&x,const int&y){return 1ll*x*y%mod;}
inline int ksm(const int&ba,const int&p){
	int ret=1;
	for(int t=p,b=ba;t;t>>=1,b=MOD(b,b))
		if(t&1) ret=MOD(ret,b);
	return ret;
}

namespace poly{
	const int g=3,gi=(mod+1)/3;
	void NTT(int*a,const int&len,const int&tag){
		static int r[maxn];
		for(int t=0;t<len;++t)
			if((r[t]=r[t>>1]>>1|(t&1?len>>1:0))<t)
				swap(a[t],a[r[t]]);
		for(int t=1,wn,s=tag==1?g:gi;t<len;t<<=1){
			wn=ksm(s,(mod-1)/(t<<1));
			for(int i=0;i<len;i+=t<<1)
				for(int j=0,w=1,p;j<t;++j,w=MOD(w,wn))
					p=MOD(w,a[i+j+t]),a[i+j+t]=MOD(a[i+j]-p+mod),a[i+j]=MOD(a[i+j]+p);
		}
		if(tag!=1)
			for(int t=0,i=mod-(mod-1)/len;t<len;++t)
				a[t]=MOD(a[t],i);
	}
	void INV(int*a,int*b,const int&len){
		if(len==1) return b[0]=ksm(a[0],mod-2),void();
		INV(a,b,len>>1);
		static int A[maxn],B[maxn];
		memset(A,0,len<<3); memset(B,0,len<<3);
		memcpy(A,a,len<<2); memcpy(B,b,len<<2);
		NTT(A,len<<1,1); NTT(B,len<<1,1);
		for(int t=0;t<len<<1;++t) A[t]=MOD(A[t],MOD(B[t],B[t]));
		NTT(A,len<<1,0);
		for(int t=0;t<len;++t) b[t]=MOD(MOD(B[t]+B[t])-A[t]+mod);
	}
	int sav[maxn],invsav[maxn],n,m;//2^
	void MOD(int*a){
		static int A[maxn],B[maxn];
		int len=n;
		memset(A,0,len<<3); memset(B,0,len<<3);
		memcpy(A,a,len<<2); reverse(A,A+len);
		NTT(A,len<<1,1);
		for(int t=0;t<len<<1;++t) A[t]=::MOD(A[t],invsav[t]);
		NTT(A,len<<1,0);
		reverse(A,A+len); memset(A+len,0,len<<2);
		NTT(A,len<<1,1);
		for(int t=0;t<len<<1;++t) A[t]=::MOD(A[t],sav[t]);
		NTT(A,len<<1,0);
		for(int t=0;t<=m;++t) a[t]=::MOD(a[t]-A[t]+mod);
		//F*G
	}
	void KSM(int*a,int*Mod,int p){
		static int ret[maxn],base[maxn];
		memset(ret,0,sizeof ret); memset(base,0,sizeof base);
		ret[0]=1;
		for(int t=0;t<n;++t) base[t]=a[t];
		while(p){
			p>>=1;
		}
	}
}

int main(){
#ifndef ONLINE_JUDGE
	freopen("in.in","r",stdin);
	freopen("out.out","w",stdout);
#endif
      
	return 0;
}