跳转至

F - Social Distance

题意

長さ N+1 の整数列 X_0,X_1,\ldots,X_N が与えられます。ここで、0=X_0<X_1<\ldots<X_N です。

これから、1 から N までの番号のついた N 人の人が、数直線上に現れます。人 i は、区間 [X_{i−1},X_i] の中から一様ランダムに選ばれた実数座標に出現します。

人と人の距離の最小値の期待値を \bmod 998244353 で求めてください。

解析

考虑求出人与人距离的最小值的互补累积分布函数 f(z) (即 f(z) 表示人与人的距离的最小值大于等于 z 的概率),这样我们就可以得到答案为 \displaystyle\int_{0}^{+\infty}f(z)\mathrm{d}z

先尝试对固定的 z 求出 f(z)。设第 i 个人的位置为 x_i,则对于固定的 z,我们需要对 i=1\ldots n-1 满足 x_i+z\le x_{i+1}。设 y_i=x_i-iz,则我们的要求可以改写为 y_1\le y_2\le\ldots\le y_n。于是我们只需求出在 n 个区间均匀随机分布的变量 y_i\in[X_{i-1}-iz,X_i-iz] 满足 y_1\le y_2\le\ldots\le y_n 的概率。

这可以用动态规划 (严格来讲,这不是动态规划,而是递推) 在 O(n^3) 的时间复杂度内求出。将所有区间 [X_{i-1}-iz,X_i-iz] 的左右端点排序,从小到大的第 i 个端点记做 v_i。设 dp_{i,j} 表示 [v_1,v_i] 这段区间内分布了 y_1\le y_2\le\ldots\le y_j 的概率,我们如下转移:

  • \forall p\in[j+1,k],[v_{i-1},v_i]\in[X_{p-1}-pz,X_p-pz],使 dp_{i,k}\gets dp_{i,k}+dp_{i-1,j}\dfrac{1}{(k-j+1)!}\prod\limits_{j<p\le k}\dfrac{v_i-v_{i-1}}{X_p-X_{p-1}}

现在考虑变化的 z。由于每个 v_i 都可视作关于 z 的线性函数,在所有区间端点的大小关系不变时,观察递推形式,发现 f(z) 是一个关于 z 的不超过 n 次的多项式函数。我们求出区间内 n+1 个不同的点插值得到 f(z) 在该区间的解析式后,再用微积分基本定理求 f(z) 在这段区间的积分即可。

由于所有区间的端点的大小关系只有 O(n^2) 种,我们可以在 O(n^6) 时间解决本题。

实现

#include <bits/stdc++.h>

typedef long long int ll;
const int maxn = 23, mod = 998244353;

ll gcd(ll a, ll b){
    return b ? gcd(b, a % b) : a;
}

int qpow(int a, int b){
    int res = 1;
    while(b){
        if(b & 1) res = (ll)res * a % mod;
        a = (ll)a * a % mod, b >>= 1;
    }
    return res;
}

struct frac{
    ll x, y;
    frac(ll a = 0){x = a, y = 1; }
    frac(ll a, ll b){
        ll g = gcd(a, b);
        a /= g, b /= g;
        if(b < 0) x = -a, y = -b;
        else x = a, y = b;
    }
    void out(void){
        std::printf("%lld/%lld", x, y);
    }
    int val(void){
        return (x % mod) * qpow(y % mod, mod - 2) % mod;
    }
};
frac operator-(const int &a, const frac &b){
    return frac(a * b.y - b.x, b.y);
}
frac operator+(const frac &a, const frac &b){
    return frac(a.x * b.y + b.x * a.y, a.y * b.y);
}
frac operator*(const frac &a, const int &b){
    return frac(a.x * b, a.y);
}
frac operator/(const frac &a, const int &b){
    return frac(a.x, a.y * b);
}
bool operator==(const frac &a, const frac &b){
    return a.x == b.x && a.y == b.y;
}
bool operator!=(const frac &a, const frac &b){
    return a.x != b.x || a.y != b.y;
}
bool operator<(const frac &a, const frac &b){
    return a.x * b.y < b.x * a.y;
}

int eval(const std::vector<int> &f, int z){
    int res = 0, x = 1;
    for(int i = 0; i < (int)f.size(); ++i)
        res = (res + (ll)f[i] * x) % mod, x = (ll)x * z % mod;
    return res;
}

int inv[maxn];
std::vector<int> integral(const std::vector<int> &f){
    std::vector<int> g(f.size() + 1);
    for(int i = f.size(); i; --i)
        g[i] = (ll)f[i - 1] * inv[i] % mod;
    return g;
}

