Codeforces 551 DIV2 F - Serval and Bonus Problem (R2800)

こういうのに慣れて行きたい。

問題へのリンク

問題概要

長さ の線分上に、ランダムな区間を 個とったときの、区間が 本以上重なっている部分の長さの期待値を求めよ (998244353 で割った余りの形式で)。

なお、区間のランダムな選び方とは、線分から 2 点ずつを一様ランダムに選び、その 2 点を繋ぐものとする。

制約

• 

考えたこと

まず の値に大きな意味はなくて、 としてよい。 として期待値を求めて、それを最後に 倍すればよい。

問題は、線分上の各点 に対して、 が 本以上に被覆される確率密度関数 を求めて、それを で積分すればよいと考えたくなる...が、よくわからない式になってしまう。今回実は点 の座標値自体の情報を直接扱わずに解くことができる。

単純に、 個の区間を一様ランダムに並べたとき、一般には 個の部分区間に分割されるが、 本以上に被覆される部分区間が何個あるかの期待値を求めれば OK。

個の値の順列をすべて走査して、それぞれの「 本以上に被覆された部分区間の個数」の総和を求めて、最後に で割ることにする。さらに、 個の区間のうち左端が の順に登場するものに限定して最後に を掛けることにする。さらに、各区間に対して左端と右端を swap できるので を掛ける。

さて、これは箱根駅伝みたいな DP でできる！！！

• dp1[ i ][ j ] := N 区間のうち、左端が i 個登場し、そのうちの j 個が右端終了状態であるような個数

• dp2[ i ][ j ] := N 区間のうち、左端が i 個登場し、そのうちの j 個が右端終了状態である段階での、被覆 以上の部分区間の個数の総和

とすると

• dp1[ 0 ][ 0 ] = 1
• dp1[ i + 1 ][ j ] += dp1[ i ][ j ]
• dp1[ i ][ j + 1 ] += (i - j) × dp1[ i ][ j ]
• dp2[ i + 1 ][ j ] += dp2[ i ][ j ]
• dp2[ i + 1 ][ j ] += (i - j >= K ? dp1[ i ][ j ] : 0)
• dp2[ i ][ j + 1 ] += (i - j) × dp2[ i ][ j ]
• dp2[ i ][ j + 1 ] += (i - j) × (i - j >= K ? dp1[ i ][ j ] : 0)

という感じになる。習い覚えた modint を使いまくった。

#include <iostream>
#include <vector>
using namespace std;


template<int MODULO> struct Fp {
    long long val;

    constexpr Fp(long long v = 0) noexcept : val(v % MODULO) {
        if (val < 0) v += MODULO;
    }
    constexpr Fp operator - () const noexcept {
        return val ? MODULO - 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 >= MODULO) val -= MODULO;
        return *this;
    }
    constexpr Fp& operator -= (const Fp& r) noexcept {
        val -= r.val;
        if (val < 0) val += MODULO;
        return *this;
    }
    constexpr Fp& operator *= (const Fp& r) noexcept {
        val = val * r.val % MODULO;
        return *this;
    }
    constexpr Fp& operator /= (const Fp& r) noexcept {
        long long a = r.val, b = MODULO, 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 % MODULO;
        if (val < 0) val += MODULO;
        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;
    }
};

template<int MODULO> constexpr ostream& operator <<
(ostream &os, const Fp<MODULO>& x) noexcept {
    return os << x.val;
}
template<int MODULO> constexpr istream& operator >>
(istream &is, Fp<MODULO>& x) noexcept {
    return is >> x.val;
}

template<int MODULO> constexpr Fp<MODULO> modpow
(const Fp<MODULO> &a, long long n) noexcept {
    if (n == 0) return 1;
    auto t = modpow(a, n / 2);
    t = t * t;
    if (n & 1) t = t * a;
    return t;
}


template<int MODULO> struct BiCoef {
    vector<Fp<MODULO> > fac, inv, finv;
    constexpr BiCoef(int n = 210000) noexcept : fac(n, 1), inv(n, 1), finv(n, 1) {
        for(int i = 2; i < n; i++){
            fac[i] = fac[i-1] * i;
            inv[i] = -inv[MODULO%i] * (MODULO/i);
            finv[i] = finv[i-1] * inv[i];
        }
    }
    constexpr Fp<MODULO> com(int n, int k) const noexcept {
        if (n < k || n < 0 || k < 0) return 0;
        return fac[n] * finv[k] * finv[n-k];
    }
};




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

int main() {
    BiCoef<MOD> bc;
    long long N, K, L;
    cin >> N >> K >> L;

    vector<vector<mint> > dp1(N+1, vector<mint>(N+1, 0));
    auto dp2 = dp1;
    dp1[0][0] = 1;

    for (int i = 0; i <= N; ++i) {
        for (int j = 0; j <= i; ++j) {
            if (i + 1 <= N) {
                dp1[i+1][j] += dp1[i][j];
                dp2[i+1][j] += dp2[i][j];
                if (i-j >= K) dp2[i+1][j] += dp1[i][j];
            }
            if (j + 1 <= i) {
                dp1[i][j+1] += dp1[i][j] * (i - j);
                dp2[i][j+1] += dp2[i][j] * (i - j);
                if (i-j >= K) dp2[i][j+1] += dp1[i][j] * (i - j);
            }
        }
    }
    mint res = dp2[N][N] * bc.fac[N] * bc.finv[N*2+1] * modpow(mint(2), N) * L;
    cout << res << endl;
}