2020-10-11

HHKB プログラミングコンテスト 2020 F - Random Max (橙色, 600 点)

時々やってくる積分ゲー。同時に多項式ゲーでもあった。

問題へのリンク

問題概要

$N$ 個の区間 $\lbrack L_{i}, R_{i} \rbrack$ が与えられる。これらの区間から一様分布にしたがって点をとってくる (連続値)。

各点の座標の最大値の期待値を $(N+1)! \times \prod_{i=1}^{N}(R_{i} - L_{i})$ をかけた値 (整数値になる) を 1000000007 で割ったあまりを求めよ。

制約

$1 \le N \le 1000$
$0 \le L_{i} \lt R_{i} \le 10^{9}$

考えたこと

最初は $\prod_{i=1}^{N}(R_{i} - L_{i})$ はともかく、なぜ $(N+1)!$ が...と思った。でも、次のことを思い出すと「確かに...」となった。なお、この性質を生かした解法もある模様 (これやこれ)。

高校数学の美しい物語：最大値と最小値の分布（一般論と例）

f:id:drken1215:20201011115209p:plain

よって、 $(N+1)!$ をかけておけば確かに整数になりそう (区間が $N$ 個被っている範囲もあれば、 $1, 2, \dots, N-1$ 個だけが被っている範囲もありうるので、 $(N+1)$ をかけるだけでは不十分で $(N+1)!$ をかける必要がある)。

それはそうと、もしこの問題が離散値の問題だったら、次のように考える典型テクニックがある

$F(x)$ := 最大値が $x$ 以下となる確率

としてあげると、最大値が $x$ となる確率は $F(x) - F(x-1)$ で表せるので、求める期待値は

$\sum_{x} x(F(x) - F(x-1))$

となる。 $x$ が連続量の場合もこれを参考にして、期待値は

$\int xF(x)' dx$

で求められる！！！これをどうにかする。

解法 (1)：形式的冪級数で殴る

まず区間の左端と右端をかき集めると $2N$ 個になる。よって、 $2N$ 個の区間それぞれについて求めて合算していけば OK。

区間 $\lbrack l, r \rbrack$ について、 $F(x)$ は次のように求められる。なおここでは、 $\prod_{i=1}^{N}(R_{i} - L_{I})$ はすでに掛けられた状態で考えることとし、 $(N-1)!$ を掛けるのは最後に行うこととする。

はじめは $F(x) = 1$ と初期化しておく
与えられる各区間に対して、
- $L \ge r$ となっていた場合には $F(x) = 0$ と更新
- $R \le l$ となっていた場合には $F(x) *= (R - L)$ と更新
- それ以外の場合には $F(x) *= (x - L)$ と更新

こうして求めた $F(x)$ に対して、形式的冪級数を用いて $\int xF(x)' dx$ を計算して得た不定積分関数を $G(x)$ としたとき、 $G(r) - G(l)$ を加算すれば OK。

愚直にやると計算量は $O(N^{3} \log N)$ となるけど、それでも間に合った。 $F(x)$ を求める部分で掛け算の順序を少し工夫すると $O(N^{2} (\log N)^{2})$ にはなる。

atcoder.jp

解法 (2)：後ろからやる

実装的にも計算量的にも、x 座標が大きい方から見ていくととても楽になる。スタート時点では $F(x)$ は定数になっている。そしてある区間の右端にぶつかるたびに、そこの区間に対して $F(x) *= \frac{x - L}{R - L}$ をしていけば OK。そしてどこかいずれかの区間の左端に触れた時点で終了。

このようにした場合の計算量は $O(N^{2} \log N)$ となる。

atcoder.jp

解法 (3)：DP にする

解法 (1)(2) では、次の積分計算を直接 FPS ライブラリで殴り倒した。しかし FPS 計算は

$(x - l)$ という形の式をかける (DP で表現できる)
微分する
$x$ をかける
積分する

というだけなので、上手くやれば DP だけでできる。たとえば $x-l$ をかけるのは次のような遷移 (dp が ndp になる想定) でできる。

ndp[ i ] = dp[ i - 1 ] - l × dp[ i ]

そして一つ一つの操作は $O(N)$ でできるので、解法 (2) の考え方と合わせれば、全体の計算量も $O(N^{2})$ におさえられる。

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

