[数学记录]CF923E Perpetual Subtraction
command_block
2020-02-12 22:55:53
填坑。去年省选就看到学长做这题了。
**题意** : 有一个$[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!}$
$=\sum\limits_{i=k}A[n-i]B[i-k]$,那么$A*B$取后$n$项再反过来即可。
我们还需要逆变换。
$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$之间变换了,这道题也就做完了。
```cpp
#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;
}
```