题意

给定正整数 $n,m\le 2\times 10^5$，有一枚质地均匀的，$n$ 个面的骰子。对于 $i=1,2\dots m$ 分别求出，不断扔这个骰子，直到每一面都出现过，所需步数的 $i$ 次方的期望。模数 $998244353$。

题解

考虑 $k$ 次方的组合意义，即任意选 $k$ 个时刻的方案数，满足：对于任意一个时刻，都至少存在一个面，没有被扔出来过。

考虑枚举这 $k$ 个时刻的最大值 $i$

此时枚举最大值第一次出现的位置 $t+1$，方案数即为 $\displaystyle \sum_{t=0}^{k-1}i^t\left(i+1\right)^{k-1-t}=\sum_{t=0}^{k-1}\binom{k}{t}i^t$

由于 $0,1,\dots,i$ 时刻至少存在一个面没有被扔出来过，因此概率为 $1-\left\{\begin{matrix}i\\n\end{matrix}\right\}/n^i$

因此步数的 $k$ 次方的期望 $\displaystyle a_k=\sum_{i\ge 0}\left(1-\left\{\begin{matrix}i\\n\end{matrix}\right\}/n^i\right)\sum_{t=0}^{k-1}\binom{k}{t}i^t$

令 $\displaystyle b_t=\sum_{i\ge 0}\left(1-\left\{\begin{matrix}i\\n\end{matrix}\right\}/n^i\right)i^t$

展开可得

$$ \begin{aligned} b_t&=\sum_{i\ge 0}\sum_{j=0}^{n-1}\left(-1\right)^{n-j-1}\binom{n}{j}\left(\frac{j}{n}\right)^ii^t\\ &=\sum_{j=0}^{n-1}\left(-1\right)^{n-j-1}\binom{n}{j}\left(\left[\dfrac{x^t}{t!}\right]\dfrac{1}{1-\frac{j}{n}e^x}\right) \end{aligned} $$

令 $b_t$ 对应的 EGF 为 $B_0$，OGF 为 $B_1$，$G=\sum_{j=0}^{n-1}\left(-1\right)^{n-j-1}\binom{n}{j}\frac{1}{1-jx/n}$，则 $B_0=G\circ e^x$，令 $G_0=G\circ \left(x+1\right)$

则 $B=G_0\circ \left(e^x-1\right)$，此时可以直接分治乘求出 $G_1\equiv G_0\pmod{x^m}$

那么 $B_0\equiv \left(G_1\circ \left(x-1\right) \circ e^x\right)\pmod {x^m}$，令 $G_1\circ \left(x-1\right)=\sum_{i=0}^{m-1}g_ix^i$

则 $\displaystyle B_1\equiv \sum_{i=0}^{m-1}\dfrac{g_i}{1-ix}\pmod{x^m}$，再用一次分治乘即可得到 $b_t$

由于 $a_k=\sum_{t=0}^{k-1}\binom{k}{t}b_t$，再做一次卷积即可。

Code

