【bzoj3625】【CF438E/round250E】小朋友和二叉树【FFT/NTT】【多项式求逆】【多项式开根】

题目链接
题解：我们设f[i]表示权值为i的二叉树的数目，g[i]为i这个值是否在C集合中出现。则很容易得到f[x]=∑xi−1g[i]∑x−ij=0f[j]f[x−i−j]$f[x]=\sum _{i-1}^{x}g[i]\sum _{j=0}^{x-i}f[j]f[x-i-j]$，且f[0]=1$f[0]=1$。
因此，观察可得f=gf2+1$f=g{f}^2+1$
=>gf2−f+1=0$g{f}^2-f+1=0$
=>f=1±1−4g√2g$f=\frac{1\pm \sqrt{1-4g}}{2g}$
=>f=4g2g(1±1−4g√)$f=\frac{4g}{2g(1\pm \sqrt{1-4g})}$
=>f=21±1−4g√$f=\frac{2}{1\pm \sqrt{1-4g}}$
因为g常数项为0，所以1−4g−−−−−√$\sqrt{1-4g}$常数项必定为1。如果f=21±1−4g√$f=\frac{2}{1\pm \sqrt{1-4g}}$取负号，无法进行多项式求逆。但是蒟蒻博主暂时没有搞明白取负号为什么一定不行。
所以f=21+1−4g√$f=\frac{2}{1+\sqrt{1-4g}}$。
因此我们可以通过多项式开根和多项式求逆得到f。
它们的大体思路都是倍增。
多项式求逆：
A(x)B(x)≡1(mod xn)$A(x)B(x)\equiv 1(mod{x}^n)$
A(x)B(x)−1≡0(mod xn)$A(x)B(x)-1\equiv 0(mod{x}^n)$
(A(x)B(x)−1)2≡0(mod x2n)$(A(x)B(x)-1{)}^2\equiv 0(mod{x}^{2n})$
A(x)(2B(x)−B(x)2A(x))≡1(mod x2n)$A(x)(2B(x)-B(x{)}^{2}A(x))\equiv 1(mod{x}^{2n})$
就这样一层一层推上去就可以了。
多项式开根：
B(x)2≡A(x)(mod xn)$B(x{)}^2\equiv A(x)(mod{x}^n)$
B(x)2−A(x)≡0(mod xn)$B(x{)}^2-A(x)\equiv 0(mod{x}^n)$
B(x)4−2B(x)2A(x)+A(x)2≡0(mod x2n)$B(x{)}^4-2B(x{)}^{2}A(x)+A(x{)}^2\equiv 0(mod{x}^{2n})$
B(x)4+2B(x)2A(x)+A(x)2≡4B(x)2A(x)(mod x2n)$B(x{)}^4+2B(x{)}^{2}A(x)+A(x{)}^2\equiv 4B(x{)}^{2}A(x)(mod{x}^{2n})$
(B(x)2+A(x))2≡(2B(x))2A(x)(mod x2n)$(B(x{)}^2+A(x){)}^2\equiv (2B(x){)}^{2}A(x)(mod{x}^{2n})$
(B(x)2+A(x)2B(x))2≡A(x)(mod x2n)$(\frac{B(x{)}^2+A(x)}{2B(x)}{)}^2\equiv A(x)(mod{x}^{2n})$
(B(x)2+A(x)2B(x))≡A(x)(mod x2n)$(\frac{B(x)}{2}+\frac{A(x)}{2B(x)})\equiv A(x)(mod{x}^{2n})$
就这样一层一层推上去就可以了。
最后一步可以减少进行DFT和IDFT的次数，以免被卡常数。
刚开始我自己实现了一个开根，没有化简到最后一步，结果被丧心病狂的B站卡成狗，后来重新学习做人。。。
开根和求逆都是我自己实现的迭代版本，感觉丑到炸。为什么常数却那么大！
具体代码实现：

