LOJ#3120. 珍珠 容斥+生成函数+NTT

求：有多少种序列满足 $a[i] subseteq [1,D]$ 且 $m leqslant sum_{i=1}^{D} frac{cnt[i]}{2}$
裸做的话就是一个背包：$f[i][j]$ 表示前 $i$ 个数匹配了 $j$ 对的方案数，然后由于没有匹配的肯定是单个出现，所以转移的话比较容易.
再加上几个点的特判，和不知道能不能过的矩阵乘法，大概能拿 72 分.
然后分析上面那个 $sum_{i=1}^{D} frac{cnt[i]}{2}$，由于这个是向下取整，所以一个数出现次为奇数的话就会少一对贡献.
那么出现次数为奇数的个数就不超过 $n-2 imes m.$
不妨枚举出现次数为奇数的有多少个：
$f_{i}=inom{D}{i}n! (frac{e^x-e^{-x}}{2})^i (frac{e^x+e^{-x}}{2}) ^{D-i}[x^n]$
但是这个式子不好推，因为前后有两个 $frac{e^x + e^{-x}}{2}$ 这样的生成函数.
不妨考虑二项式反演，即钦定一些位置是奇数，然后其余位置随便填.
$Rightarrow f_{i}=inom{D}{i}n! (frac{e^x-e^{-x}}{2})^i e^{x(D-i)} [x^n]$
考虑对 $frac{e^x-e^{-x}}{2}$ 用二项式定理，然后消掉后面的 $e^{-xi}$
会得到 $f_{i}=frac{inom{D}{i}n!}{2^i} sum_{j=0}^{i} (-1)^j e^{xD-2xj}[x^n]$
根据 $e^x=sum_{i=0}^{infty} frac{x^i}{i!}$，第 $n$ 项系数为 $frac{x^n(D-2j)^n}{n!}$
那么就得到 $f_{i}=frac{inom{D}{i}n!}{2^i} sum_{j=0}^{i} (-1)^j frac{x^n (D-2j)^n}{n!}$
这个可以写成卷积的形式，然后用 NTT 算就行，最后那个二项式反演也可以用 NTT 来加速.
在程序的开头可以特判答案为 0 和答案为 $D^n$ 的情况.
时间复杂度为 $O(D log D)$.

#include <cstdio>  
#include <vector>
#include <cstring>
#include <algorithm>    
#define N 100009  
#define ll long long 
#define mod 998244353
#define setIO(s) freopen(s".in","r",stdin) 
using namespace std;   
int fac[N],inv[N],A[N<<2],B[N<<2],f[N],g[N],lw[N];    
int ADD(int x,int y) { 
    return (ll)(x+y)%mod;  
} 
int DEC(int x,int y) { 
    return (ll)(x-y+mod)%mod; 
} 
int MUL(int x,int y) {  
    return (ll)x*y%mod; 
}
int qpow(int x,int y) { 
    int tmp=1; 
    for(;y;y>>=1,x=(ll)x*x%mod) 
        if(y&1) {   
            tmp=(ll)tmp*x%mod; 
        } 
    return tmp; 
} 
int get_inv(int x) { 
    return qpow(x,mod-2); 
}  
void NTT(int *a,int len,int op) { 
    for(int i=0,k=0;i<len;++i) { 
        if(i>k) { 
            swap(a[i],a[k]); 
        } 
        for(int j=len>>1;(k^=j)<j;j>>=1); 
    }  
    for(int l=1;l<len;l<<=1) { 
        int wn=qpow(3,(mod-1)/(l<<1)); 
        if(op==-1) {    
            wn=get_inv(wn); 
        } 
        for(int i=0;i<len;i+=l<<1) { 
            int w=1,x,y;  
            for(int j=0;j<l;++j) {  
                x=a[i+j],y=(ll)w*a[i+j+l]%mod;  
                a[i+j]=(ll)(x+y)%mod; 
                a[i+j+l]=(ll)(x-y+mod)%mod;   
                w=(ll)w*wn%mod;  
            }
        }
    }
    if(op==-1) { 
        int iv=get_inv(len); 
        for(int i=0;i<len;++i) {    
            a[i]=(ll)a[i]*iv%mod;  
        }
    }
}
void init() {
    fac[0]=lw[0]=1; 
    lw[1]=get_inv(2);  
    for(int i=1;i<N;++i) { 
        fac[i]=(ll)fac[i-1]*i%mod; 
        lw[i]=(ll)lw[i-1]*lw[1]%mod;  
    }  
    inv[1]=1; 
    for(int i=2;i<N;++i) { 
        inv[i]=(ll)(mod-mod/i)*inv[mod%i]%mod; 
    } 
    inv[0]=1; 
    for(int i=1;i<N;++i) inv[i]=(ll)inv[i-1]*inv[i]%mod; 
}
int C(int x,int y) { 
    return (ll)fac[x]*inv[y]%mod*inv[x-y]%mod;   
}
int main() { 
    // setIO("input");      
    int n,m,D,lim;  
    scanf("%d%d%d",&D,&n,&m);     
    if(n-2*m<0) {  
        printf("0
");  
        return 0; 
    } 
    if(n-2*m>=D) { 
        printf("%d
",qpow(D,n)); 
        return 0; 
    }  
    init();  
    for(lim=1;lim<=(D<<1);lim<<=1);  
    for(int i=0;i<=D;++i) { 
        int d=(i&1)?mod-1:1;   
        A[i]=(ll)d*inv[i]%mod*qpow(DEC(D,2*i),n)%mod;   
        B[i]=inv[i];  
    }
    NTT(A,lim,1),NTT(B,lim,1);  
    for(int i=0;i<lim;++i) {    
        A[i]=(ll)A[i]*B[i]%mod; 
    } 
    NTT(A,lim,-1);  
    for(int i=0;i<=D;++i) { 
        f[i]=(ll)A[i]*C(D,i)%mod*fac[i]%mod*lw[i]%mod;  
    }    
    for(int i=0;i<lim;++i) {    
        A[i]=B[i]=0; 
    }
    for(int i=0;i<=D;++i) {  
        int d=(i&1)?mod-1:1;  
        A[i]=(ll)f[D-i]*fac[D-i]%mod;   
        B[i]=(ll)d*inv[i]%mod;   
    }
    NTT(A,lim,1),NTT(B,lim,1); 
    for(int i=0;i<lim;++i) { 
        A[i]=(ll)A[i]*B[i]%mod;  
    }   
    NTT(A,lim,-1);    
    int ans=0;   
    for(int i=0;i<=n-2*m;++i) { 
        ans=ADD(ans,(ll)inv[i]*A[D-i]%mod);  
    }         
    printf("%d
",ans); 
    return 0; 
}