P4723 【模板】常系数齐次线性递推
题目描述
求一个满足$k$阶齐次线性递推数列${a_i}$的第$n$项。
即：$a_n=sumlimits_{i=1}^{k}f_i imes a_{n-i}$
输入输出格式
输入格式：第一行两个数$n$,$k$，如题面所述。
第二行$k$个数，表示$f_1 f_2 cdots f_k$
第三行$k$个数，表示$a_0 a_1 cdots a_{k-1}$
输出格式：一个数，表示 $a_n \% 998244353$ 的值
输入输出样例
说明
$N = 10^{9} , K = 32000$
题解
先修《数学选修4-2：矩阵与变换》。
常系数线性递推
给出递推式(f_n=sum_{i=1}^ka_if_{n-i})，和初值条件(f_1,f_2,dots ,f_k)，求(f_n)。
可以用生成函数解一下，然后多项式求逆，(O(nlog n))。当然这对于(n=10^9)的数据范围是不行的。
矩阵快速幂解法
由递推式，构造矩阵和列向量
计算(A^{n-k}F)即可，(O(k^3log n))。然而这对于(k=32000)的数据范围还是不行。
矩阵的特征值和特征向量
(A^nF)是线性变换的形式，自然也可以用特征值与特征向量来求解。
特征多项式是(f_A(lambda)=|A-lambda I|)，于是特征方程为
手动高斯消元，消成上三角可得
然后呢，把特征值解出来然后用特征向量那套理论吗？虽然这种方法可行，但是只能做低次的（4次及以下），所以我们需要新科技。
矩阵的多项式
对于(n)次多项式(f(x))，将矩阵(A)看做自变量带入，得
记(f(A))为(A)的(n)次多项式。与另一个(A)的(m)次多项式(g(A))，其乘法运算满足交换律，即
Cayley-Hamilton 定理
特征多项式(f_A(x)=sum_{i=0}^{k-1}a_{k-i}x^i-x^k)，这里注意下标。
定理：(f_A(A)=0)，即矩阵被自己的特征多项式化零。记忆方法：(f_A(A)=|A-AI|=0)。
推论：(A^n=q(A)f_A(A)+r(A)=r(A))，其中(r(A)=A^nmod f_A(A))。
所以(A^nF=r(A)F=sum_{i=0}^{k-1}r_iA^iF)，(O(k^4))。？？？
多项式优化
将Cayley-Hamilton定理的推论变成普通多项式：(x^n=q(x)f_A(x)+r(x))，其中(r(x)=x^nmod f_A(x))。
于是我们可以对多项式(x)做快速幂并取模(f_A(x))，即可得出(r(x))的系数，即(r(A))的系数，(O(klog klog n))。
而(A^iF=egin{bmatrix}f_{1+i} & f_{2+i} & dots & f_{k+i}end{bmatrix}^T)，所以(f_{n+k}=sum_{i=0}^{k-1}r_if_{k+i})。
问题转化成了如何求(f)的前(2k)项。这时使用生成函数和多项式求逆，(O(klog k))。
于是我们便得到了(O(klog klog n+klog k+k))的优秀解法。
小优化
刚才说的快速幂是指(A^{n-k})，如果我们做到(A^{n-1})，就只需要前(k)项了，那么就用不着什么生成函数了。(O(klog klog n+k))。
UPD：然后我发现(a_i)系数要翻转，还要取相反数，是我把特征方程写反了？的确数学书上的应该是(|lambda I-A|)，但这没有影响。我看差异是我和其他人转移矩阵列的不一样，我把他们的做了一个中心对称（这是讲义的风格），然后特征方程都不一样了？我觉得这不是我以现有的水平能搞明白的问题。但是大众的写法的优势在于他们的最高次项系数为1，在做暴力多项式取模的时候很好办。
UPD：最近我发现最开始的那个矩阵含有 (a) 的最后一行写反了……
void num_trans(polynomial&a,int inverse){
int limit=a.size(),len=log2(limit);
static vector<int> bit_rev;
if(bit_rev.size()!=limit){
bit_rev.resize(limit);
for(int i=0;i<limit;++i) bit_rev[i]=bit_rev[i>>1]>>1|(i&1)<<(len-1);
}
for(int i=0;i<limit;++i)if(i<bit_rev[i]) swap(a[i],a[bit_rev[i]]);
for(int step=1;step<limit;step<<=1){
int gn=fpow(inverse==1?g:g_inv,(mod-1)/(step<<1));
for(int even=0;even<limit;even+=step<<1){
int odd=even+step,gk=1;
for(int k=0;k<step;++k,gk=mul(gk,gn)){
int t=mul(gk,a[odd+k]);
a[odd+k]=add(a[even+k],mod-t),a[even+k]=add(a[even+k],t);
}
}
}
if(inverse==-1){
int lim_inv=fpow(limit,mod-2);
for(int i=0;i<limit;++i) a[i]=mul(a[i],lim_inv);
}
}
polynomial poly_inv(polynomial a,int n){ // mod x^n
polynomial b(1,fpow(a[0],mod-2));
if(n==1) return b;
int limit;
for(limit=2;limit<n;limit<<=1){
polynomial a1(a.begin(),a.begin()+limit);
a1.resize(limit<<1),num_trans(a1,1);
b.resize(limit<<1),num_trans(b,1);
for(int i=0;i<limit<<1;++i) b[i]=mul(add(2,mod-mul(a1[i],b[i])),b[i]);
num_trans(b,-1),b.resize(limit);
}
a.resize(limit<<1),num_trans(a,1);
b.resize(limit<<1),num_trans(b,1);
for(int i=0;i<limit<<1;++i) b[i]=mul(add(2,mod-mul(a[i],b[i])),b[i]);
num_trans(b,-1),b.resize(n);
return b;
}
polynomial poly_div(polynomial f,polynomial g){ // return the quotient
int n=f.size()-1,m=g.size()-1;
reverse(g.begin(),g.end()),g.resize(n-m+1),g=poly_inv(g,n-m+1);
reverse(f.begin(),f.end()),f.resize(n-m+1);
int limit=1<<int(ceil(log2(2*(n-m)+1)));
f.resize(limit),g.resize(limit);
num_trans(f,1),num_trans(g,1);
for(int i=0;i<limit;++i) f[i]=mul(f[i],g[i]);
num_trans(f,-1),f.resize(n-m+1);
return reverse(f.begin(),f.end()),f;
}
polynomial poly_mod(polynomial f,polynomial g){ // return the reminder
int n=f.size()-1,m=g.size()-1;
polynomial q=poly_div(f,g);
int limit=1<<int(ceil(log2(n+1)));
g.resize(limit),q.resize(limit);
num_trans(g,1),num_trans(q,1);
for(int i=0;i<limit;++i) g[i]=mul(g[i],q[i]);
num_trans(g,-1);
for(int i=0;i<m;++i) f[i]=add(f[i],mod-g[i]);
return f.resize(m),f;
}
int n,k;
void mul(polynomial&a,polynomial b,co polynomial&p){
static co int limit=1<<int(ceil(log2(2*k-1)));
a.resize(limit),b.resize(limit);
num_trans(a,1),num_trans(b,1);
for(int i=0;i<limit;++i) a[i]=mul(a[i],b[i]);
num_trans(a,-1),a.resize(2*k-1);
a=poly_mod(a,p);
}
int main(){
read(n),read(k);
polynomial a(k),f(k);
for(int i=1;i<=k;++i) a[k-i]=(mod-read<int>()%mod)%mod; // [-1e9,1e9]
for(int i=0;i<k;++i) f[i]=(read<int>()%mod+mod)%mod;
if(n<k) return printf("%d
",f[n]),0;
a.push_back(1);
polynomial rmd(1,1),tmp(2);tmp[1]=1;
for(;n;n>>=1,mul(tmp,tmp,a))
if(n&1) mul(rmd,tmp,a);
int ans=0;
for(int i=0;i<k;++i) ans=add(ans,mul(rmd[i],f[i]));
printf("%d
",ans);
return 0;
}