知识点简单总结——常系数齐次线性递推
知识点简单总结——常系数齐次线性递推
Cayley-Hamilton定理
参考自：https://www.cnblogs.com/regenwald/articles/7804744.html
令 $ A $ 为一个 $ n \times n $ 阶的矩阵， $ p( \lambda ) $ 为其特征方程，即 $ \det{ ( A - \lambda I ) } $ 。
则有 $ p( A ) = 0 $ 。证明参见上述博客。
之后证明会用到这个定理。
常系数齐次线性递推
给定 $ f_{0...k-1} , c_{1...k} $ ，求 $ f_{n} = \sum\limits_{ i = 1 }^{ k } f_{n-i} c_{i} $ ， $ n \le 10^{9} , k \le 32000 $ 。
很明显有 $ O(nk) $ 的暴力和 $ O( k^{ 3 } \log { n } ) $ 的矩阵乘法。
考虑下矩阵乘法的实现过程，
设转移矩阵为 $ A $ ，初始向量为 $ S $ 。
毫无疑问对于只需要求一个 $ f_{n} $ 的情况下矩阵里的东西太多了，考虑省掉一部分。
考虑构造序列 $ g $ 有 $ A^{ n } = \sum\limits_{ i = 0 }^{ k - 1 } g_{ i } A^{ i } $ 。
那么由于仅仅需要求出数列一项的答案，我们直接把上式改成求
进而答案就是
...不对啊，这 $ ( S \times A^{ i } )_{ 0 } $ ...
不就是 $ S_{i} $ 。
好了答案出来了： $ ans = \sum\limits_{ i = 0 }^{ k - 1 } g_{ i } S_{ i } $ 。
准备考虑求 $ g $ 数组。
转移是一个矩阵乘一个向量得到一个新向量，考虑写出特征方程然后变为特征多项式：
由Cayley-Hamilton定理有 $ H(A) = 0 $ ， $ 0 $ 指 $ 0 $ 矩阵。
集合上述结论考虑将求 $ A^{N} \ mod \ H(A) = G(A) $ 。
很巧的求出结果就是要求的 $ g $ 数组。
上式直接当 $ x^{N} $ 的多项式取模就好。
（是不是有点细节没说清诶。。。算了咕了吧反正主体思路到位了（被打））
最终时间复杂度为 $ O( k \cdot \log{ k } \cdot \log{ m } ) $ 。
多项式取模注意预处理下模数数组的的逆和其反转数组的逆，不然会T死（
代码常数略大，开 $ O2 $ 在 $ 1s $ 左右。
#include<bits/stdc++.h>
using namespace std;
typedef long long lint;
template<typename TP>inline void read(TP &tar)
{
TP ret=0,f=1;char ch=getchar();
while(ch<'0'||ch>'9'){if(ch=='-') f=-1;ch=getchar();}
while(ch>='0'&&ch<='9'){ret=ret*10+ch-'0';ch=getchar();}
tar=ret*f;
}
template<typename TP,typename... Args>inline void read(TP& t,Args&... args){read(t),read(args...);}
namespace RKK
{
const int N=140011,maxn=131072;
const int mo=998244353,G=3;
void doadd(int &a,int b){if((a+=b)>=mo) a-=mo;}int add(int a,int b){return (a+=b)>=mo?a-mo:a;}
void domul(int &a,int b){a=1ll*a*b%mo;}int mul(int a,int b){return 1ll*a*b%mo;}
int fpow(int a,int p){int ret=1;while(p){if(p&1) domul(ret,a);domul(a,a),p>>=1;}return ret;}
int wg[N],iwg[N];int rev[N],lastlen;
inline void init(){for(int i=1;i<maxn;i<<=1) wg[i]=fpow(G,(mo-1)/(i<<1)),iwg[i]=fpow(wg[i],mo-2);}
inline void getrev(int len){for(int i=1;i<len;i++)rev[i]=(rev[i>>1]>>1)|((i&1)*(len>>1));}
inline void ntt(int *a,int len,int tp)
{
if(lastlen!=len) getrev(len),lastlen=len;
for(int i=0;i<len;i++)if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int i=1;i<len;i<<=1)
{
int w0=(~tp)?wg[i]:iwg[i];
for(int j=0;j<len;j+=(i<<1))
{
int w=1;
for(int k=0;k<i;k++,domul(w,w0))
{
int w1=a[j+k],w2=mul(w,a[j+k+i]);
a[j+k]=add(w1,w2),a[j+k+i]=add(w1,mo-w2);
}
}
}
int ilen=fpow(len,mo-2);
if(tp==-1)for(int i=0;i<len;i++) domul(a[i],ilen);
}
namespace poly
{
void clear(int *a,int n){memset(a,0,n*4);}
void copy(int *a,int *b,int n){memcpy(a,b,n*4);}
namespace inv
{
int a[N],b[N];
void work(int *f,int *g,int n)
{
g[0]=fpow(f[0],mo-2);
for(int len=1;len<n<<1;len<<=1)
{
copy(a,f,len),copy(b,g,len);
ntt(a,len<<1,1),ntt(b,len<<1,1);
for(int i=0;i<len<<1;i++) g[i]=mul(b[i],(2ll-mul(a[i],b[i])+mo));
ntt(g,len<<1,-1),clear(g+len,len);
}
clear(a,n*2),clear(b,n*2);
}
}
}
int t[N];
int n,m,ans,len;
int f[N],g[N],a[N],b[N],c[N],ic[N];
void mod(int *a,int n,int m)
{
int nn=n-m+1;for(int i=0;i<nn;i++) t[i]=a[n-1-i];
ntt(t,len,1);for(int i=0;i<len;i++) domul(t[i],ic[i]);ntt(t,len,-1);
reverse(t,t+nn),poly::clear(t+nn,len-nn);
ntt(t,len,1);for(int i=0;i<len;i++) domul(t[i],c[i]);ntt(t,len,-1);
for(int i=0;i<m-1;i++) a[i]=add(a[i],mo-t[i]);poly::clear(a+m-1,len-(m-1));
poly::clear(t,len);
}
int main()
{
init();
read(n,m);for(int i=1;i<=m;i++) read(c[i]),doadd(c[i],mo);
for(int i=0;i<m;i++) read(f[i]),doadd(f[i],mo);if(n<m) return printf("%d\n",f[n]),0;
len=1;while(len<=m) len<<=1;
reverse(c+1,c+1+m);for(int i=0;i<m;i++) c[i]=add(0,mo-c[i+1]);c[m]=1;
for(int i=0;i<=m;i++) b[i]=c[m-i];poly::inv::work(b,ic,len);
g[0]=1;a[1]=1;len<<=1;
ntt(c,len,1),ntt(ic,len,1);
while(n)
{
if(n&1)
{
ntt(g,len,1),ntt(a,len,1);for(int i=0;i<len;i++) domul(g[i],a[i]);ntt(g,len,-1),ntt(a,len,-1);
mod(g,len,m+1);
}
ntt(a,len,1);for(int i=0;i<len;i++) domul(a[i],a[i]);ntt(a,len,-1);
mod(a,len,m+1);
n>>=1;
}
for(int i=0;i<m;i++) doadd(ans,mul(f[i],g[i]));printf("%d\n",ans);
return 0;
}
}
int main(){return RKK::main();}
应用
不是很明显了嘛。
