洛谷 P2791 幼儿园篮球题

---

https://www.luogu.org/problemnew/show/P2791

我喜欢唱♂跳♂rap♂篮球

要求的是：\(\sum_{i=0}^kC_m^iC_{n-m}^{k-i}i^L\)

这个\(i^L\)很烦，就把第二类斯特林数的式子套进去

\(\sum_{i=0}^kC_m^iC_{n-m}^{k-i}i^L\)

\(\sum_{i=0}^kC_m^iC_{n-m}^{k-i}\sum_{j=0}^iC_{i}^j\begin{Bmatrix}L\\j\end{Bmatrix}j!\)

\(\sum_{j=0}^k\begin{Bmatrix}L\\j\end{Bmatrix}j!\sum_{i=j}^kC_m^iC_{n-m}^{k-i}C_{i}^j\)

后面\(\sum\)三个组合数好像很不好搞，但是\(C_{m}^{i}C_{i}^{j}=C_{m}^{j}C_{m-j}^{i-j}\)，可以拆出一个与\(i\)无关的组合数

\(\sum_{j=0}^k\begin{Bmatrix}L\\j\end{Bmatrix}j!C_{m}^{j}\sum_{i=j}^kC_{m-j}^{i-j}C_{n-m}^{k-i}\)

把式子化的好看一点，发现可以套范德蒙德卷积（\(\sum_{i=0}^kC_{n}^{i}C_{m}^{k-i}=C_{n+m}^k\)）

\(\sum_{j=0}^k\begin{Bmatrix}L\\j\end{Bmatrix}j!C_{m}^{j}\sum_{i=0}^kC_{m-j}^{k-i-j}C_{n-m}^{i}\)

\(\sum_{j=0}^k\begin{Bmatrix}L\\j\end{Bmatrix}j!C_{m}^{j}C_{n-j}^{k-j}\)

注意上面循环\(i\)的上界实际上是\(\min(k,L,m)=O(L)\)

求出\(n=L\)的一行第二类斯特林数，每次询问就可以\(O(L)\)了

把组合数全拆出来约分就洛谷rk1了= =