// modint
template<int MOD> struct Fp {
    long long val;
    constexpr Fp(long long v = 0) noexcept : val(v % MOD) {
        if (val < 0) val += MOD;
    }
    constexpr int getmod() const { return MOD; }
    constexpr Fp operator - () const noexcept {
        return val ? MOD - val : 0;
    }
    constexpr Fp operator + (const Fp& r) const noexcept { return Fp(*this) += r; }
    constexpr Fp operator - (const Fp& r) const noexcept { return Fp(*this) -= r; }
    constexpr Fp operator * (const Fp& r) const noexcept { return Fp(*this) *= r; }
    constexpr Fp operator / (const Fp& r) const noexcept { return Fp(*this) /= r; }
    constexpr Fp& operator += (const Fp& r) noexcept {
        val += r.val;
        if (val >= MOD) val -= MOD;
        return *this;
    }
    constexpr Fp& operator -= (const Fp& r) noexcept {
        val -= r.val;
        if (val < 0) val += MOD;
        return *this;
    }
    constexpr Fp& operator *= (const Fp& r) noexcept {
        val = val * r.val % MOD;
        return *this;
    }
    constexpr Fp& operator /= (const Fp& r) noexcept {
        long long a = r.val, b = MOD, u = 1, v = 0;
        while (b) {
            long long t = a / b;
            a -= t * b, swap(a, b);
            u -= t * v, swap(u, v);
        }
        val = val * u % MOD;
        if (val < 0) val += MOD;
        return *this;
    }
    constexpr bool operator == (const Fp& r) const noexcept {
        return this->val == r.val;
    }
    constexpr bool operator != (const Fp& r) const noexcept {
        return this->val != r.val;
    }
    friend constexpr istream& operator >> (istream& is, Fp<MOD>& x) noexcept {
        is >> x.val;
        x.val %= MOD;
        if (x.val < 0) x.val += MOD;
        return is;
    }
    friend constexpr ostream& operator << (ostream& os, const Fp<MOD>& x) noexcept {
        return os << x.val;
    }
    friend constexpr Fp<MOD> modpow(const Fp<MOD>& r, long long n) noexcept {
        if (n == 0) return 1;
        auto t = modpow(r, n / 2);
        t = t * t;
        if (n & 1) t = t * r;
        return t;
    }
    friend constexpr Fp<MOD> modinv(const Fp<MOD>& r) noexcept {
        long long a = r.val, b = MOD, u = 1, v = 0;
        while (b) {
            long long t = a / b;
            a -= t * b, swap(a, b);
            u -= t * v, swap(u, v);
        }
        return Fp<MOD>(u);
    }
};

const int MOD = 1000000007;
using mint = Fp<MOD>;

// 代入
mint eval(const vector<mint> &f, const mint& v){
    mint res = 0;
    for (int i = (int)f.size()-1; i >= 0; --i) {
        res *= v;
        res += f[i];
    }
    return res;
}

// 微分
vector<mint> diff(const vector<mint>& f) {
    int n = (int)f.size();
    vector<mint> res(n-1);
    for (int i = 1; i < n; ++i) res[i-1] = f[i] * i;
    return res;
}

// 積分
vector<mint> integral(const vector<mint>& f) {
    int n = (int)f.size();
    vector<mint> res(n+1, 0);
    for (int i = 0; i < n; ++i) res[i+1] = f[i] / (i+1);
    return res;
}

// x - v をかける
vector<mint> mul(const vector<mint> &f, const mint &v) {
    int n = (int)f.size();
    vector<mint> res(n+1, 0);
    for (int i = 0; i <= n; ++i) {
        if (i > 0) res[i] += f[i-1];
        if (i < n) res[i] += f[i] * (-v);
    }
    return res;
}

int main() {
    int N;
    cin >> N;
    vector<long long> L(N), R(N);
    long long finL = 0; // max(L)
    for (int i = 0; i < N; ++i) {
        cin >> L[i] >> R[i];
        finL = max(finL, L[i]);
    }

    vector<pair<long long, int>> alt;
    mint all = 1;
    for (int i = 0; i < N; ++i) {
        alt.emplace_back(R[i], i);
        all *= R[i] - L[i];
    }
    alt.emplace_back(finL, -1); // 最後の区間
    sort(alt.begin(), alt.end(), greater<pair<long long, int>>());
    vector<mint> F(1, all);

    mint res = 0;
    for (int i = 0; i + 1 < alt.size(); ++i) {
        long long p = alt[i+1].first, q = alt[i].first; // 積分区間
        if (q <= finL) break;
        long long l = L[alt[i].second], r = q; // 考えている区間
        F = mul(F, l);
        for (auto &v : F) v /= r - l;
        auto G = integral(mul(diff(F), 0));
        auto add = eval(G, q) - eval(G, p);
        res += add;
    }
    for (int n = 1; n <= N + 1; ++n) res *= n;
    cout << res << endl;
}