题意

有三个变量，初始为 $a,b,c$，每次随机选择一个变量 $+1$，求出操作 $1,2,\cdots n$ 次后最大值的期望。$\bmod\;998244353$

题解

以下定义 $F_{k}=\exp\left(x\right)\bmod x^{k+1}$

首先转化为：操作 $n$ 次最大值小于等于 $k$ 的概率：

$$ 3^{-n}\left[z^n/n!\right]F_{k-a}F_{k-b}F_{k-c} $$

那么需要对所有 $k$ 求和，即求：

$$ \sum_{i=0}^LF_{i+a^{\prime}}F_{i+b^{\prime}}F_{i+c^{\prime}} $$

以下定义：

$$ \begin{aligned} H\left(a,b,c,L\right)&=\sum_{i=0}^LF_{i+a}F_{i+b}F_{i+c}\\ G\left(a,b,c,L\right)&=\sum_{i=0}^L\dfrac{x^{i+a}}{\left(i+a\right)!}F_{i+b}F_{i+c}\\ f\left(a,b,c,L\right)&=\sum_{i=0}^m\dfrac{x^{i+a}}{\left(i+a\right)!}\dfrac{x^{i+b}}{\left(i+b\right)!}F_{i+c}\\ g\left(a,b,c,L\right)&=\sum_{i=0}^L\dfrac{x^{i+a}}{\left(i+a\right)!}\dfrac{x^{i+b}}{\left(i+b\right)!}iF_{i+c}\\ T\left(a,b\right)&=F_aF_b \end{aligned} $$

那么有 d-finite 递推：

