题意

给定分数 $r$ 和正整数 $d$（$0\le r\lt 1,d\le 10^7$），求 $\displaystyle\sum_{i\ge 0}r^ii^d\pmod{998244353}$

我的做法

首先用斯特林数把后面那个拆开

$$ \begin{aligned} &\sum_{i\ge 0}r^ii^d\\ =&\sum_{i\ge 0}r^i\sum_{j}\left\{\begin{matrix}d\\j\end{matrix}\right\}i^{\underline{j}}\\ =&\sum_{j\ge 0}\left\{\begin{matrix}d\\j\end{matrix}\right\}\sum_{i\ge 0}\binom{i}{j}j!r^i\\ =&\sum_{j\ge 0}\left\{\begin{matrix}d\\j\end{matrix}\right\}j!\dfrac{r^{j}}{\left(1-r\right)^{j+1}}\\ =&\sum_{j= 0}^d\sum_{t=0}^j\left(-1\right)^{j-t}\binom{j}{t}t^d\dfrac{r^j}{\left(1-r\right)^{j+1}}\\ =&\dfrac{1}{1-r}\sum_{t=0}^dt^d\left(\dfrac{r}{1-r}\right)^t\sum_{j=t}^d\binom{j}{t}\left(\dfrac{r}{r-1}\right)^{j-t} \end{aligned} $$

记 $\dfrac{r}{r-1}=k$，那么后面一项求的是 $\displaystyle f\left(m,p\right)=\sum_{i=0}^{m}\binom{i+p}{i}k^i$，有如下递推式

$$f\left(m,p\right)=\begin{cases}f\left(m-1,p\right)\cdot k+f\left(m,p-1\right) & p\ge 1\\f\left(m-1,p\right)\cdot k+1 & p=0\end{cases}$$

记 $\displaystyle F_{q}\left(z\right)=\sum_{m+p=q}f\left(m,p\right)z^p$

可以得到 $\displaystyle F_q\left(z\right)=F_{q-1}\left(z\right)\cdot \left(z+k\right)+1$

因此 $\displaystyle F_q\left(z\right)=\dfrac{\left(z+k\right)^{q+1}-1}{\left(z+k\right)-1}$

可以线性算出所有系数，即原问题中的 $f\left(d,0\right),f\left(d-1,1\right),\cdots f\left(d-t,t\right)\cdots f\left(0,d\right)$

线性筛出所有的 $t^d$ 即可在 $\mathcal{O}\left(d\right)$ 时间内解决原问题

Code

#include <cstdio>
#include <algorithm>
using namespace std;
constexpr int N=10000010,p=998244353;
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;}
int r,d,inv[N],f[N],ps[N],pc,pd[N];bool vis[N];
int main(){
    scanf("%d%d",&r,&d);
    if(r==0) return printf("%d\n",d==0),0;
    inv[1]=1;
    for(int i=2;i<=d+1;++i) inv[i]=1ll*(p-p/i)*inv[p%i]%p;
    int k=1ll*r*fp(r-1,p-2)%p;
    for(int i=0,j=fp(k,d+1),ik=fp(k,p-2);i<=d+1;++i,j=1ll*ik*j%p*(d+2-i)%p*inv[i]%p) f[i]=j;
    k=(p+1-k)%p;
    for(int i=d+1;i;--i) f[i-1]=(f[i-1]+1ll*f[i]*k)%p;
    pd[0]=(d==0);
    pd[1]=1;
    for(int i=2;i<=d;++i){
        if(!vis[i]) pd[i]=fp(i,d),ps[++pc]=i;
        for(int j=1;j<=pc && i*ps[j]<=d;++j){
            vis[i*ps[j]]=true;
            pd[i*ps[j]]=1ll*pd[i]*pd[ps[j]]%p;
            if((i%ps[j])==0) break;
        }
    }
    int ans=0;k=1ll*r*fp(p+1-r,p-2)%p;
    for(int i=0,j=1;i<=d;++i,j=1ll*j*k%p){
        ans=(ans+1ll*pd[i]*f[i+1]%p*j)%p;
    }
    ans=1ll*ans*fp(p+1-r,p-2)%p;
    printf("%d\n",ans);
    return 0;
}

UPD: Binomial Sum 做法

求 $\;\;\displaystyle \sum_{i=0}^{\infty}q^ii^n$

令 $F\left(x\right)=\frac{1}{1-qx}$

满足 $qF+\left(qx-1\right)F^{\prime}=0$

令 $\mathscr{F}\left(x+1\right)$ 为 $F\left(x+1\right)$ 的 $\bmod\;x^{n+1}$ 截断

