常系数齐次线性递推快速算法学习笔记

最新推荐文章于 2023-01-29 13:54:33 发布

angzuo8655

最新推荐文章于 2023-01-29 13:54:33 发布

阅读量273

点赞数

原文链接：http://www.cnblogs.com/1000Suns/p/11266855.html

版权

今天集训被线代狠狠的虐了一发。

不过还有一点收获的，比如这个。

数列 \(f\) 满足 \(f_n=\sum\limits_{i=1}^ka_if_{n-i}(n\ge k)\)，其中 \(a_1\dots a_k,f_0\dots f_{k-1}\) 均给出。求 \(f_n\)。

\(n\le 10^9,k\le 30000\)。

先要弄懂一些基本定义，

矩阵，行列式，高斯消元这些基本的东西就自己看别的东西去吧，我也不知道有什么好的资料，不妨去洛谷搜模板看题解，那里面的都不错。

然后讲一下特征值，特征多项式和 Hamilton-Cayley 定理，就能做这题了。

对于 \(n\times n\) 的矩阵 \(A\)，如果存在数 \(\lambda\) 和非零列向量 \(x\) 满足 \(Ax=\lambda x\)，那么 \(\lambda\) 是 \(A\) 的特征值，\(x\) 是 \(A\) 的特征向量。

那么 \(Ax=\lambda I x\)。\((\lambda I-A)x=0\)。因为 \(x\) 不为零向量，所以 \(\det(\lambda I-A)=0\)，也就是 \(\lambda I-A\) 不满秩。

我们称 \(\det(\lambda I-A)=0\) 为 \(A\) 的特征多项式，元是 \(\lambda\)。特征多项式是 \(n\) 次的，他的 \(n\) 个根就是 \(A\) 的所有特征值。（可能有相等的根）

对于上三角矩阵，所有的特征值就是主对角线上的所有值。

如果有 \(n\) 个线性无关的特征向量 \(x_i\)（当且仅当 \(A\) 满秩），那么有 \(A\begin{bmatrix}x_1&x_2&\cdots&x_n\end{bmatrix}=\begin{bmatrix}x_1&x_2&\cdots&x_n\end{bmatrix}\begin{bmatrix}\lambda_1&0&0&\cdots&0\\0&\lambda_2&0&\cdots&0\\\cdots\\0&0&0&\cdots&\lambda_n\end{bmatrix}\)。

Hamilton-Cayley 定理：对于矩阵 \(A\) 的特征多项式 \(f(\lambda)=\sum\limits_{i=0}^nc_i\lambda^i\)，有 \(f(A)=0\)，即 \(\sum\limits_{i=0}^nc_iA^i=0\)。

会了这些，就可以开始了。

\(O(k^3\log n)\) 的相信大家都会。（什么？不会？赶快去学矩阵快速幂）

首先考虑写出转移矩阵 \(A\) 和初始行向量 \(f\)，我们要求的是 \(f\times A^n\) 的第 \(0\) 维。

假如我们构造出了一个序列 \(c\) 使得

\[A^n=\sum\limits_{i=0}^{k-1}c_iA^i\]

那么有：

\[f\times A^n=\sum\limits_{i=0}^{k-1}c_i(f\times A^i)\]

\[(f\times A^n)_0=\sum\limits_{i=0}^{k-1}c_i(f\times A^i)_0\]

\[f_n=\sum\limits_{i=0}^{k-1}c_if_i\]

那么就可以 \(O(k)\) 计算了。

那么 \(c\) 怎么弄呢？

令 \(R(A)=\sum\limits_{i=0}^{k-1}c_iA^i\)

假如存在 \(k\) 次多项式 \(G(A)\)，使得 \(A^n=F(A)G(A)+R(A)\)。(标准多项式除法形式)

当 \(G(A)=0\) 时就有 \(A^n=R(A)\)，所以要求的就是 \(A^n\bmod G(A)\)，快速幂+多项式除法 \(O(k\log k\log n)\) 解决。

那么如何构造 \(G(A)=0\) 的多项式呢？

看到上面的 Hamilton-Cayley 定理，令 \(G(\lambda)=\det(\lambda I-A)\) 即可。

手玩一下，发现 \(\det(\lambda I-A)=-\sum\limits_{i=0}^{k-1}\lambda^ia_{k-i}+\lambda^k\)。（注意交换行的时候取相反数！！！）

所以 \(g_i=-a_{k-i}(i\ne k),g_k=1\)。

