常系数齐次线性递推

现在终于把线性递推常数减小了许多。（实际上我原先写的多项式取模才这么慢）

Fiduccia 算法

一个方法是求 \(x^n\) 对递推数列的特征多项式

\[p(x)=x^k-p_1x^{k-1}-p_2x^{k-2}-\cdots-p_k \]

取模的结果，然后逐项代入即可求得答案。这个取模可以像普通的快速幂一样，每次对 \(p(x)\) 取模。证明不会。

但是这玩意跑的非常慢，卡了半天常仍然只能做到每个点 1.7s 的水平。当然我本身大常数选手。

这玩意当然不是重点。放个代码仅供参考。

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <iostream>
#include <cmath>
using namespace std;
const int mod=998244353,g=3;
int n,k,wl,ans,a[150010],b[150010],r[150010],inv[150010],res[150010];
int q[150010],f[150010],c[150010],p[150010],tmp[150010];
int val[150010];
void get(int n){
    wl=1;
    while(n>=wl)wl<<=1;
    for(int i=0;i<=wl;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(__lg(wl)-1));
}
int qpow(int a,int b){
    int ans=1;
    while(b){
        if(b&1)ans=1ll*a*ans%mod;
        a=1ll*a*a%mod;
        b>>=1;
    }
    return ans;
}
const int invg=qpow(g,mod-2);
void ntt(int a[],int n,int tp){
    for(int i=1;i<n;i++)if(i<r[i])swap(a[i],a[r[i]]);
    for(int mid=1;mid<n;mid<<=1){
        int wn=qpow(tp==1?g:invg,(mod-1)/(mid<<1));
        for(int j=0;j<n;j+=mid<<1){
            int w=1;
            for(int k=0;k<mid;k++,w=1ll*w*wn%mod){
                int x=a[j+k],y=1ll*w*a[j+mid+k]%mod;
                a[j+k]=(x+y)%mod;a[j+mid+k]=(x-y+mod)%mod;
            }
        }
    }
    if(tp^1){
        int inv=qpow(n,mod-2);
        for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
    }
}
void getinv(int n,int a[],int b[]){
    if(n==1){
        b[0]=qpow(a[0],mod-2);
        return;
    }
    getinv((n+1)>>1,a,b);
    get(n<<1);
    for(int i=0;i<n;i++)c[i]=a[i];
    for(int i=n;i<wl;i++)c[i]=0;
    ntt(b,wl,1);ntt(c,wl,1);
    for(int i=0;i<wl;i++)b[i]=1ll*(2-1ll*b[i]*c[i]%mod+mod)%mod*b[i]%mod;
    ntt(b,wl,-1);
    for(int i=n;i<wl;i++)b[i]=0;
}
void getdiv(int n,int m,int a[],int b[]){
    get(n<<1);
    for(int i=0;i<wl;i++)q[i]=res[i]=c[i]=0;
    for(int i=0;i<=n;i++)q[n-i]=a[i];
    ntt(q,wl,1);
    for(int i=0;i<wl;i++)q[i]=1ll*q[i]*inv[i]%mod;
    ntt(q,wl,-1);
    for(int i=n-m+1;i<wl;i++)q[i]=0;
    reverse(q,q+n-m+1);
    ntt(b,wl,1);ntt(q,wl,1);
    for(int i=0;i<wl;i++)b[i]=1ll*b[i]*q[i]%mod;
    ntt(b,wl,-1);
    for(int i=0;i<m;i++)res[i]=(a[i]-b[i]+mod)%mod;
}
int main(){
    scanf("%d%d",&n,&k);
    for(int i=1;i<=k;i++)scanf("%d",&p[i]),p[i]=(p[i]+mod)%mod;
    reverse(p,p+k+1);p[k]=1;
    for(int i=0;i<k;i++)p[i]=(mod-p[i])%mod;
    for(int i=0;i<=k;i++)tmp[i]=f[k-i]=p[i];
    getinv(k-1,f,inv);
    get(k<<2);ntt(inv,wl,1);
    for(int i=0;i<k;i++)scanf("%d",&val[i]),val[i]=(val[i]+mod)%mod;
    b[1]=a[0]=1;
    while(n){
        if(n&1){
            get(k<<1);
            ntt(a,wl,1);ntt(b,wl,1);
            for(int i=0;i<wl;i++)a[i]=1ll*a[i]*b[i]%mod;
            ntt(a,wl,-1);ntt(b,wl,-1);
            getdiv((k-1)<<1,k,a,p);
            memcpy(a,res,sizeof(a));
            memcpy(p,tmp,sizeof(p));
        }
        get(k<<1);
        ntt(b,wl,1);
        for(int i=0;i<wl;i++)b[i]=1ll*b[i]*b[i]%mod;
        ntt(b,wl,-1);
        getdiv((k-1)<<1,k,b,p);
        memcpy(b,res,sizeof(b));
        memcpy(p,tmp,sizeof(p));
        n>>=1;
    }
    for(int i=0;i<k;i++)ans=(ans+1ll*val[i]*a[i])%mod;
    printf("%d\n",ans);
    return 0;
}