那么所求即为 $\mathscr{F}\left(x\right)$ 的 $0,1,\dots,n$ 次项

而 $qF\left(x+1\right)+\left(qx+q-1\right)F^{\prime}\left(x+1\right)=0$ 中受到干扰的项为 $qf_n+qnf_n+\left(q-1\right)\left(n+1\right)f_{n+1}=0$ （其中 $f_n=\left[x^n\right]F\left(x+1\right)$ ）

$$ \begin{aligned} q\mathscr{F}\left(x+1\right)+\left(qx+q-1\right)\mathscr{F}^{\left(1\right)}\left(x+1\right)&=-f_{n+1}\left(n+1\right)\left(q-1\right)x^n \\ q\mathscr{F}\left(x\right)+\left(qx-1\right)\mathscr{F}^{\left(1\right)}\left(x\right)&=-f_{n+1}\left(n+1\right)\left(q-1\right)\left(x-1\right)^n \end{aligned} $$

于是可得递推式

$$ qf^{\prime}_k+qkf^{\prime}_{k}-\left(k+1\right)f^{\prime}_{k+1}=-f_{n+1}\left(n+1\right)\left(q-1\right)\left(-1\right)^{n-k}\binom{n}{k} $$

其中 $f_{n+1}=\dfrac{q^{n+1}}{\left(1-q\right)^{n+2}}$

Code

#include <algorithm>
#include <iostream>
constexpr int N(1e7+10),p(998244353);
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;}
int main()
{
    static int fac[N],ifac[N],inv[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)-1);i>=0;--i) ifac[i]=1ll*(i+1)*ifac[i+1]%p,inv[i+1]=1ll*ifac[i+1]*fac[i]%p;
    int r,d;std::cin>>r>>d;
    if(!r)
    {
        if(d) std::cout<<0<<std::endl;
        else std::cout<<1<<std::endl;
        return 0;
    }
    int fnp1(1ll*fp(r,d+1)*fp(fp(p+1-r,p-2),d+2)%p);
    int val(1ll*(d+1)*fnp1%p*(p+1-r)%p);
    int invr(fp(r,p-2));
    static int res[N];
    for(int i(d),rval(val);i>=0;--i,rval=1ll*rval*(p-i-1)%p*inv[d-i]%p)
    {
        int cur((1ll*(i+1)*res[i+1]+rval)%p);
        res[i]=1ll*invr*inv[i+1]%p*cur%p;
        // std::cerr<<res[i]<<" ";
    }
    // std::cerr<<std::endl;
    static int pw[N];
    {
        static int ps[N],pc;static bool vis[N];
        pw[1]=1;pw[0]=fp(0,d);
        for(int i(2);i<=d;++i)
        {
            if(!vis[i])
            {
                pw[i]=fp(i,d);
                ps[++pc]=i;
            }
            for(int j(1);j<=pc && ps[j]*i<=d;++i)
            {
                pw[ps[j]*i]=1ll*pw[ps[j]]*pw[i]%p;
                vis[ps[j]*i]=true;
                if((i%ps[j])==0) break;
            }
        }
        // for(int i(0);i<=d;++i) std::cerr<<pw[i]<<" ";
        // std::cerr<<std::endl;
    }
    int ans(0);
    for(int i(0);i<=d;++i) ans=(ans+1ll*res[i]*pw[i])%p;
    std::cout<<ans<<std::endl;
    return 0;
}

推广

求 $\;\;\displaystyle \sum_{i=0}^{n-1}q^ii^m$

令 $F\left(x\right)=\frac{1}{1-qx}\bmod x^n$，则 $F$ 满足 $qF+\left(qx-1\right)F^{\prime}=nq^nx^{n-1}$

令 $\mathscr{F}\left(x+1\right)$ 为 $F\left(x+1\right)$ 的 $\bmod\;x^{m+1}$ 截断

那么 $qF\left(x+1\right)+\left(qx+q-1\right)F^{\prime}=nq^n\left(x+1\right)^{n-1}$ 中收到干扰的项为 $qf_m+qmf_m+\left(q-1\right)\left(m+1\right)f_{m+1}=nq^n\binom{n-1}{m}$

那么

