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

题目链接
题解:我们设f[i]表示权值为i的二叉树的数目,g[i]为i这个值是否在C集合中出现。则很容易得到f[x]=i1xg[i]j=0xif[j]f[xij],且f[0]=1
因此,观察可得f=gf2+1
=>gf2f+1=0
=>f=1±14g2g
=>f=4g2g(1±14g)
=>f=21±14g
因为g常数项为0,所以14g常数项必定为1。如果f=21±14g取负号,无法进行多项式求逆。但是蒟蒻博主暂时没有搞明白取负号为什么一定不行。
所以f=21+14g
因此我们可以通过多项式开根和多项式求逆得到f。
它们的大体思路都是倍增。
多项式求逆:
A(x)B(x)1(mod xn)
A(x)B(x)10(mod xn)
(A(x)B(x)1)20(mod x2n)
A(x)(2B(x)B(x)2A(x))1(mod x2n)
就这样一层一层推上去就可以了。
多项式开根:
B(x)2A(x)(mod xn)
B(x)2A(x)0(mod xn)
B(x)42B(x)2A(x)+A(x)20(mod x2n)
B(x)4+2B(x)2A(x)+A(x)24B(x)2A(x)(mod x2n)
(B(x)2+A(x))2(2B(x))2A(x)(mod x2n)
(B(x)2+A(x)2B(x))2A(x)(mod x2n)
(B(x)2+A(x)2B(x))A(x)(mod x2n)
就这样一层一层推上去就可以了。
最后一步可以减少进行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\n",b[i]*2%mod);
    }
    return 0;
}
posted @ 2018-05-13 18:02  ez_2016gdgzoi471  阅读(...)  评论(...编辑  收藏