[数学记录]CF923E Perpetual Subtraction

· · 个人记录

填坑。去年省选就看到学长做这题了。

题意 : 有一个[0,n]内的离散随机变量,给出其分布的概率生成函数。

然后进行m轮如下的操作:

设当前这个数为x,将x随机替换为[0,x]中的任意一个数。

问操作后其分布的概率生成函数。

n\leq10^5;\ m\leq 10^{18}

我们设F_r(x)为r轮后的概率生成函数,特殊地,F_0(x)即为开始时给出的幂级数。

我们考虑一次操作之后会有何影响 :

F_{r+1}(x)=\sum\limits_{i=1}^n\sum\limits_{j=i}^n\dfrac{F_r[j]}{j+1}x^i =\sum\limits_{j=1}^n\dfrac{F_r[j]}{j+1}\sum\limits_{i=1}^jx^i =\sum\limits_{j=1}^n\dfrac{F_r[j]}{j+1}\dfrac{x^{j+1}-1}{x-1}

将\frac{1}{x-1}提取出来

=\frac{1}{x-1}\sum\limits_{j=1}^nF_r[j]\dfrac{x^{j+1}-1}{j+1} =\frac{1}{x-1}\sum\limits_{j=1}^n∫_1^xF_r[j]t^j\text{dt} =\frac{1}{x-1}∫_1^xF_r(t)\text{dt}

这里是积分到1令人很不爽,我们尝试使用z=x-1替换x

F(x)=F(z+1)=\frac{1}{z}∫_1^{z+1}F_r(t)\text{dt}

对积分也进行位移:

=\frac{1}{z}∫_0^{z}F_r(t+1)\text{d}(t+1)

现在我们的F_r又很不好受了,方便起见设G_r(t)=F_r(t+1),现在就是不定积分了。

=\frac{1}{z}∫G_r(z) =\frac{1}{z}\sum\limits_{i=0}^nG_r[i]\dfrac{x^{i+1}}{i+1} =\sum\limits_{i=0}^n\dfrac{G_r[i]}{i+1}x^i

我们观察系数,有G_{r+1}[i]=\dfrac{G_r[i]}{i+1},这是个等比数列,随便powM就好了。

问题在于怎么得到G_0(x)?

G(x)=F(x+1)=\sum\limits_{i=0}F[i](x+1)^i =\sum\limits_{i=0}F[i]\sum\limits_{j=0}^i\dbinom{i}{j}x^j =\sum\limits_{j=0}x^i\sum\limits_{i=j}\dbinom{i}{j}F[i] =\sum\limits_{j=0}x^j\sum\limits_{i=j}\dfrac{i!F[i]}{j!(i-j)!} G[k]=\sum\limits_{i=k}\dfrac{i!F[i]}{k!(i-k)!}

设A[m]=(n-m)!F[n-m];B[m]=\frac{1}{m!}

我们还需要逆变换。 $G[k]=\sum\limits_{i=k}\dbinom{i}{k}F[i]\Rightarrow F[k]=\sum\limits_{i=k}(-1)^{k-i}\dbinom{i}{k}G[i]

众所周知二项式反演能够NTT,现在我们能在F,G之间变换了,这道题也就做完了。

#include<algorithm>
#include<cstdio>
#define ll long long
#define mod 998244353
#define G 3
#define Maxn 135000
using namespace std;
inline int read()
{
  register int X=0;
  register char ch=0;
  while(ch<48||ch>57)ch=getchar();
  while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
  return X;
}
ll powM(ll a,int t=mod-2)
{
  ll ans=1;
  while(t){
    if(t&1)ans=ans*a%mod;
    a=a*a%mod;
    t>>=1;
  }return ans;
}
int r[Maxn<<2];
ll invG=powM(G),fac[Maxn],inv[Maxn];
void NTT(ll *f,bool op,int n)
{
  for (int i=0;i<n;i++)
    if (r[i]<i)swap(f[r[i]],f[i]);
  for (int len=1;len<n;len<<=1){
    int w=powM(op==1 ? G:invG,(mod-1)/len/2);
    for (int p=0;p<n;p+=len+len){
      ll buf=1;
      for (int i=p;i<p+len;i++){
        int sav=f[i+len]*buf%mod;
        f[i+len]=f[i]-sav;
        if (f[i+len]<0)f[i+len]+=mod;
        f[i]=f[i]+sav;
        if (f[i]>=mod)f[i]-=mod;
        buf=buf*w%mod;
        }
    }
  }
}
ll g[Maxn<<1];
void times(ll *f,ll *gg,int len,int lim)
{
  int m=len+len,n;
  for(int i=0;i<len;i++)g[i]=gg[i];
  for(n=1;n<m;n<<=1);
  for(int i=len;i<n;i++)g[i]=0;
  for(int i=0;i<n;i++)
    r[i]=(r[i>>1]>>1)|((i&1)?n>>1:0);
  NTT(f,1,n);NTT(g,1,n);
  for(int i=0;i<n;++i)f[i]=f[i]*g[i]%mod;
  NTT(f,0,n);ll invn=powM(n);
  for(int i=0;i<lim;++i)f[i]=f[i]*invn%mod;
  for(int i=lim;i<n;++i)f[i]=0;
}
ll f[Maxn<<1],p[Maxn<<1],m;
int n;
int main()
{
  n=read()+1;
  scanf("%I64d",&m);m%=(mod-1);
  fac[0]=1;
  for (int i=1;i<=n;i++)fac[i]=fac[i-1]*i%mod;
  inv[n]=powM(fac[n]);
  for (int i=n;i;i--)inv[i-1]=inv[i]*i%mod;
  for (int i=0;i<n;i++)f[i]=read();
  for (int i=0;i<n;i++)p[n-i-1]=f[i]*fac[i]%mod;
  for (int i=0;i<n;i++)f[i]=inv[i]%mod;
  times(p,f,n,n);
  for (int i=0;i<n;i++)f[n-i-1]=p[i]*inv[n-i-1]%mod;
  for (int i=0;i<n;i++)f[i]=f[i]*powM(i+1,mod-1-m)%mod;
  for (int i=0;i<n;i++)
    if (i&1)p[n-i-1]=mod-(f[i]*fac[i]%mod);
    else p[n-i-1]=f[i]*fac[i]%mod;
  for (int i=0;i<n;i++)f[i]=inv[i]%mod;
  times(p,f,n,n);
  for (int i=0;i<n;i++)f[n-i-1]=p[i]*inv[n-i-1]%mod;
  for (int i=0;i<n;i++)
    printf("%I64d ",(i&1) ? (mod-f[i])%mod : f[i]);
  puts("");return 0;
}