BJ模拟 生日禮物【NTT+斯特林數+組合數學】

題目描述：

今天是Jane的生日。Alice和Bob都有一些糖果，于是這兩個人就去買了N個白色的盒子去包裝這些糖果作為Jane的生日禮物。他們将随機地把這些盒子分成兩堆，一堆給Alice，一堆給Bob（每堆至少有一個盒子）。

我們知道Alice有 N1 個不同的糖果，Bob有 N2 個相同的糖果（由于Bob很懶，是以他直接買了相同的糖果），然後Alice買的糖果和Bob買的糖果是完全不一樣的。他們兩人各自把糖果放進盒子裡，并且任意一個盒子都不會是空的。Alice和Bob各自有各自包裝糖果的規則：

Alice很有個性，她會根據盒子的數量來決定包裝糖果的方法。記Alice擁有的盒子數量為 k 。若 k 是個質數，則她會随機地包裝糖果（每個糖果都有可能進入任意一個盒子裡）。若 k 不是質數，她會往其中的 k−1 個盒子裡分别随機放入1個糖果，将剩下的糖果全部放入剩下那個盒子裡。

Bob有M種顔料（這M種顔料中沒有白色），他将随機地給盒子上色，使得每個盒子都有上色且顔色不一樣。然後Bob會随機地包裝糖果（每個糖果都可能進入任意一個盒子裡）。

問題來了，Jane會收到多少種禮物？答案對786433取模。

兩人把包裝好的盒子湊在一起算一種禮物，二種禮物不同當且僅當有兩種禮物中至少有一個盒子的包裝不同。

N1、N2、N、M<=100000

解題思路：

設Alice拿了 num n u m 個盒子，那麼有

ans=∑num=1N−1Alice(num)∗Bob(n−num) a n s = ∑ n u m = 1 N − 1 A l i c e ( n u m ) ∗ B o b ( n − n u m )

。

這明顯是卷積形式。

分類讨論，得出：

Alice(i)=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪SiN1,i∈P(N1i−1),i∉p,i≠N11,i∉p,i=N1 A l i c e ( i ) = { S N 1 i , i ∈ P ( N 1 i − 1 ) , i ∉ p , i ≠ N 1 1 , i ∉ p , i = N 1

Bob(i)=(Mi)(N2−1i−1) B o b ( i ) = ( M i ) ( N 2 − 1 i − 1 )

其中 SiN1 S N 1 i 為第二類斯特林數，也可以用NTT快速求出一行的值，詳見 這裡。

#include<bits/stdc++.h>
#define ll long long
using namespace std;
int getint()
{
    int i=,f=;char c;
    for(c=getchar();c!='-'&&(c<'0'||c>'9');c=getchar());
    if(c=='-')f=-,c=getchar();
    for(;c>='0'&&c<='9';c=getchar())i=(i<<)+(i<<)+c-'0';
    return i*f;
}
const int N=,mod=,g=;
int T,n1,n2,m,Q;
int pn,pri[N],notpri[N],pos[N<<];
ll fac[N],fac_inv[N],w[N],w_inv[N];
ll a[N<<],b[N<<],S[N],ans[N];
ll Pow(ll x,int y)
{
    ll res=;
    for(;y;y>>=,x=x*x%mod)
        if(y&)res=res*x%mod;
    return res;
}
void sieve()
{
    fac[]=;
    for(int i=;i<N;i++)fac[i]=fac[i-]*i%mod;
    fac_inv[N-]=Pow(fac[N-],mod-);
    for(int i=N-;i>=;i--)fac_inv[i]=fac_inv[i+]*(i+)%mod;
    int len=,num=;
    while(len<=N+N)
    {
        len<<=;
        w[++num]=Pow(g,(mod-)/len);
        w_inv[num]=Pow(w[num],mod-);
    }
    notpri[]=;
    for(int i=;i<N;i++)
    {
        if(!notpri[i])pri[++pn]=i;
        for(int j=;j<=pn;j++)
        {
            ll k=i*pri[j];
            if(k>=N)break;
            notpri[k]=;
            if(i%pri[j]==)break;
        }
    }
}
ll C(int x,int y)
{
    return x>=y?fac[x]*fac_inv[y]%mod*fac_inv[x-y]%mod:;
}
void NTT(ll *f,int len,int on)
{
    for(int i=;i<len;i++)pos[i]=(i&)?pos[i>>]>>|(len>>):pos[i>>]>>;
    for(int i=;i<len;i++)if(i<pos[i])swap(f[i],f[pos[i]]);
    for(int i=,num=;i<len;i<<=,num++)
    {
        ll wi=on==?w[num]:w_inv[num];
        for(int j=;j<len;j+=(i<<))
        {
            ll wn=;
            for(int k=j;k<j+i;k++)
            {
                ll u=f[k],v=wn*f[k+i]%mod;
                f[k]=(u+v)%mod,f[k+i]=(u-v+mod)%mod;
                wn=wn*wi%mod;
            }
        }
    }
    if(on==-)
        for(int i=,len_inv=Pow(len,mod-);i<len;i++)
            f[i]=f[i]*len_inv%mod;

}
void pre_S()
{
    memset(a,,sizeof(a)),memset(b,,sizeof(b)),memset(S,,sizeof(S));
    for(int i=;i<=n1;i++)a[i]=((i&)?mod-:)*fac_inv[i]%mod,b[i]=Pow(i,n1)*fac_inv[i]%mod;
    int len=;while(len<=*n1)len<<=;
    NTT(a,len,),NTT(b,len,);
    for(int i=;i<len;i++)a[i]=a[i]*b[i]%mod;
    NTT(a,len,-);
    for(int i=;i<=n1;i++)S[i]=a[i];
}
void solve()
{
    memset(a,,sizeof(a)),memset(b,,sizeof(b)),memset(ans,,sizeof(ans));
    for(int i=;i<=min(n1,m);i++)a[i]=notpri[i]?(i==n1?:C(n1,i-)):S[i];
    for(int i=;i<=min(n2,m);i++)b[i]=C(m,i)*C(n2-,i-)%mod;
    int len=;while(len<=*m)len<<=;
    NTT(a,len,),NTT(b,len,);
    for(int i=;i<len;i++)a[i]=a[i]*b[i]%mod;
    NTT(a,len,-);
    for(int i=;i<=m;i++)ans[i]=a[i];
}
int main()
{
    //freopen("lx.in","r",stdin);
    sieve();
    for(T=getint();T;T--)
    {
        n1=getint(),n2=getint(),m=getint(),Q=getint();
        pre_S();solve();
        while(Q--)printf("%lld\n",ans[getint()]);
    }
    return ;
}