uoj#593 新年的军队 题解

发布时间 2023-06-01 16:19:17作者: gtm1514

后天南大营,这个趣味编程(注意不是算法)整的人很慌,于是 Delov 在怂恿人写猪国杀。不好评价。

去年写猪国杀的时候我在干嘛来着?哦和 joke3579 加训多项式啊那没事了。他老是说这个题然而没补,现在我补一下。

感觉不如寄希望于微积分和离散数学能拼过一点人。虽然也就是民科水平。线性代数甚至还不如民科水平。

于是决定今天补一道经典老题。话说牛子老师什么时候个人主页重度福瑞控改成轻度了。

乐了,看懂题解花了差不多五个小时写的常数还巨大。

题意

对于所有 \(n\)\(m\) 个降低的排列 \(a\)\(a_k\)\(1-n\) 的个数。

题解

巨大逆天。感觉不光是概率密度的转化还是推式子都是人力所不能及的。EI 的题解确实是人类可以看懂的,写的相当清楚。

首先我们可以把排列映射到 \([0,1]\) 上均匀分布的实数序列,显然每个排列都是等概率出现的。现在我们需要求得 \(k\) 位置的实数分布,因此可以使用概率密度函数扩充状态。

对于 \(n\) 阶排列中排名 \(i\) 的数,显然其概率密度为

\[\frac{x^{j-1}(1-x)^{n-j}}{(j-1)!(n-j)!} \]

设之前的概率密度为 \(f(x)\),那如果添加一个新的数使得大于上一个数,这个位置的概率密度就是 \(\int_0^xf(z)\text dz\)。可以考虑成在这里加一个 \(x\),那么上一个数更小的概率就是这个。反之如果小于,概率显然是 \(\int_x^1f(z)\text dz\)。我们发现我们可以轻松地以如上方法刻画排列的升高或降低。

现在我们列出式子。设 \(F(u,v,z)\) 为有 \(u\)\(<\)\(v\)\(>\),排列结尾处的概率密度。那么考虑后边加一个 \(>\) 还是 \(<\)

\[F(u,v,z)=1+u\int_0^zF(u,v,t)\text dt+v\int_z^1F(u,v,t)\text dt \]

两边求导(注意到 \(\int_z^1F(u,v,t)\text dt=C-\int_0^zF(u,v,t)\text dt\)):

