けんちょんの競プロ精進記録

競プロの精進記録や小ネタを書いていきます

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 と初期化しておく
  • 与えられる各区間  \lbrack L, R \rbrack に対して、
    •  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;
}