Project Euler Problem 389: Platonic Dice

Mathematical Analysis

Single Die Properties

For a single fair n-sided die (faces 1 through n):

• Expected value: E[X] = (n + 1) / 2
• Variance: Var(X) = (n^2 - 1) / 12

Law of Total Variance

The key identity used at each layer is the law of total variance:

Var(X) = E[Var(X | Y)] + Var(E[X | Y])

If X is the sum of Y independent n-sided dice, then:

• E[X | Y] = Y * (n + 1) / 2
• Var(X | Y) = Y * (n^2 - 1) / 12

Substituting into the law of total variance:

• E[X] = E[Y] * (n + 1) / 2
• Var(X) = E[Y] * (n^2 - 1) / 12 + Var(Y) * ((n + 1) / 2)^2

Derivation

We apply the recurrence layer by layer. Let E_k and V_k denote the mean and variance after layer k.

Layer 0: T ~ d4

Layer 1: C = sum of T d6 dice

• E[C] = E[T] * (7/2) = 2.5 * 3.5 = 8.75
• Var(C) = E[T] * (35/12) + Var(T) * (7/2)^2
• Var(C) = 2.5 * 2.91667 + 1.25 * 12.25 = 7.29167 + 15.3125 = 22.6042

Layer 2: O = sum of C d8 dice

• E[O] = 8.75 * 4.5 = 39.375
• Var(O) = 8.75 * (63/12) + 22.6042 * (4.5)^2
• Var(O) = 8.75 * 5.25 + 22.6042 * 20.25 = 45.9375 + 457.7344 = 503.6719

Layer 3: D = sum of O d12 dice

• E[D] = 39.375 * 6.5 = 255.9375
• Var(D) = 39.375 * (143/12) + 503.6719 * (6.5)^2
• Var(D) = 39.375 * 11.91667 + 503.6719 * 42.25 = 469.2188 + 21280.0365 = 21749.2552

Layer 4: I = sum of D d20 dice

• E[I] = 255.9375 * 10.5 = 2687.3438
• Var(I) = 255.9375 * (399/12) + 21749.2552 * (10.5)^2
• Var(I) = 255.9375 * 33.25 + 21749.2552 * 110.25
• Var(I) = 8509.9219 + 2397866.4404 = 2406376.3623

Proof of Correctness

The derivation rests on two standard results:

1. Law of total expectation: E[X] = E[E[X|Y]] — used to propagate means.
2. Law of total variance: Var(X) = E[Var(X|Y)] + Var(E[X|Y]) — used to propagate variances.

Both hold for any random variables X, Y with finite second moments. Since each layer is a sum of i.i.d. dice conditioned on the count from the previous layer, the conditional mean and variance are linear in the count, making the formulas exact.

Correctness

Theorem. The method described above computes exactly the quantity requested in the problem statement.

Proof. The preceding analysis identifies the admissible objects and derives the formula, recurrence, or exhaustive search carried out by the algorithm. The computation evaluates exactly that specification, so every valid contribution is included once and no invalid contribution is counted. Therefore the returned value is the required answer. $\square$

Complexity Analysis

• Time: O(1) — five layers of constant-time arithmetic.
• Space: O(1) — only two running values (mean and variance) are maintained.

No simulation, enumeration, or generating functions are needed. The recursive variance formula gives the exact answer in closed form.

Answer

#include <bits/stdc++.h>
using namespace std;

/*
 * Project Euler Problem 389: Platonic Dice
 *
 * Chain: T~d4 -> C=sum(T x d6) -> O=sum(C x d8) -> D=sum(O x d12) -> I=sum(D x d20)
 *
 * Law of total variance applied recursively:
 *   E[sum]   = E[Y] * (n+1)/2
 *   Var[sum] = E[Y] * (n^2-1)/12 + Var[Y] * ((n+1)/2)^2
 *
 * Answer: 2406376.3623
 */

int main() {
    // die_mean(n) = (n+1)/2,  die_var(n) = (n*n - 1)/12
    int sides[] = {4, 6, 8, 12, 20};

    // Layer 0: single d4
    double e = (sides[0] + 1.0) / 2.0;
    double v = (sides[0] * sides[0] - 1.0) / 12.0;

    // Layers 1 through 4
    for (int i = 1; i < 5; i++) {
        int n = sides[i];
        double mu    = (n + 1.0) / 2.0;
        double sigma2 = (1.0 * n * n - 1.0) / 12.0;

        double new_e = e * mu;
        double new_v = e * sigma2 + v * mu * mu;

        e = new_e;
        v = new_v;
    }

    // Round to 4 decimal places
    double answer = round(v * 10000.0) / 10000.0;

    cout << fixed << setprecision(4);
    cout << "Var(I) = " << answer << endl;

    // Verify
    assert(abs(answer - 2406376.3623) < 1e-4);
    cout << "Assertion passed." << endl;

    return 0;
}

"""
Project Euler Problem 389: Platonic Dice
=========================================

Find the variance of I, where the chain is:
    T ~ d4
    C = sum of T d6
    O = sum of C d8
    D = sum of O d12
    I = sum of D d20

Uses the law of total variance applied recursively:
    E[sum]   = E[Y] * (n+1)/2
    Var[sum] = E[Y] * (n^2-1)/12 + Var[Y] * ((n+1)/2)^2

Answer: 2406376.3623
"""

from __future__ import annotations

def die_mean(n: int) -> float:
    """Expected value of a single fair n-sided die."""
    return (n + 1) / 2

def die_var(n: int) -> float:
    """Variance of a single fair n-sided die."""
    return (n ** 2 - 1) / 12

def propagate(e_y: float, v_y: float, n: int):
    """
    Given Y with mean e_y and variance v_y, compute the mean and variance
    of the sum of Y independent fair n-sided dice.

    Parameters
    ----------
    e_y : float
        E[Y], the expected number of dice.
    v_y : float
        Var(Y), the variance of the number of dice.
    n : int
        Number of sides on each die.

    Returns
    -------
    (E[sum], Var(sum))
    """
    mu = die_mean(n)
    sigma2 = die_var(n)
    e_sum = e_y * mu
    v_sum = e_y * sigma2 + v_y * mu ** 2
    return e_sum, v_sum

def solve() -> float:
    """Compute Var(I) through the five Platonic dice layers."""
    dice_sides = [4, 6, 8, 12, 20]

    # Layer 0: single d4
    e = die_mean(dice_sides[0])
    v = die_var(dice_sides[0])

    # Layers 1-4: each layer rolls previous-sum dice of the next type
    for n in dice_sides[1:]:
        e, v = propagate(e, v, n)

    return round(v, 4)

def collect_layer_data():
    """Collect mean and variance at each layer for visualization."""
    dice_sides = [4, 6, 8, 12, 20]
    labels = ["T (d4)", "C (d6)", "O (d8)", "D (d12)", "I (d20)"]

    e = die_mean(dice_sides[0])
    v = die_var(dice_sides[0])
    layers = [{"label": labels[0], "mean": e, "variance": v}]

    for i, n in enumerate(dice_sides[1:], start=1):
        e, v = propagate(e, v, n)
        layers.append({"label": labels[i], "mean": e, "variance": v})

    return layers

def visualize(output_path: str) -> None:
    """Create a visualization of how mean and variance grow across layers."""
    layers = collect_layer_data()
    labels = [d["label"] for d in layers]
    means = [d["mean"] for d in layers]
    variances = [d["variance"] for d in layers]