题解：P16044 [ICPC 2022 NAC] Double Sort

· · 题解

题意简述

从 1\sim m 中等概率选择 n 个不同整数。将其排序、差分、再次排序并求前缀和。求最终序列每个位置的期望值。

解题思路

设排序后的原序列为 a_1<a_2<\dots<a_n。定义 x_1=a_1，并令 x_i=a_i-a_{i-1}（2\le i\le n）。

每个 x_i 都是正整数，且总和不超过 m。反过来，这样的差分序列唯一确定一个子集。因此共有 \binom{m}{n} 种等概率的差分序列。

将差分排序为 y_1\le\dots\le y_n。第 k 个答案是前 k 项之和。对正整数按高度拆分，记 N_h 为不少于 h 的差分个数，则：

\sum_{i=1}^k y_i=\sum_{h\ge1}(N_h-n+k)_+

其中 (z)_+=\max(z,0)。固定 r 个差分，并要求它们都不少于 h。每个指定差分减去 h-1 后，方案数为 \binom{m-r(h-1)}{n}。交换枚举差分序列和下标集合的顺序，可得：

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

令 c=n-k。当 c\ge1 时，Newton 展开给出：

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

预处理：

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

交换 h 与 r 的求和顺序后，c\ge1 时的答案为：

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

当 c=0 时，答案为 nT_1/\binom{m}{n}。

交错求和会严重损失浮点精度。代码用大整数精确计算正项和与负项和，仅在最后做一次除法。组合数按相邻两项递推。总计需要 O(m\log n+n^2) 次大整数运算。

参考代码

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

using ll=long long;
using ull=unsigned long long;
using ui=unsigned int;
using u128=unsigned __int128;
using ld=long double;
const int N=55;
const int M=10005;
const ull base=1000000000;
struct big
{
    vector<ui> a;
    big(ull x=0)
    {
        while(x)
        {
            a.push_back(x%base);
            x/=base;
        }
    }
    void trim()
    {
        while(!a.empty()&&a.back()==0)a.pop_back();
    }
    big &operator+=(const big &x)
    {
        int n=a.size(),m=x.a.size();
        if(n<m)
        {
            a.resize(m);
            n=m;
        }
        ull cur=0;
        for(int i=0;i<n;i++)
        {
            cur+=a[i]+(i<m?x.a[i]:0);
            a[i]=cur%base;
            cur/=base;
        }
        if(cur)a.push_back(cur);
        return *this;
    }
    big &operator-=(const big &x)
    {
        int n=a.size(),m=x.a.size();
        ll cur=0;
        for(int i=0;i<n;i++)
        {
            cur=(ll)a[i]-(i<m?x.a[i]:0)-cur;
            if(cur<0)
            {
                cur+=base;
                a[i]=cur;
                cur=1;
            }
            else
            {
                a[i]=cur;
                cur=0;
            }
        }
        trim();
        return *this;
    }
    void mul(ull x)
    {
        u128 cur=0;
        for(auto &v:a)
        {
            cur+=(u128)v*x;
            v=cur%base;
            cur/=base;
        }
        while(cur)
        {
            a.push_back(cur%base);
            cur/=base;
        }
    }
    void div(ui x)
    {
        int n=a.size();
        ull cur=0;
        for(int i=n-1;i>=0;i--)
        {
            cur=cur*base+a[i];
            a[i]=cur/x;
            cur%=x;
        }
        trim();
    }
    ld val()const
    {
        int n=a.size();
        ld res=0;
        for(int i=n-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==0)
        {
            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;
}