Problem 982: The Third Dice
View on Project EulerProject Euler Problem 982 Solution
EulerSolve provides an optimized solution for Project Euler Problem 982, The Third Dice, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Three fair six-sided dice are rolled. A player who sees all three faces chooses one die to hide. The second player then sees the two remaining faces and chooses the option with the larger expected payoff: either take the hidden die, whose value is uncertain, or take the better of the two visible dice, whose value is known exactly. The task is to determine the optimal expected payoff when the hider plays as well as possible. The local implementations do not guess a simple rule such as “always hide the largest die”. They solve the game exactly. The key observation is that once a visible pair has been fixed, the second player only needs two quantities: the conditional mean of the hidden die and the larger visible face. That turns the full game into a compact linear program. Mathematical Approach Write an ordered state as \(s=(d_1,d_2,d_3)\in\{1,\ldots,6\}^3\), and let \(a\in\{1,2,3\}\) denote which die is hidden. All \(216\) ordered states are equally likely, so \(p_s=1/216\). The observations that actually matter The state space is ordered, but the observation is not: after one die is hidden, the second player only sees the unordered pair of visible faces....
Detailed mathematical approach
Problem Summary
Three fair six-sided dice are rolled. A player who sees all three faces chooses one die to hide. The second player then sees the two remaining faces and chooses the option with the larger expected payoff: either take the hidden die, whose value is uncertain, or take the better of the two visible dice, whose value is known exactly. The task is to determine the optimal expected payoff when the hider plays as well as possible.
The local implementations do not guess a simple rule such as “always hide the largest die”. They solve the game exactly. The key observation is that once a visible pair has been fixed, the second player only needs two quantities: the conditional mean of the hidden die and the larger visible face. That turns the full game into a compact linear program.
Mathematical Approach
Write an ordered state as \(s=(d_1,d_2,d_3)\in\{1,\ldots,6\}^3\), and let \(a\in\{1,2,3\}\) denote which die is hidden. All \(216\) ordered states are equally likely, so \(p_s=1/216\).
The observations that actually matter
The state space is ordered, but the observation is not: after one die is hidden, the second player only sees the unordered pair of visible faces. Therefore the observation space is
$$\mathcal{O}=\{\{u,v\}:1\le u\le v\le 6\},\qquad |\mathcal{O}|=\binom{7}{2}=21.$$
For each state-action pair \((s,a)\), let \(o(s,a)\in\mathcal{O}\) be that visible multiset, let \(h(s,a)\) be the hidden face, and define
$$m(o)=\max(o)$$
as the best visible die the second player could take immediately. The smaller visible die is irrelevant to optimal play, because choosing it is always dominated by choosing the larger visible die.
A hiding strategy written as probabilities
The hider may randomize. Let \(x_{s,a}\ge 0\) be the probability of hiding die \(a\) in state \(s\). For every state, these probabilities must sum to \(1\):
$$\sum_{a=1}^3 x_{s,a}=1\qquad(s\in\{1,\ldots,6\}^3).$$
For a fixed observation \(o\), the total probability of producing that observation is
$$q_o=\sum_{\substack{(s,a)\\o(s,a)=o}} p_s\,x_{s,a}.$$
Conditioned on seeing \(o\), the expected value of the hidden die is
$$\mu_o=\frac{1}{q_o}\sum_{\substack{(s,a)\\o(s,a)=o}} p_s\,h(s,a)\,x_{s,a}\qquad(q_o>0).$$
So after observation \(o\), the second player compares the sure visible value \(m(o)\) with the conditional hidden expectation \(\mu_o\). The contribution of \(o\) to the overall expected payoff is therefore
$$q_o\max(\mu_o,m(o)).$$
Linearizing the opponent's best response
The expression above is piecewise linear, but it can be written exactly with one auxiliary variable \(T_o\) for each observation class. First note that
$$q_o\mu_o=\sum_{\substack{(s,a)\\o(s,a)=o}} p_s\,h(s,a)\,x_{s,a},$$
and also
$$q_om(o)=m(o)\sum_{\substack{(s,a)\\o(s,a)=o}} p_s\,x_{s,a}.$$
Hence it is enough to impose
$$T_o\ge \sum_{\substack{(s,a)\\o(s,a)=o}} p_s\,h(s,a)\,x_{s,a},$$
$$T_o\ge m(o)\sum_{\substack{(s,a)\\o(s,a)=o}} p_s\,x_{s,a}.$$
If the objective is to minimize \(\sum_{o\in\mathcal{O}}T_o\), every \(T_o\) is forced down to the larger of those two quantities. Therefore
$$\sum_{o\in\mathcal{O}}T_o=\sum_{o\in\mathcal{O}} q_o\max(\mu_o,m(o)),$$
which is exactly the second player’s expected payoff. The whole problem is thus the linear program
$$\min \sum_{o\in\mathcal{O}}T_o$$
subject to the state-normalization equations and the two inequalities above for every observation.
Worked example: seeing \(\{2,5\}\)
Take the observation \(o=\{2,5\}\). Then \(m(o)=5\). This observation can come from many ordered state-action pairs, for example by hiding the first die in \((4,2,5)\) or by hiding the second die in \((2,6,5)\). The strategy decides how much probability mass each such pair contributes.
If the resulting conditional mean of the hidden die is \(\mu_o=29/6\), then the sure visible \(5\) is better, so this observation contributes \(5q_o\). If instead \(\mu_o=16/3\), then taking the hidden die is better, and the contribution becomes \((16/3)q_o\). The two inequalities for \(T_o\) cover both situations without any case split inside the solver.
How the Code Works
Building the tableau
The C++, Python, and Java implementations enumerate all \(216\) ordered three-die states and all \(3\) hide choices, giving \(648\) state-action columns. For each column they precompute the observation class, the hidden face, and the larger visible face. They also keep the smaller two-dice version of the same model, which has \(36\) states, \(2\) hide choices, and \(6\) observation classes.
The linear program uses one normalization row per state and two rows per observation class. For the target instance that is \(216+2\cdot 21=258\) rows. It also uses \(21\) auxiliary observation variables, \(42\) slack variables for the inequalities, and \(216\) artificial variables for Phase I.
Phase I, Phase II, and the validation case
The state-normalization equations do not come with an obvious feasible basis, so the implementation starts with a standard two-phase simplex. Phase I minimizes the sum of artificial variables; feasibility is certified when that optimum reaches \(0\). The artificial columns are then removed, any row that has become redundant is dropped, and Phase II minimizes \(\sum T_o\) by maximizing its negative in tableau form.
Before solving the three-dice target, the same machinery is run on the two-dice instance and checked against the known value \(145/36\). The C++ implementation also parallelizes the coefficient fill for the observation rows, while the Python and Java implementations perform the same arithmetic serially.
Complexity Analysis
Tableau construction is linear in the number of state-action pairs, because each column contributes to exactly one hidden-value row and one visible-value row. For the target case this means \(648\) such contributions, plus fixed row and column setup.
The main cost is the simplex pivoting. In the three-dice instance, Phase I starts from a \(258\times 927\) tableau, and Phase II continues with \(711\) columns after the artificial variables are removed. As with any simplex implementation, the theoretical worst case is exponential in the number of pivots, but at this matrix size the exact solve is practical. Memory usage is \(O(mn)\) for the tableau.
Footnotes and References
- Problem page: https://projecteuler.net/problem=982
- Linear programming: Wikipedia - Linear programming
- Simplex algorithm: Wikipedia - Simplex algorithm
- Expected value: Wikipedia - Expected value
- Conditional expectation: Wikipedia - Conditional expectation
Problem 982 source code
C++
#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <thread>
#include <vector>
using namespace std;
static constexpr double EPS = 1e-10;
template <typename Func>
static void parallel_for(int total, int threads, Func func) {
if (threads <= 1 || total < 128) {
for (int i = 0; i < total; ++i) func(i);
return;
}
threads = min(threads, total);
int chunk = (total + threads - 1) / threads;
vector<thread> pool;
pool.reserve(threads);
for (int t = 0; t < threads; ++t) {
int start = t * chunk;
int end = min(total, start + chunk);
if (start >= end) break;
pool.emplace_back([=, &func]() {
for (int i = start; i < end; ++i) func(i);
});
}
for (auto &th : pool) th.join();
}
struct Simplex {
int m;
int n;
vector<vector<double>> T;
vector<int> basis;
Simplex(int rows, int cols)
: m(rows), n(cols), T(rows + 1, vector<double>(cols + 1, 0.0)), basis(rows, -1) {}
void pivot(int r, int c) {
double inv = 1.0 / T[r][c];
for (int j = 0; j <= n; ++j) T[r][j] *= inv;
for (int i = 0; i <= m; ++i) {
if (i == r) continue;
double factor = T[i][c];
if (fabs(factor) <= EPS) continue;
for (int j = 0; j <= n; ++j) T[i][j] -= factor * T[r][j];
}
basis[r] = c;
}
void set_objective(const vector<double> &c) {
fill(T[m].begin(), T[m].end(), 0.0);
for (int j = 0; j < n; ++j) T[m][j] = -c[j];
T[m][n] = 0.0;
for (int i = 0; i < m; ++i) {
int var = basis[i];
if (var < 0) continue;
double coeff = c[var];
if (fabs(coeff) <= EPS) continue;
for (int j = 0; j <= n; ++j) T[m][j] += coeff * T[i][j];
}
}
bool solve() {
while (true) {
int enter = -1;
for (int j = 0; j < n; ++j) {
if (T[m][j] < -EPS) {
enter = j;
break;
}
}
if (enter == -1) return true;
double min_ratio = 0.0;
int leave = -1;
for (int i = 0; i < m; ++i) {
double a = T[i][enter];
if (a > EPS) {
double ratio = T[i][n] / a;
if (leave == -1 || ratio < min_ratio - 1e-12 ||
(fabs(ratio - min_ratio) <= 1e-12 && basis[i] < basis[leave])) {
min_ratio = ratio;
leave = i;
}
}
}
if (leave == -1) return false;
pivot(leave, enter);
}
}
double objective_value() const { return T[m][n]; }
};
struct GameData {
int dice = 0;
int actions = 0;
int states = 0;
int obs = 0;
double prob = 0.0;
vector<int> visible_max;
vector<int> action_obs;
vector<int> action_hidden;
};
static GameData build_game(int dice) {
GameData g;
g.dice = dice;
g.actions = dice;
if (dice == 2) {
g.states = 36;
g.obs = 6;
g.prob = 1.0 / 36.0;
g.visible_max.resize(g.obs);
for (int v = 0; v < g.obs; ++v) g.visible_max[v] = v + 1;
int n_action = g.states * g.actions;
g.action_obs.resize(n_action);
g.action_hidden.resize(n_action);
for (int a = 1; a <= 6; ++a) {
for (int b = 1; b <= 6; ++b) {
int s = (a - 1) * 6 + (b - 1);
int base = s * 2;
g.action_obs[base] = b - 1;
g.action_hidden[base] = a;
g.action_obs[base + 1] = a - 1;
g.action_hidden[base + 1] = b;
}
}
return g;
}
if (dice == 3) {
g.states = 216;
g.obs = 21;
g.prob = 1.0 / 216.0;
int obs_index[7][7];
for (int i = 0; i <= 6; ++i) for (int j = 0; j <= 6; ++j) obs_index[i][j] = -1;
g.visible_max.resize(g.obs);
int idx = 0;
for (int u = 1; u <= 6; ++u) {
for (int v = u; v <= 6; ++v) {
obs_index[u][v] = idx;
obs_index[v][u] = idx;
g.visible_max[idx] = v;
idx++;
}
}
int n_action = g.states * g.actions;
g.action_obs.resize(n_action);
g.action_hidden.resize(n_action);
for (int a = 1; a <= 6; ++a) {
for (int b = 1; b <= 6; ++b) {
for (int c = 1; c <= 6; ++c) {
int s = ((a - 1) * 6 + (b - 1)) * 6 + (c - 1);
int base = s * 3;
g.action_obs[base] = obs_index[b][c];
g.action_hidden[base] = a;
g.action_obs[base + 1] = obs_index[a][c];
g.action_hidden[base + 1] = b;
g.action_obs[base + 2] = obs_index[a][b];
g.action_hidden[base + 2] = c;
}
}
}
return g;
}
cerr << "Unsupported dice count: " << dice << "\n";
return g;
}
struct LPData {
Simplex lp;
int n_action;
int n_T;
int n_slack;
int n_art;
int art_start;
int slack_start;
LPData(int rows, int cols, int actions, int tvars, int slack, int art, int art_s, int slack_s)
: lp(rows, cols), n_action(actions), n_T(tvars), n_slack(slack), n_art(art),
art_start(art_s), slack_start(slack_s) {}
};
static LPData build_lp(const GameData &g, unsigned threads) {
int n_action = g.states * g.actions;
int n_T = g.obs;
int n_orig = n_action + n_T;
int n_slack = g.obs * 2;
int n_art = g.states;
int n_total = n_orig + n_slack + n_art;
int rows = g.states + n_slack;
int slack_start = n_orig;
int art_start = n_orig + n_slack;
LPData data(rows, n_total, n_action, n_T, n_slack, n_art, art_start, slack_start);
Simplex &lp = data.lp;
for (int s = 0; s < g.states; ++s) {
int row = s;
int base = s * g.actions;
for (int a = 0; a < g.actions; ++a) lp.T[row][base + a] = 1.0;
lp.T[row][art_start + s] = 1.0;
lp.T[row][n_total] = 1.0;
lp.basis[row] = art_start + s;
}
for (int o = 0; o < g.obs; ++o) {
for (int t = 0; t < 2; ++t) {
int row = g.states + o * 2 + t;
int tvar = n_action + o;
lp.T[row][tvar] = -1.0;
lp.T[row][slack_start + o * 2 + t] = 1.0;
lp.T[row][n_total] = 0.0;
lp.basis[row] = slack_start + o * 2 + t;
}
}
auto worker = [&](int idx) {
int obs = g.action_obs[idx];
int hidden = g.action_hidden[idx];
int row_hidden = g.states + obs * 2;
int row_visible = row_hidden + 1;
lp.T[row_hidden][idx] += g.prob * hidden;
lp.T[row_visible][idx] += g.prob * g.visible_max[obs];
};
parallel_for(n_action, static_cast<int>(threads), worker);
return data;
}
static void remove_artificial(Simplex &lp, int art_start, int art_count) {
if (art_count == 0) return;
vector<int> rows_to_remove;
vector<char> is_art(lp.n, 0);
for (int i = 0; i < art_count; ++i) is_art[art_start + i] = 1;
for (int i = 0; i < lp.m; ++i) {
int var = lp.basis[i];
if (var < 0 || !is_art[var]) continue;
int pivot_col = -1;
for (int j = 0; j < lp.n; ++j) {
if (is_art[j]) continue;
if (fabs(lp.T[i][j]) > EPS) {
pivot_col = j;
break;
}
}
if (pivot_col != -1) {
lp.pivot(i, pivot_col);
} else {
rows_to_remove.push_back(i);
}
}
if (!rows_to_remove.empty()) {
vector<char> drop(lp.m, 0);
for (int r : rows_to_remove) drop[r] = 1;
int new_m = lp.m - static_cast<int>(rows_to_remove.size());
vector<vector<double>> newT(new_m + 1, vector<double>(lp.n + 1, 0.0));
vector<int> new_basis(new_m, -1);
int r2 = 0;
for (int i = 0; i < lp.m; ++i) {
if (drop[i]) continue;
newT[r2] = lp.T[i];
new_basis[r2] = lp.basis[i];
r2++;
}
newT[new_m] = lp.T[lp.m];
lp.T.swap(newT);
lp.basis.swap(new_basis);
lp.m = new_m;
}
int new_n = lp.n - art_count;
vector<vector<double>> newT(lp.m + 1, vector<double>(new_n + 1, 0.0));
for (int i = 0; i <= lp.m; ++i) {
int col_new = 0;
for (int j = 0; j < lp.n; ++j) {
if (is_art[j]) continue;
newT[i][col_new++] = lp.T[i][j];
}
newT[i][new_n] = lp.T[i][lp.n];
}
lp.n = new_n;
lp.T.swap(newT);
}
static double solve_game(int dice, unsigned threads) {
GameData g = build_game(dice);
if (g.states == 0) return nan("");
LPData data = build_lp(g, threads);
Simplex &lp = data.lp;
vector<double> c1(lp.n, 0.0);
for (int i = 0; i < data.n_art; ++i) c1[data.art_start + i] = -1.0;
lp.set_objective(c1);
if (!lp.solve()) {
cerr << "Phase I unbounded\n";
return nan("");
}
double phase1 = lp.objective_value();
if (phase1 < -1e-8) {
cerr << "Phase I infeasible: " << phase1 << "\n";
return nan("");
}
remove_artificial(lp, data.art_start, data.n_art);
vector<double> c2(lp.n, 0.0);
for (int i = 0; i < data.n_T; ++i) c2[data.n_action + i] = -1.0;
lp.set_objective(c2);
if (!lp.solve()) {
cerr << "Phase II unbounded\n";
return nan("");
}
return -lp.objective_value();
}
static bool run_validation(unsigned threads) {
double expected = 145.0 / 36.0;
double got = solve_game(2, threads);
if (!isfinite(got) || fabs(got - expected) > 1e-7) {
cerr << "Validation failed: expected " << fixed << setprecision(9) << expected
<< ", got " << got << "\n";
return false;
}
return true;
}
int main(int argc, char **argv) {
ios::sync_with_stdio(false);
cin.tie(nullptr);
unsigned threads = thread::hardware_concurrency();
if (threads == 0) threads = 1;
threads = min(threads, 8u);
bool validate = true;
if (argc >= 2) threads = max(1u, static_cast<unsigned>(stoul(argv[1])));
if (argc >= 3) validate = (stoul(argv[2]) != 0);
if (validate && !run_validation(min(threads, 4u))) {
return 1;
}
double answer = solve_game(3, threads);
if (!isfinite(answer)) return 1;
cout << fixed << setprecision(6) << answer << "\n";
return 0;
}
Python
import math
import sys
EPS = 1e-10
class Simplex:
def __init__(self, m, n):
self.m = m
self.n = n
self.T = [[0.0] * (n + 1) for _ in range(m + 1)]
self.basis = [-1] * m
def pivot(self, r, c):
inv = 1.0 / self.T[r][c]
for j in range(self.n + 1):
self.T[r][j] *= inv
for i in range(self.m + 1):
if i == r:
continue
factor = self.T[i][c]
if abs(factor) <= EPS:
continue
for j in range(self.n + 1):
self.T[i][j] -= factor * self.T[r][j]
self.basis[r] = c
def set_objective(self, c_obj):
for j in range(self.n + 1):
self.T[self.m][j] = 0.0
for j in range(self.n):
self.T[self.m][j] = -c_obj[j]
for i in range(self.m):
var = self.basis[i]
if var < 0:
continue
coeff = c_obj[var]
if abs(coeff) <= EPS:
continue
for j in range(self.n + 1):
self.T[self.m][j] += coeff * self.T[i][j]
def solve(self):
while True:
enter = -1
for j in range(self.n):
if self.T[self.m][j] < -EPS:
enter = j
break
if enter == -1:
return True
min_ratio = 0.0
leave = -1
for i in range(self.m):
a = self.T[i][enter]
if a > EPS:
ratio = self.T[i][self.n] / a
if leave == -1 or ratio < min_ratio - 1e-12 or (abs(ratio - min_ratio) <= 1e-12 and self.basis[i] < self.basis[leave]):
min_ratio = ratio
leave = i
if leave == -1:
return False
self.pivot(leave, enter)
def objective_value(self):
return self.T[self.m][self.n]
class GameData:
def __init__(self):
self.dice = 0
self.actions = 0
self.states = 0
self.obs = 0
self.prob = 0.0
self.visible_max = []
self.action_obs = []
self.action_hidden = []
def build_game(dice):
g = GameData()
g.dice = dice
g.actions = dice
if dice == 2:
g.states = 36
g.obs = 6
g.prob = 1.0 / 36.0
g.visible_max = [v + 1 for v in range(g.obs)]
n_action = g.states * g.actions
g.action_obs = [0] * n_action
g.action_hidden = [0] * n_action
for a in range(1, 7):
for b in range(1, 7):
s = (a - 1) * 6 + (b - 1)
base = s * 2
g.action_obs[base] = b - 1
g.action_hidden[base] = a
g.action_obs[base + 1] = a - 1
g.action_hidden[base + 1] = b
return g
if dice == 3:
g.states = 216
g.obs = 21
g.prob = 1.0 / 216.0
obs_index = [[-1] * 7 for _ in range(7)]
g.visible_max = [0] * g.obs
idx = 0
for u in range(1, 7):
for v in range(u, 7):
obs_index[u][v] = idx
obs_index[v][u] = idx
g.visible_max[idx] = v
idx += 1
n_action = g.states * g.actions
g.action_obs = [0] * n_action
g.action_hidden = [0] * n_action
for a in range(1, 7):
for b in range(1, 7):
for c in range(1, 7):
s = ((a - 1) * 6 + (b - 1)) * 6 + (c - 1)
base = s * 3
g.action_obs[base] = obs_index[b][c]
g.action_hidden[base] = a
g.action_obs[base + 1] = obs_index[a][c]
g.action_hidden[base + 1] = b
g.action_obs[base + 2] = obs_index[a][b]
g.action_hidden[base + 2] = c
return g
return g
class LPData:
def __init__(self, rows, cols, actions, tvars, slack, art, art_s, slack_s):
self.lp = Simplex(rows, cols)
self.n_action = actions
self.n_T = tvars
self.n_slack = slack
self.n_art = art
self.art_start = art_s
self.slack_start = slack_s
def build_lp(g):
n_action = g.states * g.actions
n_T = g.obs
n_orig = n_action + n_T
n_slack = g.obs * 2
n_art = g.states
n_total = n_orig + n_slack + n_art
rows = g.states + n_slack
slack_start = n_orig
art_start = n_orig + n_slack
data = LPData(rows, n_total, n_action, n_T, n_slack, n_art, art_start, slack_start)
lp = data.lp
for s in range(g.states):
row = s
base = s * g.actions
for a in range(g.actions):
lp.T[row][base + a] = 1.0
lp.T[row][art_start + s] = 1.0
lp.T[row][n_total] = 1.0
lp.basis[row] = art_start + s
for o in range(g.obs):
for t in range(2):
row = g.states + o * 2 + t
tvar = n_action + o
lp.T[row][tvar] = -1.0
lp.T[row][slack_start + o * 2 + t] = 1.0
lp.T[row][n_total] = 0.0
lp.basis[row] = slack_start + o * 2 + t
for idx in range(n_action):
obs = g.action_obs[idx]
hidden = g.action_hidden[idx]
row_hidden = g.states + obs * 2
row_visible = row_hidden + 1
lp.T[row_hidden][idx] += g.prob * hidden
lp.T[row_visible][idx] += g.prob * g.visible_max[obs]
return data
def remove_artificial(lp, art_start, art_count):
if art_count == 0:
return
rows_to_remove = []
is_art = [False] * lp.n
for i in range(art_count):
is_art[art_start + i] = True
for i in range(lp.m):
var = lp.basis[i]
if var < 0 or not is_art[var]:
continue
pivot_col = -1
for j in range(lp.n):
if is_art[j]: continue
if abs(lp.T[i][j]) > EPS:
pivot_col = j
break
if pivot_col != -1:
lp.pivot(i, pivot_col)
else:
rows_to_remove.append(i)
if rows_to_remove:
drop = [False] * lp.m
for r in rows_to_remove:
drop[r] = True
new_m = lp.m - len(rows_to_remove)
newT = [[0.0] * (lp.n + 1) for _ in range(new_m + 1)]
new_basis = [-1] * new_m
r2 = 0
for i in range(lp.m):
if drop[i]: continue
newT[r2] = lp.T[i][:]
new_basis[r2] = lp.basis[i]
r2 += 1
newT[new_m] = lp.T[lp.m][:]
lp.T = newT
lp.basis = new_basis
lp.m = new_m
new_n = lp.n - art_count
newT2 = [[0.0] * (new_n + 1) for _ in range(lp.m + 1)]
for i in range(lp.m + 1):
col_new = 0
for j in range(lp.n):
if not is_art[j]:
newT2[i][col_new] = lp.T[i][j]
col_new += 1
newT2[i][new_n] = lp.T[i][lp.n]
lp.n = new_n
lp.T = newT2
def solve_game(dice):
g = build_game(dice)
if g.states == 0:
return float('nan')
data = build_lp(g)
lp = data.lp
c1 = [0.0] * lp.n
for i in range(data.n_art):
c1[data.art_start + i] = -1.0
lp.set_objective(c1)
if not lp.solve():
return float('nan')
phase1 = lp.objective_value()
if phase1 < -1e-8:
return float('nan')
remove_artificial(lp, data.art_start, data.n_art)
c2 = [0.0] * lp.n
for i in range(data.n_T):
c2[data.n_action + i] = -1.0
lp.set_objective(c2)
if not lp.solve():
return float('nan')
return -lp.objective_value()
def run_validation():
expected = 145.0 / 36.0
got = solve_game(2)
assert abs(got - expected) <= 1e-7
def solve():
answer = solve_game(3)
return f"{answer:.6f}"
if __name__ == "__main__":
run_validation()
print(solve())
Java
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.Locale;
public class Euler982 {
static final double EPS = 1e-10;
static class Simplex {
int m, n;
double[][] T;
int[] basis;
Simplex(int rows, int cols) {
m = rows;
n = cols;
T = new double[m + 1][n + 1];
basis = new int[m];
Arrays.fill(basis, -1);
}
void pivot(int r, int c) {
double inv = 1.0 / T[r][c];
for (int j = 0; j <= n; ++j)
T[r][j] *= inv;
for (int i = 0; i <= m; ++i) {
if (i == r)
continue;
double factor = T[i][c];
if (Math.abs(factor) <= EPS)
continue;
for (int j = 0; j <= n; ++j)
T[i][j] -= factor * T[r][j];
}
basis[r] = c;
}
void setObjective(double[] cObj) {
Arrays.fill(T[m], 0.0);
for (int j = 0; j < n; ++j)
T[m][j] = -cObj[j];
T[m][n] = 0.0;
for (int i = 0; i < m; ++i) {
int var = basis[i];
if (var < 0)
continue;
double coeff = cObj[var];
if (Math.abs(coeff) <= EPS)
continue;
for (int j = 0; j <= n; ++j)
T[m][j] += coeff * T[i][j];
}
}
boolean solve() {
while (true) {
int enter = -1;
for (int j = 0; j < n; ++j) {
if (T[m][j] < -EPS) {
enter = j;
break;
}
}
if (enter == -1)
return true;
double minRatio = 0.0;
int leave = -1;
for (int i = 0; i < m; ++i) {
double a = T[i][enter];
if (a > EPS) {
double ratio = T[i][n] / a;
if (leave == -1 || ratio < minRatio - 1e-12 ||
(Math.abs(ratio - minRatio) <= 1e-12 && basis[i] < basis[leave])) {
minRatio = ratio;
leave = i;
}
}
}
if (leave == -1)
return false;
pivot(leave, enter);
}
}
double objectiveValue() {
return T[m][n];
}
}
static class GameData {
int dice, actions, states, obs;
double prob;
int[] visibleMax, actionObs, actionHidden;
}
static GameData buildGame(int dice) {
GameData g = new GameData();
g.dice = dice;
g.actions = dice;
if (dice == 2) {
g.states = 36;
g.obs = 6;
g.prob = 1.0 / 36.0;
g.visibleMax = new int[g.obs];
for (int v = 0; v < g.obs; ++v)
g.visibleMax[v] = v + 1;
int nAction = g.states * g.actions;
g.actionObs = new int[nAction];
g.actionHidden = new int[nAction];
for (int a = 1; a <= 6; ++a) {
for (int b = 1; b <= 6; ++b) {
int s = (a - 1) * 6 + (b - 1);
int base = s * 2;
g.actionObs[base] = b - 1;
g.actionHidden[base] = a;
g.actionObs[base + 1] = a - 1;
g.actionHidden[base + 1] = b;
}
}
return g;
}
if (dice == 3) {
g.states = 216;
g.obs = 21;
g.prob = 1.0 / 216.0;
int[][] obsIndex = new int[7][7];
for (int i = 0; i <= 6; ++i)
Arrays.fill(obsIndex[i], -1);
g.visibleMax = new int[g.obs];
int idx = 0;
for (int u = 1; u <= 6; ++u) {
for (int v = u; v <= 6; ++v) {
obsIndex[u][v] = idx;
obsIndex[v][u] = idx;
g.visibleMax[idx] = v;
idx++;
}
}
int nAction = g.states * g.actions;
g.actionObs = new int[nAction];
g.actionHidden = new int[nAction];
for (int a = 1; a <= 6; ++a) {
for (int b = 1; b <= 6; ++b) {
for (int c = 1; c <= 6; ++c) {
int s = ((a - 1) * 6 + (b - 1)) * 6 + (c - 1);
int base = s * 3;
g.actionObs[base] = obsIndex[b][c];
g.actionHidden[base] = a;
g.actionObs[base + 1] = obsIndex[a][c];
g.actionHidden[base + 1] = b;
g.actionObs[base + 2] = obsIndex[a][b];
g.actionHidden[base + 2] = c;
}
}
}
return g;
}
return g;
}
static class LPData {
Simplex lp;
int nAction, nT, nSlack, nArt, artStart, slackStart;
LPData(int rows, int cols, int actions, int tvars, int slack, int art, int artS, int slackS) {
lp = new Simplex(rows, cols);
nAction = actions;
nT = tvars;
nSlack = slack;
nArt = art;
artStart = artS;
slackStart = slackS;
}
}
static LPData buildLP(GameData g) {
int nAction = g.states * g.actions;
int nT = g.obs;
int nOrig = nAction + nT;
int nSlack = g.obs * 2;
int nArt = g.states;
int nTotal = nOrig + nSlack + nArt;
int rows = g.states + nSlack;
int slackStart = nOrig;
int artStart = nOrig + nSlack;
LPData data = new LPData(rows, nTotal, nAction, nT, nSlack, nArt, artStart, slackStart);
Simplex lp = data.lp;
for (int s = 0; s < g.states; ++s) {
int row = s;
int base = s * g.actions;
for (int a = 0; a < g.actions; ++a)
lp.T[row][base + a] = 1.0;
lp.T[row][artStart + s] = 1.0;
lp.T[row][nTotal] = 1.0;
lp.basis[row] = artStart + s;
}
for (int o = 0; o < g.obs; ++o) {
for (int t = 0; t < 2; ++t) {
int row = g.states + o * 2 + t;
int tvar = nAction + o;
lp.T[row][tvar] = -1.0;
lp.T[row][slackStart + o * 2 + t] = 1.0;
lp.T[row][nTotal] = 0.0;
lp.basis[row] = slackStart + o * 2 + t;
}
}
for (int idx = 0; idx < nAction; ++idx) {
int obs = g.actionObs[idx];
int hidden = g.actionHidden[idx];
int rowHidden = g.states + obs * 2;
int rowVisible = rowHidden + 1;
lp.T[rowHidden][idx] += g.prob * hidden;
lp.T[rowVisible][idx] += g.prob * g.visibleMax[obs];
}
return data;
}
static void removeArtificial(Simplex lp, int artStart, int artCount) {
if (artCount == 0)
return;
List<Integer> rowsToRemove = new ArrayList<>();
boolean[] isArt = new boolean[lp.n];
for (int i = 0; i < artCount; ++i)
isArt[artStart + i] = true;
for (int i = 0; i < lp.m; ++i) {
int var = lp.basis[i];
if (var < 0 || !isArt[var])
continue;
int pivotCol = -1;
for (int j = 0; j < lp.n; ++j) {
if (isArt[j])
continue;
if (Math.abs(lp.T[i][j]) > EPS) {
pivotCol = j;
break;
}
}
if (pivotCol != -1) {
lp.pivot(i, pivotCol);
} else {
rowsToRemove.add(i);
}
}
if (!rowsToRemove.isEmpty()) {
boolean[] drop = new boolean[lp.m];
for (int r : rowsToRemove)
drop[r] = true;
int newM = lp.m - rowsToRemove.size();
double[][] newT = new double[newM + 1][lp.n + 1];
int[] newBasis = new int[newM];
int r2 = 0;
for (int i = 0; i < lp.m; ++i) {
if (drop[i])
continue;
System.arraycopy(lp.T[i], 0, newT[r2], 0, lp.n + 1);
newBasis[r2] = lp.basis[i];
r2++;
}
System.arraycopy(lp.T[lp.m], 0, newT[newM], 0, lp.n + 1);
lp.T = newT;
lp.basis = newBasis;
lp.m = newM;
}
int newN = lp.n - artCount;
double[][] newT = new double[lp.m + 1][newN + 1];
for (int i = 0; i <= lp.m; ++i) {
int colNew = 0;
for (int j = 0; j < lp.n; ++j) {
if (isArt[j])
continue;
newT[i][colNew++] = lp.T[i][j];
}
newT[i][newN] = lp.T[i][lp.n];
}
lp.n = newN;
lp.T = newT;
}
static double solveGame(int dice) {
GameData g = buildGame(dice);
if (g.states == 0)
return Double.NaN;
LPData data = buildLP(g);
Simplex lp = data.lp;
double[] c1 = new double[lp.n];
for (int i = 0; i < data.nArt; ++i)
c1[data.artStart + i] = -1.0;
lp.setObjective(c1);
if (!lp.solve())
return Double.NaN;
double phase1 = lp.objectiveValue();
if (phase1 < -1e-8)
return Double.NaN;
removeArtificial(lp, data.artStart, data.nArt);
double[] c2 = new double[lp.n];
for (int i = 0; i < data.nT; ++i)
c2[data.nAction + i] = -1.0;
lp.setObjective(c2);
if (!lp.solve())
return Double.NaN;
return -lp.objectiveValue();
}
public static String solve() {
double res = solveGame(3);
return String.format(Locale.US, "%.6f", res);
}
public static void main(String[] args) {
double expected = 145.0 / 36.0;
double got = solveGame(2);
if (Math.abs(got - expected) > 1e-7) {
System.err.println("Validation failed");
return;
}
System.out.println(solve());
}
}