Luogu4723 【模板】常系数齐次线性递推

Luogu4723 【模板】常系数齐次线性递推

参考blog

我们可以求出\(a_n=\sum_{i=0}^{k-1} c_i a_i\)，然后计算答案。

设\(C=\sum_{i=0}^{k-1}c_i x^i\)，\(F=x^k-\sum_{i=0}^{k-1} a_i x^{k-1-i}\)（\(F\)被称为特征多项式）。

根据!@#%^&* 可以得出，\(C=x^n \mod F\)，然后根据\(C\)的系数直接计算即可。

注意到\(n\)很大，我们可以使用倍增法，边倍增边对\(F\)取模。时间复杂度\(O(k \log n \log k)\)。

\(Code:\)

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#define N 300005
#define ll long long
using namespace std;
const int p=998244353;
int n,k,t,clc,ans[N],x[N],a[N],c[N],d[N],F[N],o[N],e[N],h[N],u[N],v[N];
int s,l,E[2][35],rev[N];
void Add(int &x,int y)
{
    x=(x+y>=p)?(x+y-p):(x+y);
}
void Del(int &x,int y)
{
    x=(x>=y)?(x-y):(x-y+p);
}
void Mul(int &x,int y)
{
    x=(ll)x*y%p;
}
int add(int x,int y)
{
    return (x+y>=p)?(x+y-p):(x+y);
}
int del(int x,int y)
{
    return (x>=y)?(x-y):(x-y+p);
}
int mul(int x,int y)
{
    return (ll)x*y%p;
}
int ksm(int x,int y)
{
    int ans(1);
    while (y)
    {
        if (y & 1)
            Mul(ans,x);
        Mul(x,x);
        y >>=1;
    }
    return ans;
}
void Pre()
{
    E[0][23]=ksm(3,(p-1)/(1 << 23));
    E[1][23]=ksm(E[0][23],p-2);
    for (int i=22;i;--i)
    {
        E[0][i]=mul(E[0][i+1],E[0][i+1]);
        E[1][i]=mul(E[1][i+1],E[1][i+1]);
    }
}
void solve(int n)
{
    s=1,l=0;
    while (s<n)
        s <<=1,++l;
    for (int i=0;i<s;++i)
        rev[i]=(rev[i >> 1] >> 1) | ((i & 1) << l-1);
}
void Cut(int *a,int n,int s)
{
    if (n>s)
        return;
    memset(a+n,0,(s-n)*sizeof(int));
}
void NTT(int *a,int t)
{
    for (int i=0;i<s;++i)
        if (i<rev[i])
            swap(a[i],a[rev[i]]);
    for (int mid=1,o=1;mid<s;mid <<=1,++o)
        for (int j=0;j<s;j+=(mid << 1))
        {
            int g(1);
            for (int k=0;k<mid;++k,Mul(g,E[t][o]))
            {
                int x(a[j+k]),y(mul(g,a[j+k+mid]));
                a[j+k]=add(x,y),a[j+k+mid]=del(x,y);
            }
        }
}
void GetInv(int *f,int *g,int R)
{
    if (R==2)
    {
        g[0]=ksm(f[0],p-2);
        return;
    }
    GetInv(f,g,R >> 1);
    memcpy(c,g,(R >> 2)*sizeof(int)),memcpy(d,f,(R >> 1)*sizeof(int));
    solve(R),NTT(c,0),NTT(d,0);
    for (int i=0;i<s;++i)
        c[i]=del(add(c[i],c[i]),mul(d[i],mul(c[i],c[i])));
    NTT(c,1);
    int iv(ksm(s,p-2));
    for (int i=0;i<s;++i)
        Mul(c[i],iv);
    memcpy(g,c,(R >> 1)*sizeof(int)),memset(c,0,s*sizeof(int)),memset(d,0,s*sizeof(int));
}
void PolyMod(int *f,int *g,int *R,int n,int m)
{
    for (int i=0;i<n;++i)
        e[n-1-i]=f[i];
    for (int i=0;i<m;++i)
        o[m-1-i]=g[i];
    Cut(e,n-m+1,n);
    Cut(o,n-m+1,m);
    s=1;
    while (s<n-m+1)
        s <<=1;
    s <<=1;
    GetInv(o,h,s);
    memset(o,0,s*sizeof(int));
    Cut(h,n-m+1,s);
    solve((n-m+1) << 1);
    NTT(e,0),NTT(h,0);
    for (int i=0;i<s;++i)
        Mul(h[i],e[i]);
    NTT(h,1);
    int iv(ksm(s,p-2));
    for (int i=0;i<s;++i)
        Mul(h[i],iv);
    memset(e,0,s*sizeof(int));
    Cut(h,n-m+1,s);
    reverse(h,h+n-m+1);
    solve(n-m+1+m);
    NTT(h,0),NTT(g,0);
    for (int i=0;i<s;++i)
        Mul(h[i],g[i]);
    NTT(h,1),NTT(g,1);
    iv=ksm(s,p-2);
    for (int i=0;i<s;++i)
        Mul(g[i],iv),Mul(h[i],iv);
    for (int i=0;i<m-1;++i)
        R[i]=del(f[i],h[i]);
    memset(h,0,s*sizeof(int)),memset(f,0,s*sizeof(int));
}
void PolyMul(int *a,int *b)
{
    memcpy(u,a,(k+1)*sizeof(int)),memcpy(v,b,(k+1)*sizeof(int));
    solve(k+1 << 1);
    NTT(u,0),NTT(v,0);
    for (int i=0;i<s;++i)
        Mul(u[i],v[i]);
    NTT(u,1);
    int iv(ksm(s,p-2));
    for (int i=0;i<s;++i)
        Mul(u[i],iv);
    memset(v,0,s*sizeof(int));
    memset(a,0,(k+2)*sizeof(int));
    PolyMod(u,F,a,(k+1 << 1),k+1);
}
int main()
{
    Pre();
    scanf("%d%d",&n,&k);
    for (int i=0;i<k;++i)
        scanf("%d",&a[i]),a[i]=(a[i]%p+p)%p;
    for (int i=0;i<k;++i)
        F[k-1-i]=(p-a[i])%p;
    F[k]=1;
    ans[0]=1,x[1]=1;
    while (n)
    {
        if (n & 1)
            PolyMul(ans,x);
        PolyMul(x,x);
        n >>=1;
    }
    for (int i=0;i<k;++i)
    {
        scanf("%d",&t),Add(t,p);
        Add(clc,mul(t,ans[i]));
    }
    Add(clc,p);
    printf("%d\n",clc);
    return 0;
}