时间复杂度 \(O(k\log k\log n)\)。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=333333,mod=998244353;
#define FOR(i,a,b) for(int i=(a);i<=(b);i++)
#define ROF(i,a,b) for(int i=(a);i>=(b);i--)
#define MEM(x,v) memset(x,v,sizeof(x))
inline int read(){
    int x=0,f=0;char ch=getchar();
    while(ch<'0' || ch>'9') f|=ch=='-',ch=getchar();
    while(ch>='0' && ch<='9') x=x*10+ch-'0',ch=getchar();
    return f?-x:x;
}
int n,k,s,f[maxn],g[maxn],fac[maxn],ans[maxn],prod[maxn],tmp[maxn],lim,l,rev[maxn],invt[maxn],Arev[maxn],Btmp[maxn],Brev[maxn],Brevinv[maxn],C[maxn];
inline int add(int x,int y){return x+y<mod?x+y:x+y-mod;}
inline int sub(int x,int y){return x<y?x-y+mod:x-y;}
inline int mul(int x,int y){return 1ll*x*y%mod;}
inline int qpow(int a,int b){
    int ans=1;
    for(;b;b>>=1,a=mul(a,a)) if(b&1) ans=mul(ans,a);
    return ans;
}
void init(int upr){
    for(lim=1,l=0;lim<upr;lim<<=1,l++);
    FOR(i,0,lim-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
}
void NTT(int *A,int tp){
    FOR(i,0,lim-1) if(i<rev[i]) swap(A[i],A[rev[i]]);
    for(int i=1;i<lim;i<<=1)
        for(int j=0,r=i<<1,Wn=qpow(3,mod-1+tp*(mod-1)/r);j<lim;j+=r)
            for(int k=0,w=1;k<i;k++,w=mul(w,Wn)){
                int x=A[j+k],y=mul(w,A[i+j+k]);
                A[j+k]=add(x,y);
                A[i+j+k]=sub(x,y);
            }
    if(tp==-1){
        int linv=qpow(lim,mod-2);
        FOR(i,0,lim-1) A[i]=mul(A[i],linv);
    }
}
void inv(int *A,int *B,int deg){
    if(deg==1) return void(B[0]=qpow(A[0],mod-2));
    inv(A,B,(deg+1)>>1);
    init(deg<<1);
    FOR(i,0,deg-1) invt[i]=A[i];
    FOR(i,deg,lim-1) invt[i]=0;
    NTT(invt,1);NTT(B,1);
    FOR(i,0,lim-1) B[i]=mul(B[i],sub(2,mul(invt[i],B[i])));
    NTT(B,-1);
    FOR(i,deg,lim-1) B[i]=0;
}
void division(int *A,int *B,int *D,int n,int m){
    if(n<m){
        init(n<<1);
        FOR(i,n+1,lim-1) D[i]=A[i]=0;
        FOR(i,0,n) D[i]=A[i];
        return;
    }
    init(n<<1);
    FOR(i,0,lim-1) Arev[i]=Brev[i]=Brevinv[i]=C[i]=Btmp[i]=D[i]=0;
    FOR(i,n+1,lim-1) A[i]=0;
    FOR(i,m+1,lim-1) B[i]=0;
    FOR(i,0,n) Arev[i]=A[n-i];
    FOR(i,0,n-m) Brev[i]=B[m-i];
    inv(Brev,Brevinv,n-m+1);
    init(n<<1);
    NTT(Arev,1);NTT(Brevinv,1);
    FOR(i,0,lim-1) C[i]=mul(Arev[i],Brevinv[i]);
    NTT(C,-1);
    FOR(i,0,(n-m)/2) swap(C[i],C[n-m-i]);
    FOR(i,n-m+1,lim-1) C[i]=0;
    FOR(i,0,m) Btmp[i]=B[i];
    NTT(Btmp,1);NTT(C,1);
    FOR(i,0,lim-1) D[i]=mul(Btmp[i],C[i]);
    NTT(D,-1);
    FOR(i,0,m-1) D[i]=sub(A[i],D[i]);
    FOR(i,m,lim-1) D[i]=0;
}
int main(){
    n=read();k=read();
    FOR(i,1,k) g[k-i]=(mod-read())%mod;
    FOR(i,0,k-1) f[i]=(read()+mod)%mod;
    g[k]=fac[1]=ans[0]=1;
    while(n){
        if(n&1){
            init(k<<1);
            NTT(fac,1);NTT(ans,1);
            FOR(i,0,lim-1) prod[i]=mul(fac[i],ans[i]);
            NTT(fac,-1);NTT(ans,-1);NTT(prod,-1);
            division(prod,g,ans,(k-1)<<1,k);
        }
        init(k<<1);
        NTT(fac,1);
        FOR(i,0,lim-1) prod[i]=mul(fac[i],fac[i]);
        NTT(fac,-1);NTT(prod,-1);
        division(prod,g,fac,(k-1)<<1,k);
        n>>=1;
    }
    FOR(i,0,k-1) s=add(s,mul(f[i],ans[i]));
    printf("%d\n",s);
}

转载于:https://www.cnblogs.com/1000Suns/p/11266855.html