CF923E Perpetual Subtraction

Link
这题有两种推法，一种是生成函数，还有一种是线性代数。

一、生成函数

设进行\(t\)轮之后得到的数为\(i\)的概率为\(f_{t,i}\)。
我们知道有：

\[f_{0,i}=p_i,f_{t+1,i}=\sum\limits_{j=i}^n\frac{f_{t,j}}{j+1} \]

设\(P_t(x)=\sum\limits_{i=0}^nf_{t,i}x^i\)，将上式代入得到：

\[\begin{aligned} P_{t+1}(x)&=\sum\limits_{i=0}^n \sum\limits_{j=i}^n\frac{f_{t,j}}{j+1}x^i\\ &=\sum\limits_{j=0}^n\frac{f_{t,j}}{j+1}\sum\limits_{i=0}^jx^i\\ &=\sum\limits_{j=0}^n\frac{f_{t,j}}{j+1}\frac{x^{j+1}-1}{x-1}\\ &=\frac1{x-1}\sum\limits_{j=0}^nf_{t,j}\frac{x^{j+1}-1}{j+1}\\ &=\frac1{x-1}\sum\limits_{j=0}^nf_{t,j}\int_1^xy^j\mathrm dy\\ &=\frac1{x-1}\int_1^xP_t(y)\mathrm dy \end{aligned} \]

我们知道\(\int_0^xf(y)\mathrm dy=\int f(x)\mathrm dx\)，因此设\(P_t(x+1)=Q_t(x)=\sum\limits_{i=0}^ng_{t,i}x^i\)，替换得到：

\[Q_{t+1}(x)=\frac1x\int_0^x Q_t(y)\mathrm dy=\sum\limits_{i=0}^n\frac{g_{t,i}}{i+1}x^i \]

也就是说：

\[g_{t+1,i}=\frac{g_{t,i}}{i+1},g_{m,i}=\frac{g_{0,i}}{(i+1)^m} \]

利用二项式定理得到：

\[g_{t,i}=\sum\limits_{j=i}^n{j\choose i}f_{t,j}\tag1 \]

代入\(t=0,f_{0,i}=p_i\)得到：

\[g_{0,i}=\sum\limits_{j=i}^n{j\choose i}p_j \]

先利用卷积求出\(g_{0,i}\)，再利用\(g_{m,i}=\frac{g_{0,i}}{(i+1)^m}\)即可求出\(g_{m,i}\)。
然后对\((1)\)式进行二项式反演：

\[f_{t,i}=\sum\limits_{j=i}^n(-1)^{j-i}{j\choose i}g_{t,j} \]

最后卷积求出\(f_{m,i}\)即可。

二、线性代数

转移矩阵的特征值是显然的，通过人类智慧可以得到转移矩阵的特征向量。
然后利用NTT优化各种矩阵乘法的过程即可。

#include<cstdio>
#include<cctype>
#include<cstring>
#include<algorithm>
using i64=long long;
const int N=262144,P=998244353;
char ibuf[(1<<20)+1],*iS=ibuf;
int lim,rev[N],w[N],p[N],fac[N],ifac[N],f[N],g[N];
i64 read(){i64 x=0;while(isspace(*iS))++iS;while(isdigit(*iS))(x*=10)+=*iS++&15;return x;}
int inc(int a,int b){return a+=b-P,a+(a>>31&P);}
int dec(int a,int b){return a-=b,a+(a>>31&P);}
int mul(int a,int b){return 1ll*a*b%P;}
int pow(int a,int k){int r=1;for(;k;k>>=1,a=mul(a,a))if(k&1)r=mul(a,r);return r;}
void init(int n)
{
    int p=(lim=1<<(33-__builtin_clz(n)))>>1,g=pow(3,(P-1)/lim);
    w[p]=1,fac[0]=ifac[0]=fac[1]=ifac[1]=1;
    for(int i=1;i<lim;++i) rev[i]=(rev[i>>1]>>1)|(i&1? p:0);
    for(int i=p+1;i<lim;++i) w[i]=mul(w[i-1],g);
    for(int i=p-1;i;--i) w[i]=w[i<<1];
    for(int i=1;i<=n;++i) fac[i]=mul(fac[i-1],i);
    ifac[n]=pow(fac[n],P-2);
    for(int i=n;i;--i) ifac[i-1]=mul(ifac[i],i);
}
void NTT(int*a,int f)
{
    if(!~f) std::reverse(a+1,a+lim);
    for(int i=1;i<lim;++i) if(i<rev[i]) std::swap(a[i],a[rev[i]]);
    for(int i=1;i<lim;i<<=1) for(int j=0,d=i<<1;j<lim;j+=d) for(int k=0,x;k<i;++k) x=mul(a[i+j+k],w[i+k]),a[i+j+k]=dec(a[j+k],x),a[j+k]=inc(a[j+k],x);
    if(!~f) for(int i=0,x=P-(P-1)/lim;i<lim;++i) a[i]=mul(a[i],x);
}
int main()
{
    fread(ibuf,1,1<<20,stdin);
    int n=read()+1,m=read()%(P-1);init(n-1);
    for(int i=0;i<n;++i) p[i]=read();
    for(int i=0;i<n;++i) f[i]=ifac[i],g[i]=mul(fac[n-i-1],p[n-i-1]);
    NTT(f,1),NTT(g,1);
    for(int i=0;i<lim;++i) f[i]=mul(f[i],g[i]);
    NTT(f,-1);
    for(int i=0;i<n;++i) g[i]=mul(f[i],pow(pow(n-i,m),P-2)),f[i]=i&1? P-ifac[i]:ifac[i];
    memset(f+n,0,(lim-n)<<2),memset(g+n,0,(lim-n)<<2),NTT(f,1),NTT(g,1);
    for(int i=0;i<lim;++i) f[i]=mul(f[i],g[i]);
    NTT(f,-1);
    for(int i=0;i<n;++i) printf("%d ",mul(ifac[i],f[n-i-1]));
}