$$ \begin{aligned} H^{\prime}=&\left(\sum_{k=0}^{L}F_{k+a}F_{k+b}F_{k+c}\right)^{\prime}\\ =&3H-\sum_{k=0}^L\dfrac{x^{k+a}}{\left(k+a\right)!}F_{k+b}F_{k+c}-\sum_{k=0}^L\dfrac{x^{k+b}}{\left(k+b\right)!}F_{k+a}F_{k+c}-\sum_{k=0}^L\dfrac{x^{k+c}}{\left(k+c\right)!}F_{k+a}F_{k+b}\\ =&3H-G\left(a,b,c,L\right)-G\left(b,a,c,L\right)-G\left(c,a,b,L\right)\\ G^{\prime}\left(a,b,c,L\right)=&\sum_{k=0}^L\dfrac{x^{k+a-1}}{\left(k+a-1\right)!}F_{k+b}F_{k+c}+2G-f\left(a,b,c,L\right)-f\left(a,c,b,L\right)\\ =&3G+\dfrac{x^{a-1}}{\left(a-1\right)!}F_{b-1}F_{c-1}-\dfrac{x^{a+L}}{\left(a+L\right)!}F_{b+L}F_{c+L}+f\left(a-1,b,c,L\right)+f\left(a-1,c,b,L\right)\\&-f\left(a,b,c,L\right)-f\left(a,c,b,L\right)-\sum_{k=0}^L\dfrac{x^{k+a-1}}{\left(k+a-1\right)!}\dfrac{x^{k+b}}{\left(k+b\right)!}\dfrac{x^{k+c}}{\left(k+c\right)!}\\ =&3G+\dfrac{x^{a-1}}{\left(a-1\right)!}F_{b-1}F_{c-1}-\dfrac{x^{a+L}}{\left(a+L\right)!}F_{b+L}F_{c+L}+\left(a/x-1\right)f\left(a,b,c,L\right)+\left(a/x-1\right)f\left(a,c,b,L\right)\\ &+1/x\cdot g\left(a,b,c,L\right)+1/x\cdot g\left(a,c,b,L\right)-\sum_{k=0}^L\dfrac{x^{k+a-1}}{\left(k+a-1\right)!}\dfrac{x^{k+b}}{\left(k+b\right)!}\dfrac{x^{k+c}}{\left(k+c\right)!}\\ T^{\prime}\left(a,b\right)=&\left(F_aF_b\right)^{\prime}\\ =&2T-\dfrac{x^a}{a!}F_b-\dfrac{x^b}{b!}F_a\\ f^{\prime}\left(a,b,c,L\right)=&\dfrac{2}{x}\sum_{i=0}^L\dfrac{x^{i+a}}{\left(i+a\right)!}\dfrac{x^{i+b}}{\left(i+b\right)!}iF_{i+c}+\left(1+\dfrac{a+b}{x}\right)f-\sum_{i=0}^L\dfrac{x^{i+a}}{\left(i+a\right)!}\dfrac{x^{i+b}}{\left(i+b\right)!}\dfrac{x^{i+c}}{\left(i+c\right)!}\\ g^{\prime}\left(a,b,c,L\right)=&\dfrac{2}{x}\sum_{i=0}^{L}\dfrac{x^{i+a}}{\left(i+a\right)!}\dfrac{x^{i+b}}{\left(i+b\right)!}i^2F_{i+c}+\left(1+\dfrac{a+b}{x}\right)g-\sum_{i=0}^L\dfrac{x^{i+a}}{\left(i+a\right)!}\dfrac{x^{i+b}}{\left(i+b\right)!}\dfrac{x^{i+c}}{\left(i+c\right)!}i\\ =&\dfrac{2}{x}\left(\sum_{i=0}^{L}\dfrac{x^{i+a}}{\left(i+a\right)!}\dfrac{x^{i+b}}{\left(i+b\right)!}\left(i+a\right)\left(i+b\right)F_{i+c}-\left(a+b\right)g-abf\right)\\&+\left(1+\dfrac{a+b}{x}\right)g-\sum_{i=0}^L\dfrac{x^{i+a}}{\left(i+a\right)!}\dfrac{x^{i+b}}{\left(i+b\right)!}\dfrac{x^{i+c}}{\left(i+c\right)!}i\\ =&2x\cdot \left(f-\dfrac{x^{L+a}x^{L+b}}{\left(L+a\right)!\left(L+b\right)!}F_{L+c}+\dfrac{x^{a-1}}{\left(a-1\right)!}\dfrac{x^{b-1}}{\left(b-1\right)!}F_{c-1}+\sum_{i=-1}^{L-1}\dfrac{x^{i+a}}{\left(i+a\right)!}\dfrac{x^{i+b}}{\left(i+b\right)!}\dfrac{x^{i+c+1}}{\left(i+c+1\right)!}\right)\\ &+\left(1-\dfrac{a+b}{x}\right)g-\dfrac{2ab}{x}f-\sum_{i=0}^L\dfrac{x^{i+a}}{\left(i+a\right)!}\dfrac{x^{i+b}}{\left(i+b\right)!}\dfrac{x^{i+c}}{\left(i+c\right)!}i \end{aligned} $$

验证代码：

R.<z> = PowerSeriesRing(QQ)

def F(t):
    ans = 0
    for i in range(t+1):
        ans += z^i / factorial(i)
    return ans

def G(a,b,c,L):
    ans = 0
    for i in range(L+1):
        ans += z^(i+a) / factorial(i+a) * F(i + b) * F(i + c)
    return ans

def H(a,b,c,L):
    ans = 0
    for i in range(L+1):
        ans += F(i + a) * F(i + b) * F(i + c)
    return ans

def K(a,b,c,L):
    ans = 0
    for i in range(L+1):
        ans += z^(i+a) / factorial(i+a) * z^(i+b) / factorial(i+b) * z^(i+c) / factorial(i+c)
    return ans

def f(a,b,c,L):
    ans = 0
    for i in range(L+1):
        ans += z^(i+a) / factorial(i+a) * z^(i+b) / factorial(i+b) * F(i+c)
    return ans

def g(a,b,c,L):
    ans = 0
    for i in range(L+1):
        ans += i * z^(i+a) / factorial(i+a) * z^(i+b) / factorial(i+b) * F(i+c)
    return ans

