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$个盒子，那么有

$$ans=\sum\limits_{num=1}^{N-1}Alice(num)*Bob(n-num)$$ 。
这明显是卷积形式。

分类讨论，得出：

$$Alice(i)=\begin{cases}S_{N1}^i,i\in P\\\\\binom{N1}{i-1},i\notin p,i\ne N1\\\\1,i\notin p,i= N1\end{cases}$$
$$Bob(i)=\binom{M}{i}\binom{N2-1}{i-1}$$
其中 $S_{N1}^i$为第二类斯特林数，也可以用NTT快速求出一行的值，详见 这里。

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

}
void pre_S()
{
    memset(a,0,sizeof(a)),memset(b,0,sizeof(b)),memset(S,0,sizeof(S));
    for(int i=0;i<=n1;i++)a[i]=((i&1)?mod-1:1)*fac_inv[i]%mod,b[i]=Pow(i,n1)*fac_inv[i]%mod;
    int len=1;while(len<=2*n1)len<<=1;
    NTT(a,len,1),NTT(b,len,1);
    for(int i=0;i<len;i++)a[i]=a[i]*b[i]%mod;
    NTT(a,len,-1);
    for(int i=0;i<=n1;i++)S[i]=a[i];
}
void solve()
{
    memset(a,0,sizeof(a)),memset(b,0,sizeof(b)),memset(ans,0,sizeof(ans));
    for(int i=1;i<=min(n1,m);i++)a[i]=notpri[i]?(i==n1?1:C(n1,i-1)):S[i];
    for(int i=1;i<=min(n2,m);i++)b[i]=C(m,i)*C(n2-1,i-1)%mod;
    int len=1;while(len<=2*m)len<<=1;
    NTT(a,len,1),NTT(b,len,1);
    for(int i=0;i<len;i++)a[i]=a[i]*b[i]%mod;
    NTT(a,len,-1);
    for(int i=1;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 0;
}