jzoj5858 朝圣 NTT

Description

---

朝圣刚回来，国王就准备举办盛大的宴会，由于你很闲，所以你来帮国王筹备宴会。 现在需要为宴会准备水果，由于你国物质匮乏，所以全国一共只有 n 个苹果，m 个梨。
由于国王具有 “精神错乱” 的特质，所以他将会从这 n + m 个水果中随机拿出 k 个，用于这次宴会。
你听说领主们都有一个癖好，如果宴会中的苹果的数量是一个质数，那么他们就会很 高兴。但是因为国王是随机选的水果，所以你想知道苹果数量是质数的概率，对 998244353 取模。
负责收购梨的官员上周去小舅子家里了，所以现在还不能准确地知道有几个梨，不过 你通过估计可以确定最多有 M 个梨，所以你需要对 m ∈ [0, M ] 都求出答案。

Hint
子任务一（30pts）：n, M ≤ 15，k ≤ 5。
子任务二（20pts）：n, M, k ≤ 10^3。
子任务三（50pts）：无特殊限制。
对于所有的数据，n, M, k ≤ 10^5，k ≤ n。

Solution

---

考虑固定m求答案，显然有$ans[m]=\sum_{i=1}^k [isprime(i)]\frac{\binom{n}{i}\times\binom{m}{k-i}}{\binom{n+m}{k}}$
然后模数特殊考虑NTT，然后我就不会拆了

考虑拆后面那个组合数，令$f_{k-i}=[isprime(i)]\binom{n}{i}\frac{1}{(k-i)!}$，$g_i=\frac{1}{i!}$，然后两个卷起来乘上一个阶乘就是答案了

老年选手的NTT漏洞百出

Code

---

#include <stdio.h>
#include <string.h>
#include <algorithm>
#define rep(i,st,ed) for (int i=st;i<=ed;++i)
#define drp(i,st,ed) for (int i=st;i>=ed;--i)
#define fill(x,t) memset(x,t,sizeof(x))

typedef long long LL;
const int MOD=998244353;
const int N=2000005;

int prime[N],is_prime[N+5],rev[N];

LL f[N],g[N],fac[N],ny[N];

LL ksm(LL x,LL dep) {
    LL ret=1;
    for (;dep;dep>>=1) {
        (dep&1)?(ret=ret*x%MOD):0;
        x=x*x%MOD;
    }
    return ret;
}

void pre_work(int n) {
    rep(i,2,n) is_prime[i]=1;
    for (int i=2;i<=n;++i) {
        if (is_prime[i]) prime[++prime[0]]=i;
        for (int j=1;i*prime[j]<=n&&j<=prime[0];++j) {
            is_prime[i*prime[j]]=0;
            if (i%prime[j]==0) break;
        }
    }
    is_prime[1]=0;
}

void NTT(LL *a,int len,int f) {
    for (int i=1;i<len;++i) if (i<rev[i]) std:: swap(a[i],a[rev[i]]);
    for (int i=1;i<len;i*=2) {
        LL wn;
        if (f==1) wn=ksm(3,(MOD-1)/i/2);
        else wn=ksm(3,MOD-1-(MOD-1)/i/2);
        for (int j=0;j<len;j+=i*2) {
            LL w=1;
            for (int k=0;k<i;++k) {
                LL u=a[j+k],v=a[j+k+i]*w%MOD;
                a[j+k]=(u+v)%MOD; a[j+k+i]=(u-v+MOD)%MOD;
                w=w*wn%MOD;
            }
        }
    }
    if (f==-1) {
        LL tmp=ksm(len,MOD-2);
        for (int i=0;i<len;++i) a[i]=a[i]*tmp%MOD;
    }
}

LL C(LL n,LL m) {
    return (n<m)?(0):(fac[n]*ny[n-m]%MOD*ny[m])%MOD;
}

int main(void) {
    freopen("data.in","r",stdin);
    freopen("myp.out","w",stdout);
    int n,m,k; scanf("%d%d%d",&n,&m,&k);
    pre_work(n+m);
    fac[0]=1; rep(i,1,n+m) fac[i]=1LL*i*fac[i-1]%MOD;
    ny[0]=1; rep(i,1,n+m) ny[i]=ksm(fac[i],MOD-2);
    rep(i,0,m) g[i]=ny[i];
    for (int i=2;i<=k;++i) if (is_prime[i]) {
        f[k-i]=C(n,i)*ny[k-i]%MOD;
    }
    int len,lg; for (len=1,lg=0;len<=n*2;len*=2,++lg);
    for (int i=0;i<len;++i) rev[i]=(rev[i/2]/2)|((i&1)<<(lg-1));
    NTT(f,len,1); NTT(g,len,1);
    for (int i=0;i<len;++i) f[i]=f[i]*g[i]%MOD;
    NTT(f,len,-1);
    rep(i,0,m) printf("%lld\n", f[i]*fac[i]%MOD*ny[n+i]%MOD*fac[k]%MOD*fac[n+i-k]%MOD);
    return 0;
}