def J(a,b,c,L):
    ans = 0
    for i in range(L+1):
        ans += z^(i+a) / factorial(i+a) * z^(i+b) / factorial(i+b) * z^(i+c) / factorial(i+c) * i
    return ans

a,b,c,L = map(int,input().split())

_G = G(a,b,c,L)

_f = f(a,b,c,L)

_g = g(a,b,c,L)

print("G:",-derivative(_G,z) + 3*_G + z^(a-1)/factorial(a-1)*F(b-1)*F(c-1)
      -z^(a+L)/factorial(a+L)*F(b+L)*F(c+L) + (a/z-1)*f(a,b,c,L) + (a/z-1)*f(a,c,b,L)
      +g(a,b,c,L)/z + g(a,c,b,L)/z - K(a-1,b,c,L))

print("f:",-derivative(_f,z) + 2/z*_g + (1+(a+b)/z)*_f - K(a,b,c,L))

print("g:",-derivative(_g,z) + 2*z*(_f - z^(L+a)/factorial(L+a)*z^(L+b)/factorial(L+b)*F(L+c)
      +z^(a-1)/factorial(a-1)*z^(b-1)/factorial(b-1)*F(c-1) + K(a-1,b-1,c,L))
      +(1-(a+b)/z)*_g - 2*a*b/z*_f - J(a,b,c,L))

代码

