Codeforces Round #551 (Div. 2) F. Serval and Bonus Problem

https://codeforces.com/contest/1153/problem/F

Without loss of generality, we assume \(L=1\). The probability that x is covered by a single interval is \(2x(1-x)\). Let \(f(x)=2x(1-x)\), then the expectation of the length covered by exactly i intervals can be written by

\[ \int_{0}^{1} \binom{N}{i} f(x)^i (1-f(x))^{N-i} dx \]

If we have

\[ \int_{0}^{1} f(x)^i dx \]

the answer can be computed easily. The coefficients of \(f(x)^i\) can be computed by \(f(x)^{i-1}\) in O(N) time. The total time complexity is O(N^2).

#include <iostream>
#include <algorithm>
#include <vector>
#include <string>

using namespace std;

const int MOD = 998244353;

struct mint {
    int n;
    mint(int n_ = 0) : n(n_) {}
};

mint operator+(mint a, mint b) { a.n += b.n; if (a.n >= MOD) a.n -= MOD; return a; }
mint operator-(mint a, mint b) { a.n -= b.n; if (a.n < 0) a.n += MOD; return a; }
mint operator*(mint a, mint b) { return (long long)a.n * b.n % MOD; }
mint &operator+=(mint &a, mint b) { return a = a + b; }
mint &operator-=(mint &a, mint b) { return a = a - b; }
mint &operator*=(mint &a, mint b) { return a = a * b; }

mint modpow(mint a, int b) {
  mint res = 1;
  while (b > 0) {
    if (b & 1) res *= a;
    a *= a;
    b >>= 1;
  }
  return res;
}

mint modinv(mint a) {
  return modpow(a, MOD - 2);
}

mint F[100001] = {1, 1};
mint R[100001] = {1, 1};
mint I[100001] = {0, 1};

mint C(int n, int r) {
  if (n < 0 || r < 0 || n < r) return 0;
  return F[n] * R[n - r] * R[r];
}

void init() {
  for (int i = 2; i <= 100000; i++) {
    I[i] = I[MOD % i] * (MOD - MOD / i);
    F[i] = F[i - 1] * i;
    R[i] = R[i - 1] * I[i];
  }
}

mint dp[2020][4040];
mint idp[2020];

int main() {
  init();
  int N, K, L;
  cin >> N >> K >> L;
  // P(x) = 2x(1-x) = -2x^2 + 2x
  dp[0][0] = 1;
  for (int i = 0; i < N; i++) {
    for (int j = 0; j <= 4000; j++) {
      dp[i + 1][j + 2] -= 2 * dp[i][j];
      dp[i + 1][j + 1] += 2 * dp[i][j];
    }
  }
  for (int i = 0; i <= N; i++) {
    for (int j = 0; j <= 4000; j++) {
      idp[i] += I[j + 1] * dp[i][j];
    }
  }
  mint ans = 0;
  for (int i = K; i <= N; i++) {
    // \int_0^1 C(N,i) * P(x)^i * (1-P(x))^{N-i}
    for (int j = 0; j <= N - i; j++) {
      if (j % 2 == 0) {
        ans += C(N,i) * C(N-i,j) * idp[i + j];
      } else {
        ans -= C(N,i) * C(N-i,j) * idp[i + j];
      }
    }
  }
  ans *= L;
  cout << ans.n << endl;
}