ZOJ 3899 State Reversing【NTT+斯特林数+组合】

/*
*Rainto96
*Beijing University of Posts and Telecommunications School of Software Engineering
*http://blog.csdn.net/u011775691
题意
Alice有n1个可区分的糖，Bob有n2个不可区分的糖
Alice的盒子是不可区分的，当Alice有n个盒子时候，如果n是质数，那么Alice把糖果随意分但是保证每个盒子有一个糖。
如果n是合数，那么Alice先把每个盒子放一个糖，剩下的都丢到一个盒子里

Bob的盒子也是不可区分的，当Bob有n个盒子时，但是Bob有M个颜料（颜料数目保证大于等于他有的盒子数目），首先Bob先选
n个不同颜色把盒子涂上不同的颜色，然后再放糖果

现在有Q个询问，给出盒子总数N，问糖果放置方案数

解法：http://blog.csdn.net/u013007900/article/details/48489083
男神的博客写的很好

当Alice有质数个盒子时候，相当于是第二类斯特林数
当Alice有合数个盒子时候，相当于先选和盒子数相同的个数的糖

Bob的情况很简单，先选出n种颜色涂盒子，然后再插板法

很容易看出答案Ans(N) = sigma(alice(i) * bob(N-i))
然后考虑题目给的模数，用NTT加速即可
*/
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <queue>
#include <cstring>
#include <string>
#include <cmath>
#include <set>
#include <map>
#include <vector>
#include <climits>
using namespace std;
#define pb push_back
#define ALL(x) x.begin(),x.end()
#define VINT vector<ll>
#define PII pair<ll,ll>
#define MP(x,y) make_pair((x),(y))
#define ll long long
#define ull unsigned ll
#define MEM0(x)  memset(x,0,sizeof(x))
#define MEM(x,val) memset((x),val,sizeof(x))
#define scan(x) scanf("%lld",&(x))
#define scan2(x,y) scanf("%lld%lld",&(x),&(y))
#define scan3(x,y,z) scanf("%lld%lld%lld",&(x),&(y),&(z))
#define scan4(x,y,z,k) scanf("%lld%lld%lld%lld",&(x),&(y),&(z),&(k))
#define Max(a,b) a=max(a,b)
#define Min(a,b) a=min(a,b)
using namespace std;
const ll MOD = 786433;
ll n1,n2,M,Q;

const ll UP = 18888;
ll prm[UP + 10];
bool visP[UP + 10];
bool isprm[UP + 10];
ll pn;
void getPrm(){
    for(ll i=2;i<=UP;i++){
        if(!visP[i]){
            prm[pn++] = i;
            isprm[i] = true;
            for(ll j=i;j<=UP;j+=i) visP[j] = true;
        }
    }
}

ll qPow(ll a, ll n, ll MOD){
    ll sum = 1;
    ll now = a;
    while(n){
        if(n&1) sum = sum * now % MOD;
        now = now * now %MOD;
        n >>= 1;
    }
    return sum % MOD;
}
namespace NTT{
    const ll r=11,gl=20;
    ll p, rp[50], irp[50];
    void setMOD(ll _p = 786433){
        p = _p;
        for(ll i=0;i<gl;i++) rp[i] = qPow(r, (p-1)/(1<<i), p);
    }
    void FFT(ll a[], ll n, ll wt[]=rp){
        for(ll i=0, j=0; i<n;i++){
            if(j>i) swap(a[i], a[j]);
            ll k = n;
            while(j&(k>>=1)) j&=~k;
            j |= k;
        }
        for(ll m=1,b=1; m<n; m<<=1, b++){
            for(ll k=0,w=1; k<m;k++){
                for(ll i = k; i<n;i+=m<<1){
                    ll v = a[i+m] *w %p;
                    if((a[i+m] = a[i] - v) < 0) a[i+m] += p;
                    if((a[i] += v) >= p) a[i] -= p;
                }
                w = w * wt[b] % p;
            }
        }
    }
    void IFFT(ll a[], ll n){
        for(ll i=0;i<gl;i++) irp[i] = qPow(rp[i], n-1, p);
        FFT(a, n, irp);
        ll inv = qPow(n,p-2,p);
        for(ll i=0;i<n;i++) a[i] = a[i] * inv % p;
    }
    void Mul(ll a[], ll b[], ll n ,ll c[]){
        ll len = 1;
        while(len < n) len<<=1;
        len<<=1;
        FFT(a, len);
        FFT(b, len);
        for(ll i=0;i<len;i++) c[i] = a[i] * b[i] % p ;
        IFFT(c, len);
    }
}
const ll MAXN = 111111;
ll a[MAXN] , b[MAXN], c[MAXN];
ll stirling[MAXN];
ll fac[MAXN] , inv[MAXN];
void init(){
    fac[0] = 1;
    inv[0] = qPow(1,MOD-2,MOD);
    for(ll i=1;i<11111;i++) fac[i] = fac[i-1] * i % MOD;
    for(ll i=1;i<11111;i++) inv[i] = qPow(fac[i] , MOD-2 , MOD);
}
int main(){
	#ifndef ONLINE_JUDGE
		//freopen("C:/OJ/in.txt","r",stdin);
	#endif
	getPrm();
	init();



    ll T;scan(T);
	while(T--){
        scan4(n1,n2,M,Q);
        memset(a, 0, sizeof(a));
        memset(b, 0, sizeof(b));
        NTT::setMOD();
        for(ll i=0 ,fuyi=1;i<=n1;i++, fuyi*=-1){
            a[i] = (inv[i]*fuyi % MOD + MOD) %MOD;
            b[i] = qPow(i,n1,MOD) * inv[i] % MOD;
        }
        NTT::Mul(a,b,n1+1,stirling);

        memset(a, 0, sizeof(a));
        memset(b, 0, sizeof(b));
        for(ll i=1;i<=n1 && i<=M;i++){
            if(isprm[i]){
                a[i] = stirling[i];
            }else if(i == n1) a[i] = 1;当物品个数与箱子个数相同时!!!!!!!!!!!!!!!!!!!!!!
            else{
                a[i] = fac[n1] * inv[i-1]%MOD * inv[n1-i+1]%MOD;
            }
        }
        for(ll i=1;i<=n2 && i<=M;i++){
            b[i] = fac[M] * inv[i]%MOD * inv[M-i] % MOD;
            b[i] = b[i] * fac[n2-1]%MOD * inv[i-1]%MOD * inv[n2-i] % MOD;
        }
        NTT::Mul(a,b,M,c);
        while(Q--){
            ll qr;scan(qr);
            printf("%lld\n",c[qr]);
        }
	}
    return 0;
}