\[F'=(u-v)F \]

显然解为

\[F(u,v,z)=Ce^{(u-v)z} \]

那么尝试解 \(C\)。代入 \(z=0/1\),设 \(I=\int_0^1F(u,v,t)\text dt\),有方程

\[\begin{cases}1+vI=C\\1+uI=Ce^{u-v}\end{cases} \]

那么解出

\[C=\frac{u-v}{u-ve^{u-v}} \]

\[F=\frac{(u-v)e^{(u-v)z}}{u-ve^{u-v}} \]

接下来考虑把两边粘起来。使用 \(x\) 表示符号个数,\(y\) 计算 \(>\)。对于左边,\(<\) 只贡献 \(x\),因此 \(u=x\)\(>\) 同时贡献 \(x,y\),于是 \(v=xy\)。右边同理,\(u=xy,v=x\)。得到

\[L=\frac{(1-y)e^{x(1-y)z}}{1-ye^{x(1-y)}} \]

\[R=\frac{(1-y)e^{x(y-1)(z-1)}}{1-ye^{x(1-y)}} \]

\(n_1=k-1,n_2=n-k\),那么我们需要 \([x^{n_1}]L\)\([x^{n_2}]R\)。接下来有两条路。

有限微积分

因为我不会有限微积分,所以感觉十分奇怪。

提取系数得到

\[\begin{aligned} [x^n]L=&[x^n]\frac{(1-y)e^{x(1-y)z}}{1-ye^{x(1-y)}}\\ =&(1-y)[x^n]\sum_{i\ge 0}y^ie^{x(1-y)i}e^{x(1-y)z}\\ =&(1-y)^{n+1}\sum_{i\ge 0}y^i\frac{(z+i)^n}{n!} \end{aligned} \]

\[\begin{aligned} [x^n]R=&[x^n]\frac{(1-y)e^{x(y-1)(z-1)}}{1-ye^{x(1-y)}}\\ =&(1-y)[x^n]\sum_{i\ge 0}y^ie^{x(1-y)i}e^{x(1-y)(1-z)}\\ =&(1-y)^{n+1}\sum_{i\ge 0}y^i\frac{(1-z+i)^n}{n!} \end{aligned} \]

拼起来。

\[\begin{aligned} &([x^{n_1}]L)([x^{n_2}]R)\\ =&\frac 1{n_1!n_2!}(1-y)^{n+1}\sum_{i\ge 0}y^i(z+i)^{n_1}\sum_{j\ge 0}y^j((j+1)-z)^{n_2} \end{aligned} \]

好像是可以直接插值做的,但是常数巨大。不如继续推一下。

\[\begin{aligned} &[y^m]\frac 1{n_1!n_2!}(1-y)^{n+1}\sum_{i\ge 0}y^i(z+i)^{n_1}\sum_{j\ge 0}y^j((j+1)-z)^{n_2}\\ =&[y^m]\frac {(-1)^{n_2}}{n_1!n_2!}[y^m](1-y)^{n+1}\sum_{i\ge 0}\sum_{j\ge 0}y^{i+j}(z+i)^{n_1}(z-(j+1))^{n_2}\\ =&\frac{(-1)^{n_2}}{n_1!n_2!}\sum_{s=0}^m[y^m](1-y)^{n+1}y^s\sum_{i=0}^s(z+i)^{n_1}(z-(s-i+1))^{n_2} \end{aligned} \]

里边是个区间求和,可以拆成正方向的无穷和再加上一个差分 \(\Delta_xf(z)=f(z)-f(z+x)\)。设 \(f_s=[y^m](1-y)^{n+1}y^s\),则

\[\begin{aligned} =&\frac{(-1)^{n_2}}{n_1!n_2!}\sum_{s=0}^mf_s\Delta_{s+1}\sum_{i\ge 0}(z+i)^{n_1}(z-(s+1))^{n_2}\\ =&\frac{(-1)^{n_2}}{n_1!n_2!}\sum_{s=0}^mf_s\left(\sum_{i\ge 0}\Delta_{s+1}(z+i)^{n_1}(z+i-(s+1))^{n_2}\right)+C\\ \end{aligned} \]

差分之后发现常数项被消没了,于是还要补上一个 \(C\)。爆算 \(C\) 大概算不出来,但是观察到 \(C\) 是给每一个数加上相同的值,而所有的总和是 \(\left\langle\begin{matrix}n\\m\end{matrix}\right\rangle\),那么减一下除个 \(n\) 就得到 \(C\)。于是接着带着 \(C\) 推:

\[\begin{aligned} =&\frac{(-1)^{n_2}}{n_1!n_2!}\sum_{i\ge 0}\left(\sum_{s=0}^mf_s\Delta_{s+1}(z+i)^{n_1}(z+i-(s+1))^{n_2}\right)+C\\ =&\frac{(-1)^{n_2}}{n_1!n_2!}\sum_{i\ge 0}\left(\sum_{s=0}^mf_s((z+i)^{n_1}(z+i-(s+1))^{n_2}-(z+i+(s+1))^{n_1}(z+i)^{n_2})\right)+C \end{aligned} \]

泰勒展开告诉我们 \(f(z+x)=e^{x\text D}f(z)\)\(\text D\) 是求导算子。那么设 \(F(x)=\sum_{s\ge 0}f_sx^{s+1}\),把后边大括号单独拆一项出来,有:

\[\begin{aligned} &\sum_{s=0}^mf_s(z+i)^{n_1}(z+i-(s+1))^{n_2}\\ =&(z+i)^{n_1}\sum_{s=0}^mf_s(z+i-(s+1))^{n_2}\\ =&(z+i)^{n_1}\sum_{s=0}^mf_se^{-(s+1)\text D}(z+i)^{n_2}\\ =&(z+i)^{n_1}F(e^{-\text D}(z+i)^{n_2}) \end{aligned} \]

另一项类似地处理。于是带回原式:

\[\begin{aligned} =&\frac{(-1)^{n_2}}{n_1!n_2!}\sum_{i\ge 0}\left((z+i)^{n_1}F(e^{-\text D})(z+i)^{n_2}-(z+i)^{n_2}F(e^{\text D})(z+i)^{n_1}\right)+C\\ =&\frac{(-1)^{n_2}}{n_1!n_2!}\Sigma(z^{n_1}F(e^{-\text D})z^{n_2}-z^{n_2}F(e^{\text D})z^{n_1})+C \end{aligned} \]

其中 \(\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\),则答案为

\[[y^mp^{n_1}q^{n_2}](1-y)^{n+1}\frac{e^{pz+q(1-z)}}{(1-ye^p)(1-ye^q)} \]

分式分解得到:

\[(1-y)^{n+1}e^{pz+q(1-z)}\frac 1{e^p-e^q}\left(\frac{e^p}{1-ye^p}-\frac{e^q}{1-ye^q}\right) \]

这个 \(\dfrac 1{e^p-e^q}\) 求不了逆,十分奇怪。奇妙操作一下,求个导看看:

\[\begin{aligned} &(1-y)^{n+1}e^{pz+q(1-z)}\frac{p-q}{e^p-e^q}\left(\frac{e^p}{1-ye^p}-\frac{e^q}{1-ye^q}\right)\\ =&(1-y)^{n+1}e^{(p-q)z+q}\frac{(p-q)e^{-q}}{e^{p-q}-1}\left(\frac{e^p}{1-ye^p}-\frac{e^q}{1-ye^q}\right)\\ =&(1-y)^{n+1}e^{(p-q)z}\frac{p-q}{e^{p-q}+1}\left(\frac{e^p}{1-ye^p}-\frac{e^q}{1-ye^q}\right)\\ =&(1-y)^{n+1}e^{(p-q)z}B(p-q)\left(\frac{e^p}{1-ye^p}-\frac{e^q}{1-ye^q}\right)\\ \end{aligned} \]

得到了伯努利数的形式。此时对 \(y\) 提取系数就能得到刚才的 \(G(x)\)

\[e^{(p-q)z}B(p-q)(G(p)-G(q)) \]

那以计算 \(e^{(p-q)z}B(p-q)G(p)\) 为例。换元 \(q=pr\),提取 \([p^{n-1}r^{n_2}]\),得到

\[\sum_{i\ge 0}([p^{n-i-1}]G(p))([r^{n_2}](1-r)^i)[p^j]B(p)e^{pz} \]

前边两项合起来变成 \(Q(p)\)

\[\begin{aligned} =&\sum_{i\ge 0}Q_{n-i-1}[p^j]B(p)e^{pz}\\ =&[p^{n-1}]Q(p)B(p)e^{pz}\\ =&\sum_{i\ge 0}[p^{n-i-1}]Q(p)B(p)\frac{z^i}{i!} \end{aligned} \]

形式清爽了很多。

最后的问题是怎么算 \(G\)。考虑 \(F\) 的微分方程:设 \(f(x)=\sum_{i=0}^mf_ix^i\),则 \(F(x)=xf(x)\)。又有

\[f(x)=\sum_{i=0}^m\binom{n+1}{m-i}(-1)^{m-i}x^i \]

\[(n-m+1+i)f_i=-(m-i+1)f_{i-1}+[i==0]C \]

\[(n-m+1)f(x)+xf'(x)=-mxf(x)+x^2f'(x)+c \]

此处 \(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)\),得到

\[(n-m+me^x+1)g(x)+(1-e^x)g'(x)=C \]

\[(n-m+me^x+e^x)g(x)+(1-e^x)(g(x)+g'(x))=C \]

\[((n-m)e^{-x}+m+1)G(x)+(e^{-x}-1)G'(x)=C \]

提取系数:

\[(n-m)\sum_{j=0}^i\frac{(-1)^j}{j!}G_{i-j}+(m+1)G_i+\sum_{j=1}^i\frac{(-1)^j}{j!}(i-j+1)G_{i-j+1}=[i==0]C \]

\[(n-m)\sum_{j=1}^i\frac{(-1)^j}{j!}G_{i-j}+(n+1-i)G_i+\sum_{j=1}^{i-1}\frac{(-1)^j}{j!}(i-j+1)G_{i-j+1}=[i==0]C \]

左边只留下 \(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;
}