#include<bits/stdc++.h>
#define il inline
#define vd void
#define mod 998244353
#define poly std::vector<int>
typedef long long ll;
il ll gi(){
    ll x=0,f=1;
    char ch=getchar();
    while(!isdigit(ch))f^=ch=='-',ch=getchar();
    while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
    return f?x:-x;
}
il int pow(int x,int y){
    int ret=1;
    while(y){
        if(y&1)ret=1ll*ret*x%mod;
        x=1ll*x*x%mod;y>>=1;
    }
    return ret;
}
#define maxn 524289
poly pA,pB;
int rev[maxn],_lstN,P[maxn],iP[maxn];
il vd ntt(int*A,int N,int t){
    for(int i=0;i<N;++i)if(rev[i]>i)std::swap(A[i],A[rev[i]]);
    for(int o=1;o<N;o<<=1){
        int W=t?P[o]:iP[o];
        for(int*p=A;p!=A+N;p+=o<<1)
            for(int i=0,w=1;i<o;++i,w=1ll*w*W%mod){
                int t=1ll*w*p[i+o]%mod;
                p[i+o]=(p[i]-t+mod)%mod;p[i]=(p[i]+t)%mod;
            }
    }
    if(!t){
        int inv=pow(N,mod-2);
        for(int i=0;i<N;++i)A[i]=1ll*A[i]*inv%mod;
    }
}
int N,lg;
il vd setN(int n){
    N=1,lg=0;
    while(N<n)N<<=1,++lg;
    if(N!=_lstN)for(int i=0;i<N;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<lg-1);
}
il vd ntt(poly&a,int t){
    static int A[maxn];
    for(int i=0;i<a.size();++i)A[i]=a[i];memset(A+a.size(),0,4*(N-a.size()));
    ntt(A,N,t);
    a.resize(N);
    for(int i=0;i<N;++i)a[i]=A[i];
    int s=a.size();while(s&&!a[s-1])--s;
    a.resize(s);
}
il poly mul(poly a,poly b,int newn=-1){
    if(newn==-1)newn=a.size()+b.size()-1;
    setN(a.size()+b.size()-1);
    ntt(a,1),ntt(b,1);
    for(int i=0;i<N;++i)a[i]=1ll*a[i]*b[i]%mod;
    ntt(a,0);a.resize(newn);
    return a;
}
il poly operator+(poly a,const poly&b){
    if(a.size()<b.size())a.resize(b.size());
    for(int i=0;i<a.size();++i)if(i<b.size())a[i]=(a[i]+b[i])%mod;
    return a;
}
il poly operator-(poly a,const poly&b){
    if(a.size()<b.size())a.resize(b.size());
    for(int i=0;i<a.size();++i)if(i<b.size())a[i]=(a[i]-b[i]+mod)%mod;
    return a;
}
il poly operator*(poly a,int b){
    for(auto&i:a)i=1ll*i*b%mod;
    return a;
}
il poly qiudao(poly a){
    for(int i=0;i<a.size()-1;++i)a[i]=1ll*a[i+1]*(i+1)%mod;
    a.erase(a.end()-1);
    return a;
}
il poly jifen(poly a){
    a.insert(a.begin(),0);
    for(int i=1;i<a.size();++i)a[i]=1ll*a[i]*pow(i,mod-2)%mod;
    return a;
}
il poly getinv(poly a){
    if(a.size()==1)return poly(1,pow(a[0],mod-2));
    int n=a.size(),m=a.size()+1>>1;
    poly _a(m);
    for(int i=0;i<m;++i)_a[i]=a[i];
    poly b=getinv(_a);
    setN(n+m*2-2);
    ntt(a,1);ntt(b,1);
    for(int i=0;i<N;++i)a[i]=1ll*a[i]*b[i]%mod*b[i]%mod;
    ntt(a,0),ntt(b,0);
    a.resize(n);
    return b*2-a;
}
il poly getln(poly a,int n=-1){
    if(n==-1)n=a.size();
    a.resize(n);
    return jifen(mul(qiudao(a),getinv(a),n));
}
il poly getexp(poly a){
    if(a.size()==1)return a[0]=1,a;
    int n=a.size(),m=a.size()+1>>1;
    poly _a(m);
    for(int i=0;i<m;++i)_a[i]=a[i];
    poly b=getexp(_a);
    return mul(b,poly(1,1)-getln(b,a.size())+a,a.size());
}
il poly operator^(poly a,int b){
    int n=a.size();
    a=getexp(getln(a)*b);a.resize(n);
    return a;
}
il poly sqrt(poly a){
    if(a.size()==1)return a;
    int n=a.size(),m=a.size()+1>>1;
    poly _a(m);
    for(int i=0;i<m;++i)_a[i]=a[i];
    poly b=sqrt(_a);b.resize(n);
    return (b+mul(a,getinv(b),n))*(mod+1>>1);
}
il vd poly_init(){
    int G=3,iG=332748118;
    for(int i=1;i<maxn;i<<=1)P[i]=pow(G,(mod-1)/(i<<1)),iP[i]=pow(iG,(mod-1)/(i<<1));
}
il vd add(int&a,int b){a=a+b>=mod?a+b-mod:a+b;}
int fact[20000010],ifact[20000010];
int pr[4000010],Pr,PL[20000010];
bool yes[20000010];
int main(){
#ifdef XZZSB
    freopen("in.in","r",stdin);
    freopen("out.out","w",stdout);
#endif
    poly_init();
    int N=gi(),M=gi(),S=gi(),L=gi(),o=std::max(N,L);
    fact[0]=1;for(int i=1;i<=o;++i)fact[i]=1ll*i*fact[i-1]%mod;
    ifact[o]=pow(fact[o],mod-2);for(int i=o-1;~i;--i)ifact[i]=1ll*ifact[i+1]*(i+1)%mod;
    PL[0]=0,PL[1]=1;
    for(int i=2;i<=L;++i){
        if(!yes[i])pr[++Pr]=i,PL[i]=pow(i,L);
        for(int j=1;j<=Pr&&i*pr[j]<=L;++j){
            yes[i*pr[j]]=1;
            PL[i*pr[j]]=1ll*PL[pr[j]]*PL[i]%mod;
            if(i%pr[j]==0)break;
        }
    }
    poly f(L+1),g(L+1);
    for(int i=0,mul=1;i<=L;++i,mul=mod-mul)f[i]=1ll*mul*ifact[i]%mod;
    for(int i=0;i<=L;++i)g[i]=1ll*PL[i]*ifact[i]%mod;
    f=mul(f,g,L+1);
    while(S--){
        int n=gi(),m=gi(),k=gi(),ans=0,o=std::min(k,std::min(m,L));
        for(int i=0;i<=o;++i)
            ans=(ans+1ll*f[i]*ifact[m-i]%mod*fact[n-i]%mod*ifact[k-i])%mod;
        printf("%d\n",1ll*ans*fact[m]%mod*fact[k]%mod*ifact[n]%mod);
    }
    return 0;
}