#include <iostream>
#include <algorithm>
#include <cstring>
#include <vector>
using ll=long long;
namespace polynomial
{
    constexpr int P(998244353),N(200010),PR(3),APR((P+1)/PR);
    int ws[2][N<<2],fac[N<<2],ifac[N<<2],fn,fb,lgg[N<<2],_inv[N<<2],rev[N<<2];
    int fast_pow(int a,int b){int ans(1),off(a);while(b){if(b&1) ans=1ll*ans*off%P;off=1ll*off*off%P;b>>=1;}return ans;}
    void init_poly_weights(int len)
    {
        fn=1,fb=0;while(fn<(len<<1)) fn<<=1,++fb;
        _inv[1]=1;for(int i(2);i<=fn;++i) _inv[i]=1ll*(P-P/i)*_inv[P%i]%P;
        ifac[0]=fac[0]=1;for(int i(1);i<=fn;++i) fac[i]=1ll*i*fac[i-1]%P,ifac[i]=1ll*ifac[i-1]*_inv[i]%P;
        for(int i(0);i<fn;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(fb-1));
        for(int i(2);i<=fn;++i) lgg[i]=lgg[(i+1)>>1]+1;
        for(int mid((fn>>1)),j0(fast_pow(PR,(P-1)/(mid<<1))),j1(fast_pow(APR,(P-1)/(mid<<1)));mid>=1;mid>>=1,j0=1ll*j0*j0%P,j1=1ll*j1*j1%P)
        for(int i(0),w0(1),w1(1);i<mid;++i,w0=1ll*w0*j0%P,w1=1ll*w1*j1%P) ws[0][i+mid]=w0,ws[1][i+mid]=w1;
    }
    using poly=std::vector<int>;
    int redu(int k){return k>=P?k-P:k;}
    poly operator+(const poly &a,const poly &b)
    {
        poly ret(b);if(ret.size()<a.size()) ret.resize(a.size());
        for(int i(0);i<a.size();++i) ret[i]=redu(ret[i]+a[i]);
        return ret;
    }
    unsigned long long tt[N<<2];
    void NTT(poly &p,int V)
    {
        int bts(lgg[p.size()]);if(p.size()!=(1<<bts)) p.resize((1<<bts));
        int* w(ws[(V==1)?0:1]),len(1<<bts);for(int i(0);i<len;++i) tt[i]=p[rev[i]>>(fb-bts)];
        unsigned long long t1,t2;
        for(int l(2);l<=len;l<<=1)
            for(int j(0),mid=(l>>1);j<len;j+=l)
                for(int i(0);i<mid;++i) t1=tt[j+i],t2=tt[j+i+mid]*w[mid+i]%P,tt[j+i]=t1+t2,tt[j+i+mid]=P+t1-t2;
        if(V==1) for(int i(0);i<len;++i) p[i]=tt[i]%P;
        else for(int i(0),j(_inv[len]);i<len;++i) p[i]=tt[i]%P*j%P;
    }
    poly operator*(const poly &a,const poly &b)
    {
        static poly p1,p2;poly ret;
        p1=a,p2=b;int len(a.size()+b.size()-1),ff(1<<lgg[len]);
        p1.resize(ff),p2.resize(ff),ret.resize(ff);
        NTT(p1,1);NTT(p2,1);
        for(int i(0);i<ff;++i) ret[i]=1ll*p1[i]*p2[i]%P;
        NTT(ret,-1);ret.resize(len);return ret;
    }
    void print_poly(const poly &a){for(int v:a) printf("%d ",v);printf("\n");}
    poly poly_inv(const poly &a)
    {
        int l(a.size());if(l==1){poly ret(1);ret[0]=fast_pow(a[0],P-2);return ret;}
        poly g0(a);g0.resize((l+1)>>1);g0=poly_inv(g0);
        static poly p1;p1=a;int ff(2<<lgg[l]);poly ret(ff);g0.resize(ff);p1.resize(ff);
        NTT(p1,1);NTT(g0,1);for(int i(0);i<ff;++i) ret[i]=1ll*g0[i]*(2-1ll*g0[i]*p1[i]%P+P)%P;
        NTT(ret,-1);ret.resize(l);return ret;
    }
    poly poly_ln(const poly &a)
    {
        int l(a.size());static poly p1;
        p1.resize(l-1);for(int i(1);i<l;++i) p1[i-1]=1ll*i*a[i]%P;
        p1=p1*poly_inv(a);p1.resize(l-1);poly ret(l);
        for(int i(1);i<l;++i) ret[i]=1ll*_inv[i]*p1[i-1]%P;ret[0]=0;
        return ret;
    }
    poly poly_exp(const poly &a)
    {
        int l(a.size());if(l==1){poly ret(1);ret[0]=1;return ret;}
        poly g0(a);g0.resize((l+1)>>1);g0=poly_exp(g0);static poly g1,g2;g1=g0;g1.resize(l);g1=poly_ln(g1);g2=a;
        int ff(2<<lgg[l]);g0.resize(ff);g1.resize(ff);poly ret(ff);g2.resize(ff);
        NTT(g0,1);NTT(g1,1);NTT(g2,1);
        for(int i(0);i<ff;++i) ret[i]=1ll*g0[i]*(1-g1[i]+g2[i]+P)%P;
        NTT(ret,-1);ret.resize(l);return ret;
    }
    std::pair<int,int> cf[N];// u/(1+dx)
    std::pair<poly,poly> calc(int l,int r)
    {
        if(l==r)
        {
            auto [u,d]=cf[l];
            return {{u},{1,d}};
        }
        int mid((l+r)>>1);
        auto [ls,lt]=calc(l,mid);
        auto [rs,rt]=calc(mid+1,r);
        int T(1<<lgg[r-l+2]);ls.resize(T);lt.resize(T);rs.resize(T);rt.resize(T);
        poly s(T),t(T);
        NTT(ls,1);NTT(lt,1);NTT(rs,1);NTT(rt,1);
        for(int i(0);i<T;++i) s[i]=(1ll*ls[i]*rt[i]+1ll*rs[i]*lt[i])%P,t[i]=1ll*lt[i]*rt[i]%P;
        NTT(s,-1);NTT(t,-1);s.resize(r-l+1);t.resize(r-l+2);
        return {s,t};
    }
    poly solve_fs(const std::vector<std::pair<int,int>> &a,int m)
    {
        int n(a.size());
        std::copy(a.begin(),a.end(),cf+1);
        auto [s,t]=calc(1,n);
        t.resize(m+1);
        t=poly_inv(t);
        s=s*t;
        s.resize(m+1);
        return s;
    }
    const poly unit_poly={1};
}
using polynomial::poly;
using polynomial::operator*;
using polynomial::operator+;
using polynomial::init_poly_weights;
using polynomial::unit_poly;
using polynomial::print_poly;
using polynomial::poly_inv;
using polynomial::poly_exp;
using polynomial::poly_ln;
using polynomial::solve_fs;
constexpr int N(200010),p(998244353);
int main()
{
    int n,m;std::cin>>n>>m;
    init_poly_weights(std::max(n+1,m));
    static int inv[N],fac[N],ifac[N];inv[1]=1;
    for(int i(2);i<N;++i) inv[i]=1ll*(p-p/i)*inv[p%i]%p;
    fac[0]=ifac[0]=1;for(int i(1);i<N;++i) fac[i]=1ll*fac[i-1]*i%p,ifac[i]=1ll*ifac[i-1]*inv[i]%p;
    auto binom=[](int a,int b)->int{return 1ll*fac[a]*ifac[b]%p*ifac[a-b]%p;};
    std::vector<std::pair<int,int>> c(n);
    for(int j(0);j<n;++j)
    {
        c[j].second=1ll*(p-j)*inv[n-j]%p;
        c[j].first=1ll*(((n-j)&1)?(binom(n,j)):(p-binom(n,j)))*n%p*inv[n-j]%p;
    }
    auto t=solve_fs(c,m-1);
    for(int i(0);i<m;++i) t[i]=1ll*t[i]*fac[i]%p;
    std::reverse(t.begin(),t.end());
    poly df(m);
    for(int i(0);i<m;++i) df[i]=(i&1)?(p-ifac[i]):ifac[i];
    t=t*df;t.resize(m);
    std::reverse(t.begin(),t.end());
    c.resize(m);
    for(int i(0);i<m;++i)
    {
        t[i]=1ll*t[i]*ifac[i]%p;
        c[i].first=t[i];
        c[i].second=p-i;
    }
    t=solve_fs(c,m-1);
    poly tr(ifac+1,ifac+m+1);
    for(int i(0);i<m;++i) t[i]=1ll*t[i]*ifac[i]%p;
    t=t*tr;t.resize(m);
    for(int i(1);i<=m;++i) std::cout<<1ll*t[i-1]*fac[i]%p<<"\n";
    return 0;
}