P10007 Solution（套路合集？）

· · 题解

一看这题，我啪的一下就点进来了，很快啊！

设 F(x) 是题面代码中序列 f 的生成函数，那么发现每次操作其实就是

F(x) \leftarrow \frac{F(x)}{1-ax}

也就是说答案（设 s_t=\text{ans}_t）

s_t=[x^t]\frac{F(x)}{(1-ax)^t}

这是经典问题，可以做到 \Theta(n \log n) 的时间复杂度。

啊好像还需要一些常数优化，那就展开解释一下吧。首先套用模板，设 h(x) 为 x-ax^2 的复合逆，即

h(x)=\frac{1-\sqrt{1-4ax}}{2a}

答案就是

s_t=[x^0] \frac{F(x)}{(x-ax^2)^t}=[x^t]F(h(x)) \frac{xh'(x)}{h(x)}

虽然 xh'(x)/h(x) 微分有限，但是在最后还需要一次卷积。反正 F(x) 是给定的，没什么特殊性质，不如考虑把这部分揉进 F(h(x)) 里面去。根据代数方程的性质可知

h'(x)=\frac{1-2a h(x)}{1-4ax} \ , \ \frac{x}{h(x)}=1-ah(x)

然后设 \hat F(x)=F(x)(1-2ax)(1-ax)，这样处理起来就比较方便了：

s_t=[x^t]\frac{\hat F(h(x))}{1-4ax}

实际上，接下来的处理也很经典：

\hat F(h(x))=\hat F\left(\frac{1-x}{2a}\right) \circ \sqrt{1-4ax}

设 G(x)=\hat F((1-x)/2a)，然后将其分解为 G(x)=P(x^2)+xQ(x^2)，就有 G(\sqrt{1-4ax})=P(1-4ax)+\sqrt{1-4ax}Q(1-4ax)。

最后的做法似乎比出题人的常数稍大一点，但似乎容易处理一些。

update：好像实现效果比出题人还快了近一倍，很神秘。我的 DFT 部分也没有做精细的优化，补一下代码吧。

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#define N 2097158
#define ll long long
#define p 998244353
using namespace std;

inline void read(int &x){
    x = 0;
    char c = getchar();
    while(c<'0'||c>'9') c = getchar();
    while(c>='0'&&c<='9'){
        x = (x<<3)+(x<<1)+(c^48);
        c = getchar();
    }
}

inline int power(int a,int t){
    int res = 1;
    while(t){
        if(t&1) res = (ll)res*a%p;
        a = (ll)a*a%p;
        t >>= 1;
    }
    return res;
}

int rt[N],rev[N],fac[N],ifac[N];
int siz;

void init(int n){
    int w,lim = 1;
    while(lim<=n) lim <<= 1,++siz;
    for(int i=1;i!=lim;++i) rev[i] = (rev[i>>1]>>1)|((i&1)<<(siz-1));
    w = power(3,(p-1)>>siz);
    fac[0] = fac[1] = ifac[0] = ifac[1] = rt[lim>>1] = 1;
    for(int i=lim>>1|1;i!=lim;++i) rt[i] = (ll)rt[i-1]*w%p;
    for(int i=(lim>>1)-1;i;--i) rt[i] = rt[i<<1];
    for(int i=2;i<=n;++i) fac[i] = (ll)fac[i-1]*i%p;
    ifac[n] = power(fac[n],p-2);
    for(int i=n-1;i;--i) ifac[i] = (ll)ifac[i+1]*(i+1)%p;
}

inline void dft(int *f,int lim){
    static unsigned long long a[N];
    int x,shift = siz-__builtin_ctz(lim);
    for(int i=0;i!=lim;++i) a[rev[i]>>shift] = f[i];
    for(int mid=1;mid!=lim;mid<<=1)
    for(int j=0;j!=lim;j+=(mid<<1))
    for(int k=0;k!=mid;++k){
        x = a[j|k|mid]*rt[mid|k]%p;
        a[j|k|mid] = a[j|k]+p-x;
        a[j|k] += x;
    }
    for(int i=0;i!=lim;++i) f[i] = a[i]%p;
}

inline void idft(int *f,int lim){
    reverse(f+1,f+lim);
    dft(f,lim);
    int x = p-((p-1)>>__builtin_ctz(lim));
    for(int i=0;i!=lim;++i) f[i] = (ll)f[i]*x%p;
}

inline int getlen(int n){ return 1<<(32-__builtin_clz(n)); }

inline void composite(const int *f,int n,int c,int *R){
    static int g[N],h[N];
    g[0] = 1,h[0] = f[0];
    for(int i=1;i<=n;++i){
        g[i] = (ll)g[i-1]*c%p;
        h[i] = (ll)f[i]*fac[i]%p;
    }
    for(int i=2;i<=n;++i) g[i] = (ll)g[i]*ifac[i]%p;
    reverse(g,g+n+1);
    int lim = getlen(n<<1);
    memset(g+n+1,0,(lim-n)<<2);
    memset(h+n+1,0,(lim-n)<<2);
    dft(g,lim),dft(h,lim);
    for(int i=0;i!=lim;++i) g[i] = (ll)g[i]*h[i]%p;
    idft(g,lim);
    for(int i=0;i<=n;++i) R[i] = (ll)g[i+n]*ifac[i]%p;
}

inline int binom(int n,int m){
    if(n<m) return 0;
    return (ll)fac[n]*ifac[m]%p*ifac[n-m]%p;
}

int n,a,dp,dq;
int f[N],P[N],Q[N],h[N];
const int inv2 = (p+1)/2;

int main(){
    read(n),read(a);
    init(n<<1);
    for(int i=0;i<=n;++i) read(f[i]);
    for(int i=n+1;i;--i) f[i] = (f[i]-(ll)a*f[i-1])%p;
    for(int i=n+2;i;--i) f[i] = (f[i]-2ll*a*f[i-1])%p;
    int pw = 1,c = power(p-2*a%p,p-2);
    n += 2;
    for(int i=0;i<=n;++i){
        f[i] = ((ll)(f[i]+p)*pw)%p;
        pw = (ll)pw*c%p;
    }
    composite(f,n,p-1,f);
    for(int i=0;i<=n;++i){
        if(i&1) Q[dq++] = f[i];
        else P[dp++] = f[i];
    }
    composite(P,dp-1,1,P),composite(Q,dq-1,1,Q);
    pw = 1,c = p-4ll*a%p;
    for(int i=0;(i<<1)<=n+1;++i){
        P[i] = (ll)P[i]*pw%p,Q[i] = (ll)Q[i]*pw%p;
        pw = (ll)pw*c%p;
    }
    h[0] = pw = 1,c = 4ll*(p-a)%p;
    for(int i=1;i<=n;++i){
        pw = (ll)pw*c%p*(inv2-i+1)%p;
        h[i] = (ll)pw*ifac[i]%p;
    }
    int lim = getlen(n<<1);
    dft(h,lim),dft(Q,lim);
    for(int i=0;i!=lim;++i) Q[i] = (ll)Q[i]*h[i]%p;
    idft(Q,lim);
    memset(Q+n+1,0,(lim-n)<<2);
    n -= 2;
    for(int i=0;i<=n;++i) f[i] = P[i]+Q[i]>=p?P[i]+Q[i]-p:P[i]+Q[i];
    for(int i=1;i<=n;++i){
        f[i] = (f[i]+4ll*a*f[i-1])%p;
        printf("%d ",f[i]);
    }
    return 0;   
}