#include<cstdio>
#include<cmath>
#include<algorithm>
#pragma GCC optimize(3)
using namespace std;
typedef long long ll;
const int N=270005;
const ll mod=998244353;
int n,m,rev[N];
ll a[N],b[N],c[N],d[N],e[N],f[N],g[N];
char ch[N*4+1],*p;
inline int rd(){
    register int res=0;
    while(*p<'0'||*p>'9'){
        p++;
    }
    while(*p>='0'&&*p<='9'){
        res=res*10+*p-'0';
        p++;
    }
    return res;
}
ll fastpow(ll a,ll x){
    a%=mod;
    ll res=1;
    while(x){
        if(x&1){
            res=res*a%mod;
        }
        x>>=1;
        a=a*a%mod;
    }
    return res;
}
ll exgcd(ll a,ll b,ll &x,ll &y){
    if(!b){
        x=1;
        y=0;
        return a;
    }
    ll d=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return d;
}
ll getinv(ll a){
    //return fastpow(a,mod-2);
    ll x,y;
    exgcd(a,mod,x,y);
    return (x+mod)%mod;
}
void ntt(ll *a,int n,int dft){
    for(register int i=0;i<n;i++){
        rev[i]=(rev[i>>1]>>1)|((i&1)*(n>>1));
        if(i<rev[i]){
            swap(a[i],a[rev[i]]);
        }
    }
    for(int i=1;i<n;i<<=1){
        ll wn=fastpow(3,(mod-1)/i/2);
        if(dft==-1){
            wn=getinv(wn);
        }
        for(register int j=0;j<n;j+=i<<1){
            ll w=1,x,y;
            for(register int k=j;k<j+i;k++,w=w*wn%mod){
                x=a[k];
                y=w*a[k+i]%mod;
                a[k]=x+y;
                a[k]=a[k]>=mod?a[k]-mod:a[k]; 
                a[k+i]=x-y;
                a[k+i]=a[k+i]<0?a[k+i]+mod:a[k+i];
            }
        }
    }
    if(dft==-1){
        ll inv=getinv(n);
        for(register int i=0;i<n;i++){
            a[i]=a[i]*inv%mod;
        }
    }
}
void inverse(ll *a,ll *b,ll *c,ll *d,ll n){
    b[0]=getinv(a[0]);
    for(int k=2;k<=n;k<<=1){
        for(register int i=0;i<(k<<1);i++){
            if(i<k){
                c[i]=a[i];
            }else{
                c[i]=0;
            }
            if(i<(k>>1)){
                d[i]=b[i];
            }else{
                d[i]=0;
            }
        }
        ntt(c,k<<1,1);
        ntt(d,k<<1,1);
        for(register int i=0;i<(k<<1);i++){
            c[i]=(2*d[i]%mod-d[i]*d[i]%mod*c[i]%mod+mod)%mod;
        }
        ntt(c,k<<1,-1);
        for(register int i=0;i<k;i++){
            b[i]=c[i];
        }
    }
}
void sqrt(ll *a,ll *b,ll *c,ll *d,ll *e,ll *f){
    b[0]=1;
    for(int k=1;k<=n;k<<=1){
        for(register int i=0;i<k;i++){
            c[i]=2*b[i]%mod;
        }
        inverse(c,d,e,f,k);
        for(register int i=0;i<(k<<1);i++){
            if(i<k){
                c[i]=a[i];
            }else{
                c[i]=d[i]=0;
            }
        }
        ntt(c,k<<1,1);
        ntt(d,k<<1,1);
        for(int i=0;i<(k<<1);i++){
            c[i]=c[i]*d[i]%mod;
        }
        ntt(c,k<<1,-1);
        for(register int i=0;i<(k<<1);i++){
            b[i]=(b[i]*499122177+c[i])%mod;
        }
    }
}
int main(){
    fread(ch,N*4,1,stdin);
    p=ch;
    n=rd(),m=rd();
    for(int i=1;i<=n;i++){
        a[rd()]=1;
    }
    for(n=1;n<=m;n<<=1);
    for(register int i=0;i<n;i++){
        b[i]=a[i]?mod-4:0;
    }
    b[0]=1;
    sqrt(b,c,d,e,f,g);
    c[0]=2;
    inverse(c,b,d,e,n);
    for(register int i=1;i<=m;i++){
        printf("%lld
",b[i]*2%mod);
    }
    return 0;
}