[LOJ#3120][Luogu5401][CTS2019]珍珠(容斥+生成函数)

https://www.luogu.org/blog/user50971/solution-p5401

#include<cstdio>
#include<algorithm>
#define rep(i,l,r) for (int i=(l); i<=(r); i++)
using namespace std;

const int N=400010,mod=998244353,i2=499122177;
int D,n,m,ans,fac[N],inv[N],ip2[N],f[N],g[N],rev[N],a[N],b[N];

int ksm(int a,int b){
    int res=1;
    for (; b; a=1ll*a*a%mod,b>>=1)
        if (b & 1) res=1ll*res*a%mod;
    return res;
}

int C(int n,int m){ return n<m ? 0 : 1ll*fac[n]*inv[m]%mod*inv[n-m]%mod; }

void NTT(int a[],int n,bool f){
    for (int i=0; i<n; i++) if (i<rev[i]) swap(a[i],a[rev[i]]);
    for (int i=1; i<n; i<<=1){
        int wn=ksm(3,f?(mod-1)/(i<<1):(mod-1)-(mod-1)/(i<<1));
        for (int p=i<<1,j=0; j<n; j+=p){
            int w=1;
            for (int k=0; k<i; k++,w=1ll*w*wn%mod){
                int x=a[j+k],y=1ll*w*a[i+j+k]%mod;
                a[j+k]=(x+y)%mod; a[i+j+k]=(x-y+mod)%mod;
            }
        }
    }
    if (f) return;
    int inv=ksm(n,mod-2);
    for (int i=0; i<n; i++) a[i]=1ll*a[i]*inv%mod;
}

int main(){
    freopen("pearl.in","r",stdin);
    freopen("pearl.out","w",stdout);
    scanf("%d%d%d",&D,&n,&m);
    if (n-2*m<0){ puts("0"); return 0; }
    if (n-2*m>=D){ printf("%d
",ksm(D,n)); return 0; }
    fac[0]=ip2[0]=1;
    rep(i,1,D) fac[i]=1ll*fac[i-1]*i%mod,ip2[i]=1ll*ip2[i-1]*i2%mod;
    inv[D]=ksm(fac[D],mod-2);
    for (int i=D; i; i--) inv[i-1]=1ll*inv[i]*i%mod;
    
    int l=1,L=0,t=1;
    while (l<=(D<<1)) l<<=1,L++;
    for (int i=0; i<l; i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    rep(i,0,D) a[i]=1ll*t*ksm(D-2*i+mod+mod,n)%mod*inv[i]%mod,b[i]=inv[i],t=mod-t;
    NTT(a,l,1); NTT(b,l,1);
    for (int i=0; i<l; i++) a[i]=1ll*a[i]*b[i]%mod;
    NTT(a,l,0);
    rep(i,0,D) f[i]=1ll*ip2[i]*fac[i]%mod*C(D,i)%mod*a[i]%mod;
    for (int i=0; i<l; i++) a[i]=b[i]=0;
    
    t=(D&1)?mod-1:1;
    rep(i,0,D) a[i]=1ll*f[i]*fac[i]%mod,b[i]=1ll*t*inv[D-i]%mod,t=mod-t;
    NTT(a,l,1); NTT(b,l,1);
    for (int i=0; i<l; i++) a[i]=1ll*a[i]*b[i]%mod;
    NTT(a,l,0);
    rep(i,0,D) g[i]=1ll*a[D+i]*inv[i]%mod;
    rep(i,0,n-2*m) ans=(ans+g[i])%mod;
    printf("%d
",ans);
    return 0;
}