SRM626 div1 Easy FixedDiceGameDiv1

概要

Aliceはa個のb面ダイス、Bobはc個のd面ダイスを投げる。このとき、ダイスの目の合計が大きい方を勝ちとする(引き分けはどちらの勝ちでもない)。Aliceが勝つことがありえないときは-1を出力しろ。そうでないときはAliceが勝った条件下でAliceのダイスの目の合計の期待値を求めよ。

解法

まず a*b < c のときAliceが勝つことは不可能でこれは簡単に判定できるのでAliceが絶対に勝てないときは考えない。
Aliceが1を出したときにBobが1未満を出す確率、Aliceが2を出したときにBobが2未満を出す確率…の和がAliceが勝つ確率$P$でAliceが出したダイスの値をそれぞれ掛けたものの和がAliceが勝ったという前提なしでAliceが勝つときのダイスの目の期待値$E'$になる。Aliceがxを出す確率をdp1[x]、Bobがx未満を出す確率をdp2[x]と書くと、
E'= 1×dp1[1]×dp2[1] + 2×dp1[2]×dp2[2] + 3×dp1[3]×dp2[3] + …
P = dp1[1]×dp2[1] + dp1[2]×dp2[2] + dp1[3]×dp2[3] + …
と書ける。Aliceが勝っている前提での期待値$E$は E = E'/P となる。

dp1[x]はDPで計算することができる。dp[i][j] = (i回目までで合計jを出す期待値) とすると for(k=1; k<=b; ++k) dp[i+1][j+k] += dp[i][j]/b で dp[i][j] から更新できこれはO(a^2b^2)でできる。
dp2[x]はこれと同様にDPをした後、累積和を取ることで計算できる。

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
typedef vector<int> VI;
typedef vector<VI> VVI;
typedef vector<ll> VL;
typedef vector<VL> VVL;
typedef pair<int, int> PII;

#define FOR(i, a, n) for (ll i = (ll)a; i < (ll)n; ++i)
#define REP(i, n) FOR(i, 0, n)
#define ALL(x) x.begin(), x.end()
#define IN(a, b, x) (a<=x&&x<b)
#define MP make_pair
#define PB push_back
const int INF = (1LL<<30);
const ll LLINF = (1LL<<60);
const double PI = 3.14159265359;
const double EPS = 1e-12;
const int MOD = 1000000007;
//#define int ll

template <typename T> T &chmin(T &a, const T &b) { return a = min(a, b); }
template <typename T> T &chmax(T &a, const T &b) { return a = max(a, b); }

int dx[] = {0, 1, 0, -1}, dy[] = {1, 0, -1, 0};

double dp1[55][2605], dp2[55][2605], dp3[2505];
class FixedDiceGameDiv1 {
  public:
  double getExpectation(int a, int b, int c, int d)
  {
    if(a*b < c) return -1;
    memset(dp1, 0, sizeof(dp1)); memset(dp2, 0, sizeof(dp2)); memset(dp3, 0, sizeof(dp3));
    dp1[0][0] = dp2[0][0] = 1;

    // dp1
    REP(i, a) {
      REP(j, a*b+1) {
        if(!dp1[i][j]) continue;
        // dp1[i][j] から dp1[i+1][j+k] に配る
        FOR(k, 1, b+1) dp1[i+1][j+k] += dp1[i][j]/b;
      }
    }

    // dp2
    REP(i, c) {
      REP(j, c*d+1) {
        if(!dp2[i][j]) continue;
        FOR(k, 1, d+1) dp2[i+1][j+k] += dp2[i][j]/d;
      }
    }

    dp3[0] = dp2[c][0];
    FOR(i, 1, 2501) dp3[i] += dp3[i-1] + dp2[c][i];

    double ret = 0.0, res = 0.0;
    FOR(i, 1, a*b+1) {
      ret += i*dp1[a][i]*dp3[i-1];
      res += dp1[a][i]*dp3[i-1];
    }
    if(res == 0) return -1.0;
    return ret/res;
  }
};