Luogu5907

前言

我去写了一道叫做礼物的数数题，三分钟推完柿子直接做，然后 Peter 说洛谷上有个加强版。

我感觉很 Hard，所以过来学了个科技。

然后又不会插值的想法，所以去学了学传说中的 EI 科技 Binomial Sums。

---

思路

我们令 \(n\leftarrow n+1\)。

答案即 \(\sum_0^nx^ka^x\delta x\)。

\[\sum x^ka^x\delta x= x^k{a^x\over a-1} -{a\over a-1}\sum a^x\Delta x^k\delta x\]

\[\sum_0^nx^ka^x\delta x={n^ka^n\over a-1}+{a\over1-a}\sum_{p=0}^{k-1}\binom kp\sum_0^na^xx^p\delta x \]

\[f_k=\sum_0^nx^ka^x\delta x \]

算出 \(f\) 的 EGF，显然答案即

\[[{z^k\over k!}]{1-(ae^z)^n\over1-ae^z} \]

我们记

\[F(z)={1-z^n\over1-z} \]

答案即 \([{z^k\over k!}]F(ae^z)\)。

记

\[\mathcal F(z+a)=F(z+a)\bmod z^{k+1} \]

答案即 \([{z^k\over k!}]\mathcal F(ae^z)\)。

既有

\[(1-z)F'(z)=-nz^{n-1}+F(z) \]

于是

\[(1-z-a)F'(z+a)=-n(z+a)^{n-1}+F(z+a) \]

可得

\[(1-z-a)\mathcal F'(z+a)-\mathcal F(z+a)=-((z+a)^n\bmod z^{k+1})'-(k+1)z^k[z^k]\mathcal F(z+a) \]

\[(1-z)\mathcal F'(z)-\mathcal F(z)=-(((z+a)^n\bmod z^{k+1})\circ(z-a))'-(k+1)(z-a)^k[z^k]\mathcal F(z) \]

而

\[\sum_{i=0}^k\binom nia^{n-i}(z-a)^i=\sum_{p=0}^k\binom npz^pa^{n-p}(-1)^{k-p}\binom{n-p-1}{k-p} \]

不难想到提取系数转成递推。

---

然后 \([z^0]\mathcal F\) 咋算？

\[[z^0]\mathcal F=\sum_{i=0}^k(-a)^i\sum_{j=i}^{n-1}\binom jia^{j-i}=1+(-1)^ka\sum_{j=0}^{n-2}a^j\binom jk \]

\({\rm I}.a\neq1\)

\[g_k=\sum_0^{n-1}a^x\binom xk\delta x \]

\[\sum a^x\binom xk\delta x ={a^x\binom xk\over a-1}-{a\over a-1}\sum a^x\binom x{k-1}\delta x \]

\[g_k={a^{n-1}\binom{n-1}k\over a-1}-{a\over a-1}g_{k-1} \]

\(g\) 很容易线性递推。

\({\rm II}.a=1\)

\[[z^0]\mathcal F=1+(-1)^k\binom{n-1}{k+1} \]

---

记 \(w=[z^k]\mathcal F(z)\)，那么提取系数转成递推，每个数可以表示为 \(kw+b\) 形式，最后高斯消元解出 \(w\) 既得 \(\mathcal F\) 系数，于是答案就可做了。（比较奇怪的是没有出现解出无数多个 \(w\) 的情况，个中理由我也不大清楚）

\[h_p=[z^p]\mathcal F \]

\[ans=[{z^k\over k!}]\mathcal F(ae^z) =\sum_ih_ia^i[{z^k\over k!}]e^{iz} =\sum_ih_ia^ii^k \]

线性筛预处理 \(i^k\)，总复杂度 \(O(k+\log n)\)。

---

Code

代码挂一下吧，稍微卡了卡空间常数。

const ullt Mod=998244353,g=3;
typedef mod_ullt<Mod>modint;
typedef std::vector<modint>modvec;
typedef poly_NTT<Mod,g>poly;
typedef poly_eval<Mod,g>eval;
typedef poly_inter<Mod,g>inter;
typedef poly_cpd<Mod,g>cpd;
modint A[10000005],B[10000005];
modint F1[10000005],F2[10000005];
bol Ok[10000005];
int main()
{
    A[0]=1;for(uint i=1;i<=10000001;i++)A[i]=A[i-1]*i;
    B[10000001]=A[10000001].inv();for(uint i=10000001;i;i--)B[i-1]=B[i]*i;
    uint n,k;modint a;
    scanf("%u",&n),n++,a.read(),scanf("%u",&k);
    modint t1=a._power(n),t2=(-a)._power(k),q=a.inv(),w;
    modint i1=(a-1).inv(),i2=a._power(n-1);
    if(a==1)
    {
        F1[0]=B[k+1];
        for(uint i=0;i<=k;i++)F1[0]*=n-1-i;
        F1[0]=1+((k&1)?Mod-1:1)*F1[0];
    }
    else{
        modint now=(i2-1)*i1;w=n-1;
        for(uint i=1;i<=k;i++)
        {
            now=(-a*now+w*i2)*i1,w=w*(n-i-1)*B[i+1]*A[i];
        }
        F1[0]=now*a*(k&1?Mod-1:1)+1;
    }
    F2[k]=1;for(uint i=k-1;i;i--)F2[i]=F2[i+1]*(n-i-1)*B[k-i]*A[k-i-1];
    w=1;
    for(uint i=1;i<=k;i++)
    {
        t1*=q,w*=(n-i+1)*B[i]*A[i-1];
        F1[i]=F1[i-1]-w*t1*(((k-i)&1)?Mod-1:1)*F2[i],F2[i]=F2[i-1]-(k+1)*A[k]*B[i-1]*B[k-i+1]*t2*B[i]*A[i-1];
        t2*=-q;
    }
    w=F1[k]/(1-F2[k]);
    for(uint i=0;i<=k;i++)F2[i]=F1[i]+F2[i]*w;
    modint ans;
    F1[0]=0,F1[1]=1;
    std::vector<uint>Prime;
    for(uint i=2;i<=k;i++)
    {
        if(!Ok[i])
        {
            Prime.push_back(i);
            w=modint(i)._power(k);
            for(ullt j=i;j<=k;j*=i)
            {
                Ok[j]=true;
                F1[j]=F1[j/i]*w;
            }
        }
        for(ullt p:Prime)
            if(i*p<=k&&i%p)
            {
                for(ullt t=p;t*i<=k;t*=p)
                {
                    Ok[t*i]=true;
                    F1[t*i]=F1[t]*F1[i];
                }
            }
            else break;
    }
    w=1;
    for(uint i=0;i<=k;i++)
        ans+=F2[i]*w*F1[i],w*=a;
    ans.println();
    return 0;
}