LSB-first 算法

重点来了。首先这玩意只能在模数是质数的情况下使用。

它的思想在于把递推数列转为生成函数的形式：

\[F(x)=\frac{P(x)}{Q(x)} \]

其中 \(Q(x)\) 是递推式的生成函数（就是特征多项式的系数翻转），而 \(P(x)\) 是一个次数小于 \(k\) 的多项式。

一个结论是一定存在这样的一个 \(P(x)\)。证明还是不会。那么我们就可以用 \(F(x)\) 的前 \(k\) 项和 \(Q(x)\) 计算 \(P(x)\) 的前 \(k\) 项。

然后考虑如何提取 \([x^n]\frac{P(x)}{Q(x)}\)。我们使用 Bostan-Mori 算法解决这个问题。

该算法依赖于

\[\frac{P(x)}{Q(x)}=\frac{P(x)Q(-x)}{Q(x)Q(-x)} \]

这玩意分母显然只有偶数次项不是 \(0\) ，那么求逆之后也是只有偶数次项不是 \(0\)，也就是分子的奇数次项和偶数次项互不影响。那假如说我们要提取的 \(x^n\) 的 \(n\) 是奇数的话，那么分子的偶数次项就可以全清掉。反之同理。这样每次做一次多项式乘法然后对分子分母的奇/偶次项重新标号，就可以把 \(n\) 减半。同时由于多项式乘法使得次数变成 \(2k\)，这样重标号之后次数就仍然是 \(k\)。复杂度和上边的一样都是 \(O(k\log k\log n)\)。

这玩意非常好写而且跑的很快，我写的这份时间空间都是上边那个的四分之一，450ms。

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <iostream>
#include <cmath>
using namespace std;
const int mod=998244353,g=3;
int n,k,wl,ans,f[150010],r[150010];
int q[150010],p[150010];
void get(int n){
    wl=1;
    while(n>=wl)wl<<=1;
    for(int i=0;i<=wl;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(__lg(wl)-1));
}
int qpow(int a,int b){
    int ans=1;
    while(b){
        if(b&1)ans=1ll*a*ans%mod;
        a=1ll*a*a%mod;
        b>>=1;
    }
    return ans;
}
const int invg=qpow(g,mod-2);
void ntt(int a[],int n,int tp){
    for(int i=1;i<n;i++)if(i<r[i])swap(a[i],a[r[i]]);
    for(int mid=1;mid<n;mid<<=1){
        int wn=qpow(tp==1?g:invg,(mod-1)/(mid<<1));
        for(int j=0;j<n;j+=mid<<1){
            int w=1;
            for(int k=0;k<mid;k++,w=1ll*w*wn%mod){
                int x=a[j+k],y=1ll*w*a[j+mid+k]%mod;
                a[j+k]=(x+y)%mod;a[j+mid+k]=(x-y+mod)%mod;
            }
        }
    }
    if(tp^1){
        int inv=qpow(n,mod-2);
        for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
    }
}
int main(){
    scanf("%d%d",&n,&k);q[0]=1;
    for(int i=1;i<=k;i++){
        int x;scanf("%d",&x);x=(x+mod)%mod;
        q[i]=(mod-x)%mod;
    }
    for(int i=0;i<k;i++){
        int x;scanf("%d",&x);p[i]=(x+mod)%mod;
    }
    for(int i=0;i<=k;i++)f[i]=q[i];
    get(k<<1);
    ntt(p,wl,1);ntt(f,wl,1);
    for(int i=0;i<wl;i++)p[i]=1ll*p[i]*f[i]%mod;
    ntt(p,wl,-1);
    for(int i=k;i<wl;i++)p[i]=0;
    while(n){
        for(int i=0;i<=k;i++)f[i]=(i&1)?(mod-q[i])%mod:q[i];
        for(int i=k+1;i<wl;i++)f[i]=0;
        ntt(f,wl,1);ntt(p,wl,1);ntt(q,wl,1);
        for(int i=0;i<wl;i++)p[i]=1ll*p[i]*f[i]%mod,q[i]=1ll*q[i]*f[i]%mod;
        ntt(p,wl,-1);ntt(q,wl,-1);
        for(int i=0;i<=k;i++)p[i]=p[(i<<1)|(n&1)],q[i]=q[i<<1];
        for(int i=k+1;i<wl;i++)p[i]=q[i]=0;
        n>>=1;
    }
    ans=1ll*p[0]*qpow(q[0],mod-2)%mod;
    printf("%d\n",ans);
    return 0;
}