uoj#593 新年的军队 题解
后天南大营,这个趣味编程(注意不是算法)整的人很慌,于是 Delov 在怂恿人写猪国杀。不好评价。
去年写猪国杀的时候我在干嘛来着?哦和 joke3579 加训多项式啊那没事了。他老是说这个题然而没补,现在我补一下。
感觉不如寄希望于微积分和离散数学能拼过一点人。虽然也就是民科水平。线性代数甚至还不如民科水平。
于是决定今天补一道经典老题。话说牛子老师什么时候个人主页重度福瑞控改成轻度了。
致歉:牛子老师没写过重度,我太错。
乐了,看懂题解花了差不多五个小时写的常数还巨大。
题意
对于所有 \(n\) 阶 \(m\) 个降低的排列 \(a\) 求 \(a_k\) 为 \(1-n\) 的个数。
题解
巨大逆天。感觉不光是概率密度的转化还是推式子都是人力所不能及的。EI 的题解确实是人类可以看懂的,写的相当清楚。
首先我们可以把排列映射到 \([0,1]\) 上均匀分布的实数序列,显然每个排列都是等概率出现的。现在我们需要求得 \(k\) 位置的实数分布,因此可以使用概率密度函数扩充状态。
对于 \(n\) 阶排列中排名 \(i\) 的数,显然其概率密度为
设之前的概率密度为 \(f(x)\),那如果添加一个新的数使得大于上一个数,这个位置的概率密度就是 \(\int_0^xf(z)\text dz\)。可以考虑成在这里加一个 \(x\),那么上一个数更小的概率就是这个。反之如果小于,概率显然是 \(\int_x^1f(z)\text dz\)。我们发现我们可以轻松地以如上方法刻画排列的升高或降低。
现在我们列出式子。设 \(F(u,v,z)\) 为有 \(u\) 个 \(<\),\(v\) 个 \(>\),排列结尾处的概率密度。那么考虑后边加一个 \(>\) 还是 \(<\):
两边求导(注意到 \(\int_z^1F(u,v,t)\text dt=C-\int_0^zF(u,v,t)\text dt\)):
显然解为
那么尝试解 \(C\)。代入 \(z=0/1\),设 \(I=\int_0^1F(u,v,t)\text dt\),有方程
那么解出
接下来考虑把两边粘起来。使用 \(x\) 表示符号个数,\(y\) 计算 \(>\)。对于左边,\(<\) 只贡献 \(x\),因此 \(u=x\)。\(>\) 同时贡献 \(x,y\),于是 \(v=xy\)。右边同理,\(u=xy,v=x\)。得到
设 \(n_1=k-1,n_2=n-k\),那么我们需要 \([x^{n_1}]L\) 和 \([x^{n_2}]R\)。接下来有两条路。
有限微积分
因为我不会有限微积分,所以感觉十分奇怪。
提取系数得到
拼起来。
好像是可以直接插值做的,但是常数巨大。不如继续推一下。
里边是个区间求和,可以拆成正方向的无穷和再加上一个差分 \(\Delta_xf(z)=f(z)-f(z+x)\)。设 \(f_s=[y^m](1-y)^{n+1}y^s\),则
差分之后发现常数项被消没了,于是还要补上一个 \(C\)。爆算 \(C\) 大概算不出来,但是观察到 \(C\) 是给每一个数加上相同的值,而所有的总和是 \(\left\langle\begin{matrix}n\\m\end{matrix}\right\rangle\),那么减一下除个 \(n\) 就得到 \(C\)。于是接着带着 \(C\) 推:
泰勒展开告诉我们 \(f(z+x)=e^{x\text D}f(z)\),\(\text D\) 是求导算子。那么设 \(F(x)=\sum_{s\ge 0}f_sx^{s+1}\),把后边大括号单独拆一项出来,有:
另一项类似地处理。于是带回原式:
其中 \(\Sigma\) 为求和算子。然后由于移位算子是 \(e^{\text D}\),于是差分为 \(1-e^{\text D}\),求和即为其逆 \(\dfrac 1{1-e^{\text D}}=\int\dfrac {\text D}{1-e^{\text D}}\),长的类似伯努利数。
这样问题就只剩下算 \(G(x)=F(e^x)\)。这个一会再说。
生成函数
这是另一种方法,对不会有限微积分的我来说似乎更加简单(?)
先不拆系数,把两边合并起来。左边的 \(x\) 携程 \(p\),右边的写成 \(q\),则答案为
分式分解得到:
这个 \(\dfrac 1{e^p-e^q}\) 求不了逆,十分奇怪。奇妙操作一下,求个导看看:
得到了伯努利数的形式。此时对 \(y\) 提取系数就能得到刚才的 \(G(x)\):
那以计算 \(e^{(p-q)z}B(p-q)G(p)\) 为例。换元 \(q=pr\),提取 \([p^{n-1}r^{n_2}]\),得到
前边两项合起来变成 \(Q(p)\):
形式清爽了很多。
最后的问题是怎么算 \(G\)。考虑 \(F\) 的微分方程:设 \(f(x)=\sum_{i=0}^mf_ix^i\),则 \(F(x)=xf(x)\)。又有
则
此处 \(C=(n-m+1)f_0=(n-m+1)\binom{n+1}m(-1)^m\)。
那么 \(G(x)=F(e^x)=e^xf(e^x)=e^xg(x)\),得到
提取系数:
左边只留下 \(G_i\) 可以发现是半在线卷积的形式,那 \(O(n\log^2n)\) 做就好了。可以搞到 \(O(\dfrac{n\log^2n}{\log\log n})\),但是我比较摆所以不搞了。
int n,m,k;
poly f,g;
void solve(int l,int r){
if(l==r){
g[l]=1ll*g[l]*inv[n-l+1]%mod;
return;
}
int mid=(l+r)>>1;
solve(l,mid);
poly f1(mid-l+1),f2(mid-l+1),f3(r-l+1),f4(r-l+1),ans;
for(int i=0;i<=mid-l;i++){
f1[i]=1ll*g[i+l]*(n-m)%mod;
f2[i]=1ll*g[i+l]*(i+l)%mod;
}
for(int i=0;i<=r-l;i++){
f3[i]=(i&1)?(mod-Inv[i]):Inv[i];
f4[i]=(i&1)?Inv[i+1]:(mod-Inv[i+1]);
}
ans=f1*f3+f2*f4;
for(int i=mid+1;i<=r;i++)g[i]=sub(g[i],ans[i-l]);
solve(mid+1,r);
}
int Euler(int n,int m){
int ans=0;
for(int i=0;i<=m;i++)ans=(ans+1ll*(i&1?(mod-1):1)*C(n+1,i)%mod*qpow(m-i+1,n))%mod;
return ans;
}
int main(){
n=read();m=read();k=read();init(n+1<<1);f.resize(n+1);
int n1=k-1,n2=n-k;
poly B(n+1);
for(int i=0;i<=n;i++)B[i]=Inv[i+1];
B=getinv(B);
g.resize(n);g[0]=1ll*(m&1?(mod-1):1)*C(n+1,m)%mod*(n-m+1)%mod;
solve(0,n-1);
poly P(n),Q(n);
for(int i=0;i<n;i++){
P[n-i-1]=1ll*(n1-i&1?(mod-1):1)*g[n-i-1]%mod*C(i,n1)%mod;
Q[n-i-1]=1ll*(n2&1?(mod-1):1)*g[n-i-1]%mod*C(i,n2)%mod;
}
P=(P*B).slice(n);Q=(Q*B).slice(n);
for(int i=0;i<n;i++)f[i]=1ll*(Q[n-i-1]-P[n-i-1]+mod)*Inv[i]%mod;
f=jifen(f);
for(int i=0;i<n;i++)f[i]=1ll*f[i]*jc[n-i-1]%mod;
for(int i=0;i<n;i++)g[i]=Inv[i];
f=(f*g).slice(n);
for(int i=0;i<n;i++)f[i]=1ll*f[i]*jc[i]%mod;
int val=Euler(n,m);
for(int i=0;i<n;i++)val=sub(val,f[i]);
val=1ll*val*inv[n]%mod;
for(int i=0;i<n;i++)f[i]=add(f[i],val);
for(int i=0;i<n;i++)print(f[i]),putchar(' ');puts("");
return 0;
}