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

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

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

\[\text{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})\)，则：

\[e\cdot 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,j\in[1,m] \\ a_{i,i+1}=1,i\in[1,m-1] \]

大概是这样的

\[\begin{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,i\le m\)的和表示出来。 这意味着我们可能可以利用指数有结合律的优美性质。

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

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

然而，\(A\)是一个\(m\times m\)的方阵，\(m\le 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(m\log m\log 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;
}