题解：P16044 [ICPC 2022 NAC] Double Sort

· · 题解

题意简述

从 1 到 m 中等概率选择 n 个不同整数。将它们排序并差分，再把差分排序后求前缀和。求最终序列每一项的期望。

解题思路

设第一次排序后的数为 a_1<a_2<\dots<a_n，定义：

\begin{aligned} d_1 & = a_1 \\ d_i & = a_i-a_{i-1} \end{aligned}

第二式适用于 2\le i\le n。

所有 d_i 均为正整数，且总和不超过 m。反过来，任意满足这些条件的差分序列都唯一对应一个 n 元子集。因此，原随机过程等价于在所有满足条件的正整数序列中等概率选择一个，总方案数为 \binom{m}{n}。

将差分排序为 y_1\le y_2\le\dots\le y_n。第 k 个答案是 y_1+\dots+y_k 的期望。

固定一个正整数高度 h，记不少于 h 的差分数量为：

N_h=\#\{i:d_i\ge h\}

排序后共有 n-N_h 个数小于 h。因此，前 k 小的数中不少于 h 的恰有：

(N_h-n+k)_+

其中 (x)_+=\max(x,0)。按高度逐层计算每个正整数，便有：

y_1+\dots+y_k=\sum_{h\ge1}(N_h-n+k)_+

接下来求 N_h 的二项矩。任取 r 个下标，并要求这些位置的差分均不少于 h。将选中的差分各减去 h-1 后，仍得到 n 个正整数，但总和上界变为 m-r(h-1)。因此：

\mathbb E\binom{N_h}{r} =\binom{n}{r}\frac{\binom{m-r(h-1)}{n}}{\binom{m}{n}}

现在只需由二项矩恢复截断函数。令 c=n-k。对任意整数 X，先写成尾概率之和：

(X-c)_+=\sum_{q=c+1}^{n}[X\ge q]

二项式反演给出：

[X\ge q]=\sum_{r=q}^{n}(-1)^{r-q}\binom{r-1}{q-1}\binom{X}{r}

交换 q,r 两层求和，并使用：

\sum_{q=c+1}^{r}(-1)^{r-q}\binom{r-1}{q-1} =(-1)^{r-c-1}\binom{r-2}{c-1}

当 c\ge1 时得到：

(X-c)_+ =\sum_{r=c+1}^{n}(-1)^{r-c-1}\binom{r-2}{c-1}\binom{X}{r}

当 c=0 时，截断函数就是 X=\binom{X}{1}。

将 X=N_h 代入，再交换 h,r 的求和顺序。预处理：

T_r=\sum_{q=0}^{\lfloor(m-n)/r\rfloor}\binom{m-rq}{n}

若 c\ge1，第 k 项期望的分子为：

P_k=\sum_{r=c+1}^{n}(-1)^{r-c-1}\binom{r-2}{c-1}\binom{n}{r}T_r

若 c=0，则 P_n=nT_1。最终答案为：

\frac{P_k}{\binom{m}{n}}

式中存在巨大正负项相消，直接用浮点数累加会丢失精度。代码使用 cpp_int 精确计算组合数、正项和与负项和，完成相减后才转为 long double 做一次除法。

组合数 \binom{i}{n} 可由前一项递推。所有 T_r 共包含 O(m\log n+n) 项，全部答案还需 O(n^2) 次运算。因此时间复杂度为 O(m\log n+n^2)，空间复杂度为 O(m+n^2)。

正确性证明

差分序列中的每一项均为正数，其前缀和严格递增且不超过 m，所以它与原来的 n 元子集一一对应。两种随机模型中每个样本的概率相同。

对固定高度 h，排序后的前 k 项中，只有超过前面 n-N_h 个小数的位置才可能达到 h，数量正是 (N_h-n+k)_+。对所有高度求和，恰好还原前 k 个正整数的总和。

计算 \mathbb E\binom{N_h}{r} 时，一边枚举差分序列再选择 r 个达标下标，另一边先选择下标再对它们平移 h-1，两种计数覆盖相同对象。二项式反演恒等式随后逐点还原截断函数，因此没有使用分布近似。

最后交换的都是有限求和，正负项由大整数精确保存。故 P_k/\binom{m}{n} 正是第 k 个前缀和的期望。

参考代码

#include <bits/stdc++.h>
using namespace std;

using ll=long long;
using ull=unsigned long long;
using u128=unsigned __int128;
using ld=long double;
const int N=55;
const int M=10005;
const int L=20;
const ull base=1000000000;
struct Big
{
    ull a[L];
    Big(ull x=0)
    {
        memset(a,0,sizeof a);
        a[0]=x;
    }
    Big &operator+=(const Big &x)
    {
        ull cur=0;
        for(int i=0;i<L;i++)
        {
            cur+=a[i]+x.a[i];
            a[i]=cur%base;
            cur/=base;
        }
        return *this;
    }
    Big &operator-=(const Big &x)
    {
        int cur=0;
        for(int i=0;i<L;i++)
        {
            ll val=(ll)a[i]-x.a[i]-cur;
            if(val<0)
            {
                val+=base;
                cur=1;
            }
            else cur=0;
            a[i]=val;
        }
        return *this;
    }
    void mul(ull x)
    {
        u128 cur=0;
        for(int i=0;i<L;i++)
        {
            cur+=(u128)a[i]*x;
            a[i]=cur%base;
            cur/=base;
        }
    }
    void div(ull x)
    {
        ull cur=0;
        for(int i=L-1;i>=0;i--)
        {
            cur=cur*base+a[i];
            a[i]=cur/x;
            cur%=x;
        }
    }
    ld val()const
    {
        ld res=0;
        for(int i=L-1;i>=0;i--)res=res*base+a[i];
        return res;
    }
};
ull c[N][N];
Big comb[M],sum[N];
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int n,m;
    cin>>n>>m;
    for(int i=0;i<=n;i++)
    {
        c[i][0]=c[i][i]=1;
        for(int j=1;j<i;j++)c[i][j]=c[i-1][j-1]+c[i-1][j];
    }
    comb[n]=Big(1);
    for(int i=n+1;i<=m;i++)
    {
        comb[i]=comb[i-1];
        comb[i].mul(i);
        comb[i].div(i-n);
    }
    for(int i=1;i<=n;i++)
    {
        for(int j=m;j>=n;j-=i)sum[i]+=comb[j];
    }
    cout<<fixed<<setprecision(10);
    for(int i=1;i<=n;i++)
    {
        int k=n-i;
        Big pos,neg;
        if(!k)
        {
            pos=sum[1];
            pos.mul(n);
        }
        else
        {
            for(int j=k+1;j<=n;j++)
            {
                Big cur=sum[j];
                cur.mul(c[j-2][k-1]);
                cur.mul(c[n][j]);
                if((j-k-1)&1)neg+=cur;
                else pos+=cur;
            }
        }
        pos-=neg;
        cout<<pos.val()/comb[m].val()<<'\n';
    }
    return 0;
}