std::vector<int> operator+(const std::vector<int> &f, const std::vector<int> &g){
    std::vector<int> h(f);
    if(g.size() > h.size()) h.resize(g.size());
    for(int i = 0; i < (int)g.size(); ++i)
        h[i] = (h[i] + g[i]) % mod;
    return h;
}

std::vector<int> operator*(const std::vector<int> &f, const std::vector<int> &g){
    std::vector<int> h(f.size() + g.size() - 1);
    for(int i = 0; i < (int)f.size(); ++i)
        for(int j = 0; j < (int)g.size(); ++j)
            h[i + j] = (h[i + j] + (ll)f[i] * g[j]) % mod;
    return h;
}

std::vector<int> operator*(const std::vector<int> &f, const int &b){
    std::vector<int> h(f);
    for(int i = 0; i < (int)h.size(); ++i)
        h[i] = (ll)h[i] * b % mod;
    return h;
}

int n, X[maxn], l[maxn], r[maxn];
std::pair<int, bool> v[maxn << 1];

void init(const frac &z){
    for(int i = 1; i <= n; ++i){
        v[i * 2 - 1] = std::make_pair(i, false);
        v[i * 2] = std::make_pair(i, true);
    }
    std::sort(v + 1, v + 1 + n * 2, [&z](const std::pair<int, bool> &a, const std::pair<int, bool> &b){
        return X[a.first + a.second - 1] - z * a.first != X[b.first + b.second - 1] - z * b.first ?
            X[a.first + a.second - 1] - z * a.first < X[b.first + b.second - 1] - z * b.first :
            a.first > b.first;
    });
    for(int i = 1; i <= n * 2; ++i)
        if(v[i].second) r[v[i].first] = i;
        else l[v[i].first] = i;
}

frac next(void){
    frac res = 1e6;
    for(int i = 1; i < n * 2; ++i)
        if(v[i].first < v[i + 1].first && v[i].first + v[i].second != v[i + 1].first + v[i + 1].second)
            res = std::min(res, frac(X[v[i + 1].first + v[i + 1].second - 1] - X[v[i].first + v[i].second - 1], v[i + 1].first - v[i].first));
    return res;
}

bool check(void){
    int last = 0;
    for(int i = 1; i <= n; ++i)
        if(last < r[i]) last = std::max(last, l[i]);
        else return false;
    return true;
}

int leninv[maxn], dp[maxn << 1][maxn], y[maxn];
int f(int z){
    for(int i = 1; i <= n * 2; ++i)
        y[i] = (X[v[i].first + v[i].second - 1] - (ll)z * v[i].first) % mod;
    dp[1][0] = 1;
    for(int i = 2; i <= n * 2; ++i)
        for(int j = 0; j <= n; ++j){
            dp[i][j] = dp[i - 1][j];
            for(int k = j, mul = 1; k >= 1 && l[k] <= i - 1 && r[k] >= i; --k){
                mul = (ll)mul * (y[i] - y[i - 1]) % mod * leninv[k] % mod * inv[j - k + 1] % mod;
                dp[i][j] = (dp[i][j] + (ll)dp[i - 1][k - 1] * mul) % mod;
            }
        }
    return dp[n * 2][n];
}

int solve(int a, int b){
    std::vector<int> F;
    for(int i = 0; i <= n; ++i){
        std::vector<int> G{1};
        int mul = 1;
        for(int j = 0; j <= n; ++j)
            if(i != j)
                G = G * std::vector<int>{-j, 1}, mul = (ll)mul * (i - j) % mod;
        mul = (ll)f(i) * qpow(mul, mod - 2) % mod;
        F = F + G * mul;
    }
    F = integral(F);
    return (eval(F, b) - eval(F, a)) % mod;
}

int main(){
    std::scanf("%d", &n);
    for(int i = 0; i <= n; ++i)
        std::scanf("%d", X + i);
    for(int i = 1; i <= n + 1; ++i)
        leninv[i] = qpow(X[i] - X[i - 1], mod - 2), inv[i] = qpow(i, mod - 2);
    frac l(0), r; int ans = 0;
    while(true){
        init(l), r = next();
        if(!check()) break;
        ans = (ans + solve(l.val(), r.val())) % mod;
        l = r;
    }
    std::printf("%d\n", (ans + mod) % mod);
    return 0;
}

来源

AtCoder Regular Contest 113 F - Social Distance

评论