ZigZagK的博客
[分治NTT+倍增]2022 CCPC 桂林 D【Alice's Dolls】题解
2022年11月1日 17:35
CCPC
查看标签

题目概述

Alice's Dolls

解题报告

定义 $f_{n,k}$ 表示要得到 $n$ 个物品,$x^k$ 的期望。$n>1$ 的情况十分复杂,但是根据期望的线性性,$n>1$ 的情况可以拆成 $E(n)=E(a)+E(b)$ ,再套上 $x^k$ ,可以得到:

$$ E[(x+y)^k]=\sum_{i=0}^{k}{k\choose i}E(x^i)E(y^{k-i})\\ E(n)=E(a)+E(b)\\ f_{n,k}=\sum_{i=0}^{k}{k\choose i}f_{a,i}f_{b,k-i} $$

现在的问题就是求出 $n=1$ 的答案 $f_k$ ,考虑 $k=1$ 的期望,是用错位相减算的,$k>1$ 也可以这么考虑:

$$ f_k=\sum_{i=1}^{\infty}p(1-p)^{i-1}i^k\\ (1-p)f_{k}=\sum_{i=1}^{\infty}p(1-p)^ii^k=\sum_{i=1}^{\infty}p(1-p)^{i-1}(i-1)^k\\ i^k-(i-1)^k=-\sum_{j=0}^{k-1}{k\choose j}i^j(-1)^{k-j}\\ f_k-(1-p)f_k=-p\sum_{i=1}^{\infty}(1-p)^{i-1}\sum_{j=0}^{k-1}{k\choose j}i^j(-1)^{k-j}\\ pf_k=-\sum_{j=0}^{k-1}{k\choose j}(-1)^{k-j}\sum_{i=1}^{\infty}p(1-p)^{i-1}i^j\\ pf_k=-\sum_{j=0}^{k-1}{k\choose j}(-1)^{k-j}f_j $$

初值 $f_0=1$ ,然后就是个分治NTT的形式(也可以多项式求逆)。

得到 $f_{1,k}$ 之后,就可以通过倍增+卷积得到 $f_{n,k}$ 的答案了。

示例程序

#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long LL;
const int maxn=100000,maxt=1<<18,MOD=998244353;

int n,K,A,B,P,IP;
int INV[maxn+5],fac[maxn+5],wn[maxt+5],temA[maxt+5],temB[maxt+5];
int f[maxn+5],g[maxn+5],T,F[maxt+5],ans[maxt+5];

inline int ADD(int x,int y) {return x+y>=MOD?x+y-MOD:x+y;}
inline int MUL(int x,int y) {return (LL)x*y%MOD;}
int Pow(int w,int b) {int s;for (s=1;b;b>>=1,w=MUL(w,w)) if (b&1) s=MUL(s,w);return s;}
__attribute__((constructor)) void NTTPre(){
    int x=Pow(3,(MOD-1)/maxt);
    wn[maxt>>1]=1;
    for (int i=(maxt>>1)+1;i<maxt;i++) wn[i]=MUL(wn[i-1],x);
    for (int i=(maxt>>1)-1;i;i--) wn[i]=wn[i<<1];
    INV[0]=INV[1]=1;for (int i=2;i<=maxn;i++) INV[i]=MUL(MOD-MOD/i,INV[MOD%i]);
    fac[0]=1;for (int i=1;i<=maxn;i++) fac[i]=MUL(fac[i-1],i),INV[i]=MUL(INV[i-1],INV[i]);
}
void NTT(int *a,int n,int f){
    if (f>0){
        for (int k=n>>1;k;k>>=1)
            for (int i=0;i<n;i+=k<<1)
                for (int j=0;j<k;j++){
                    int x=a[i+j],y=a[i+j+k];
                    a[i+j+k]=MUL(x+MOD-y,wn[k+j]);
                    a[i+j]=ADD(x,y);
                }
    } else {
        for (int k=1;k<n;k<<=1)
            for (int i=0;i<n;i+=k<<1)
                for (int j=0;j<k;j++){
                    int x=a[i+j],y=MUL(a[i+j+k],wn[k+j]);
                    a[i+j+k]=ADD(x,MOD-y);
                    a[i+j]=ADD(x,y);
                }
        for (int i=0,INV=MOD-(MOD-1)/n;i<n;i++) a[i]=MUL(a[i],INV);
        reverse(a+1,a+n);
    }
}
void Solve(int L,int R){
    if (L==R) {L?f[L]=MUL(MOD-f[L],IP):f[L]=1;return;}
    int mid=L+(R-L>>1);
    Solve(L,mid);
    int t;for (t=1;t<=R-L;t<<=1);
    for (int i=L;i<=mid;i++) temA[i-L]=f[i];for (int i=mid-L+1;i<t;i++) temA[i]=0;
    temB[0]=0;for (int i=1;i<=R-L;i++) temB[i]=g[i];for (int i=R-L+1;i<t;i++) temB[i]=0;
    NTT(temA,t,1);NTT(temB,t,1);
    for (int i=0;i<t;i++) temA[i]=MUL(temA[i],temB[i]);
    NTT(temA,t,-1);for (int i=mid+1;i<=R;i++) f[i]=ADD(f[i],temA[i-L]);
    Solve(mid+1,R);
}
void Pow(int n){
    if (n==1) {for (int i=0;i<=K;i++) ans[i]=f[i];return;}
    Pow(n>>1);
    NTT(ans,T,1);
    for (int i=0;i<T;i++) ans[i]=MUL(ans[i],ans[i]);
    NTT(ans,T,-1);
    for (int i=K+1;i<T;i++) ans[i]=0;
    if (n&1){
        NTT(ans,T,1);
        for (int i=0;i<T;i++) ans[i]=MUL(ans[i],F[i]);
        NTT(ans,T,-1);
        for (int i=K+1;i<T;i++) ans[i]=0;
    }
}
int main(){
    scanf("%d%d%d%d",&n,&K,&A,&B);
    P=MUL(A,Pow(B,MOD-2));IP=Pow(P,MOD-2);
    for (int i=0,f=1;i<=K;i++,f=MOD-f) g[i]=MUL(f,INV[i]);
    Solve(0,K);
    for (int i=0;i<=K;i++) F[i]=f[i];
    for (T=1;T<=(K<<1);T<<=1);
    NTT(F,T,1);
    Pow(n);
    for (int i=0;i<=K;i++) printf("%d\n",MUL(ans[i],fac[i]));
    return 0;
}
版权声明:本博客所有文章除特别声明外,均采用 CC BY 4.0 CN协议 许可协议。转载请注明出处!