$$ \begin{aligned} q\mathscr{F}\left(x+1\right)+\left(qx+q-1\right)\mathscr{F}^{\left(1\right)}\left(x+1\right)&=\left(nq^n\left(x+1\right)^{n-1}\bmod x^{m+1}\right)-f_{m+1}\left(m+1\right)\left(q-1\right)x^m \\ q\mathscr{F}\left(x\right)+\left(qx-1\right)\mathscr{F}^{\left(1\right)}\left(x\right)&=nq^n\left(\sum_{k=0}^m h\left(n-1,m,k\right)x^k\right)-f_{m+1}\left(m+1\right)\left(q-1\right)\left(x-1\right)^m \end{aligned} $$

其中

$$ \begin{aligned} h\left(n,m,k\right)&=\sum_{i\le m}\binom{n}{i}\binom{i}{k}\left(-1\right)^{i-k}\\ &=\sum_{i}\binom{n}{k}\binom{n-k}{i-k}\left(-1\right)^{i-k} \\ &=\binom{n}{k}\binom{n-k-1}{m-k}\left(-1\right)^{m-k} \\ &=\dfrac{\left(-1\right)^{m-k}n^{\underline{k}}\left(n-m\right)^{\overline{m-k}}}{\left(m-k\right)!k!} \end{aligned} $$

考虑如何计算 $f_{m+1}$，易知 $f_{0}=\sum_{i=0}^{n-1}q^i=\frac{q^n-1}{q-1}$，那么直接用之前的递推式递推即可。

Code

#include <algorithm>
#include <iostream>
constexpr int N(1e7+10),p(998244353);
int fp(int a,long long 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;}
int redu(int k){return k>=p?k-p:k;}
int main()
{
    static int fac[N],ifac[N],inv[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)-1);i>=0;--i) ifac[i]=1ll*(i+1)*ifac[i+1]%p,inv[i+1]=1ll*ifac[i+1]*fac[i]%p;
    int r,d;long long n;std::cin>>r>>d>>n;
    if(!n)
    {
        std::cout<<0<<std::endl;
        return 0;
    }
    static int pw[N];
    {
        static int ps[N],pc;static bool vis[N];
        pw[1]=1;pw[0]=fp(0,d);
        for(int i(2);i<=d+1;++i)
        {
            if(!vis[i])
            {
                pw[i]=fp(i,d);
                ps[++pc]=i;
            }
            for(int j(1);j<=pc && ps[j]*i<=d+1;++i)
            {
                pw[ps[j]*i]=1ll*pw[ps[j]]*pw[i]%p;
                vis[ps[j]*i]=true;
                if((i%ps[j])==0) break;
            }
        }
    }
    if(!r)
    {
        if(d) std::cout<<0<<std::endl;
        else std::cout<<1<<std::endl;
        return 0;
    }
    int val;
    if(n>d+1)
    {
        int f0(1ll*(p+fp(r,n)-1)*fp(r-1,p-2)%p);
        int np((n-1)%p),modn(n%p),ivr1(fp(p+1-r,p-2));
        for(int i(0),rval(1ll*(p-modn)*fp(r,n)%p);i<=d;++i,rval=1ll*rval*inv[i]%p*(p+np-i+1)%p)
        {
            int lhs((1ll*r*(i+1)%p*f0+rval)%p);
            f0=1ll*inv[i+1]*ivr1%p*lhs%p;
        }
        val=1ll*(p-f0)*(d+1)%p*(r-1)%p;
    }
    else
    {
        val=0;
    }
    int invr(fp(r,p-2));
    static int res[N],bias[N];
    for(int i(d),rval(val);i>=0;--i,rval=1ll*rval*(p-i-1)%p*inv[d-i]%p)
    {
        bias[i]=rval;
    }
    if(n-1<=d)
    {
        bias[n-1]=(bias[n-1]+n*fp(r,n))%p;
    }
    else
    {
        static int fn[N];
        int np((n-1)%p),nm((n-1-d)%p),modn(n%p);
        fn[0]=1;
        for(int i(1);i<=d;++i) fn[i]=1ll*fn[i-1]*(p+np-i+1)%p;
        for(int i(0),rval(1ll*ifac[d]*modn%p*fp(r,n)%p);
            i<=d;++i,rval=1ll*rval*(d-i+1)%p*(p-inv[i])%p*(nm+i-1)%p)
        {
            bias[d-i]=(bias[d-i]+1ll*rval*fn[d-i])%p;
        }
    }
    for(int i(d);i>=0;--i)
    {
        int cur((1ll*(i+1)*res[i+1]+bias[i])%p);
        res[i]=1ll*invr*inv[i+1]%p*cur%p;
    }
    int ans(0);
    for(int i(0);i<=d;++i) ans=(ans+1ll*res[i]*pw[i])%p;
    std::cout<<ans<<std::endl;
    return 0;
}