#include <iostream>
#include <algorithm>
#include <assert.h>
constexpr int p(998244353),N(2e6+10),_3((p+1)/3);
int inv[N],ifac[N],fac[N];
int fp(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 K(int a,int b,int c,int L,int deg,int coef,int res[])
{
    assert(a>=0);assert(b>=0);assert(c>=0);
    for(int i(0);i<=L && 3*i+a+b+c<=deg;++i) res[3*i+a+b+c]=(res[3*i+a+b+c]+1ll*coef*ifac[i+a]%p*ifac[i+b]%p*ifac[i+c])%p;
}
void J(int a,int b,int c,int L,int deg,int coef,int res[])
{
    assert(a>=0);assert(b>=0);assert(c>=0);
    for(int i(0);i<=L && 3*i+a+b+c<=deg;++i) res[3*i+a+b+c]=(res[3*i+a+b+c]+1ll*coef*ifac[i+a]%p*ifac[i+b]%p*ifac[i+c]%p*i)%p;
}
void F(int a,int b,int c,int deg,int coef,int res[])
{
    assert(a>=0);assert(b>=0);assert(c>=0);coef=1ll*coef*ifac[a]%p*ifac[b]%p;
    for(int i(0);i<=c && i+a+b<=deg;++i) res[i+a+b]=(res[i+a+b]+1ll*coef*ifac[i])%p;
}
void F(int a,int c,int deg,int coef,int res[])
{
    assert(a>=0);assert(c>=0);coef=1ll*coef*ifac[a]%p;
    for(int i(0);i<=c && i+a<=deg;++i) res[i+a]=(res[i+a]+1ll*coef*ifac[i])%p;
}
void FF(int a,int b,int c,int deg,int coef,int res[])
{
    assert(a>=0);assert(b>=0);assert(c>=0);
    int nd(deg-a);if(nd<0) return;static int T[N];
    std::fill(T,T+nd,0);F(b,c,nd-1,p-1,T);F(c,b,nd-1,p-1,T);coef=1ll*coef*ifac[a]%p;
    res[a]=(res[a]+coef)%p;
    for(int i(1),t(1);i<=nd;++i)
    {
        t=(2ll*t+T[i-1])*inv[i]%p;
        res[a+i]=(res[a+i]+1ll*t*coef)%p;
    }
}
void fg(int a,int b,int c,int L,int deg,int fx[],int gx[])
{
    assert(a>0);assert(b>0);assert(c>0);if(deg<0) return;
    int numf(0),numg(0);
    for(int i(0);i*2+a+b<=deg && i<=L;++i)
    {
        int j(deg-i*2-a-b);
        if(j<=i+c) numf=(numf+1ll*ifac[i+a]*ifac[i+b]%p*ifac[j])%p,
                   numg=(numg+1ll*ifac[i+a]*ifac[i+b]%p*ifac[j]%p*i)%p;
    }
    fx[deg]=numf,gx[deg]=numg;
    static int t1[N],t2[N];
    std::fill(t1,t1+deg,0);
    std::fill(t2,t2+deg,0);
    F(L+a,L+b,L+c,deg-2,p-2,t1+1);
    F(a-1,b-1,c-1,deg-2,2,t1+1);
    K(a-1,b-1,c,L,deg-2,2,t1+1);
    J(a,b,c,L,deg-1,p-1,t1);
    K(a,b,c,L,deg-1,1,t2);
    int x(1ll*(p-2)*a%p*b%p);
    for(int i(deg);i;--i)
    {
        fx[i-1]=(2*(p-gx[i])+1ll*(p+i-a-b)*fx[i]+t2[i-1])%p;
        gx[i-1]=(2ll*(p+i-1-a-b)*fx[i-1]+1ll*x*fx[i]+1ll*(p-i-a-b)*gx[i]+2*(i>1?t2[i-2]:0)+t1[i-1])%p*_3%p;
    }
}
void G(int a,int b,int c,int L,int deg,int coef,int res[],int fb[],int gb[],int fc[],int gc[])
{
    assert(a>0);assert(b>0);assert(c>0);if(deg<0) return;
    static int T[N];std::fill(T,T+deg,0);
    FF(a-1,b-1,c-1,deg-1,1,T);FF(a+L,b+L,c+L,deg-1,p-1,T);
    K(a-1,b,c,L,deg-1,p-1,T);
    for(int i(0);i<deg;++i)
    {
        T[i]=(T[i]+1ll*a*(fb[i+1]+fc[i+1])+p-fb[i]+p-fc[i]+gb[i+1]+gc[i+1])%p;
    }
    for(int i(1),t(0);i<=deg;++i)
    {
        t=(3ll*t+T[i-1])*inv[i]%p;
        res[i]=(res[i]+1ll*coef*t)%p;
    }
}
void H(int a,int b,int c,int L,int deg,int coef,int res[])
{
    assert(a>0);assert(b>0);assert(c>0);if(deg<0) return;
    static int fa[N],fb[N],fc[N],ga[N],gb[N],gc[N];
    fg(b,c,a,L,deg-1,fa,ga);
    fg(a,c,b,L,deg-1,fb,gb);
    fg(a,b,c,L,deg-1,fc,gc);
    static int T[N];std::fill(T,T+deg,0);
    G(a,b,c,L,deg-1,p-1,T,fb,gb,fc,gc);
    G(b,a,c,L,deg-1,p-1,T,fa,ga,fc,gc);
    G(c,a,b,L,deg-1,p-1,T,fa,ga,fb,gb);
    res[0]=(res[0]+1ll*(L+1)*coef)%p;
    for(int i(1),t(L+1);i<=deg;++i)
    {
        t=(3ll*t+T[i-1])*inv[i]%p;
        res[i]=(res[i]+1ll*coef*t)%p;
    }
}
int main()
{
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);int a,b,c,n;std::cin>>a>>b>>c>>n;
    fac[0]=1;for(int i(1);i<N;++i) fac[i]=1ll*fac[i-1]*i%p;ifac[N-1]=fp(fac[N-1],p-2);
    for(int i(N-1);i;--i) ifac[i-1]=1ll*ifac[i]*i%p,inv[i]=1ll*ifac[i]*fac[i-1]%p;
    static int fx[N];int x[3]={a,b,c};std::sort(x,x+3);int m(x[2]);
    FF(0,m-x[0],m-x[1],n,p-1,fx);
    H(m-a+1,m-b+1,m-c+1,n-1,n,p-1,fx);
    for(int i(1),j(_3);i<=n;++i,j=1ll*j*_3%p) std::cout<<(m+n+1+1ll*fac[i]*fx[i]%p*j)%p<<"\n";
    return 0;
}