Problem 389: Platonic Dice

View on Project Euler

Project Euler Problem 389 Solution

EulerSolve provides an optimized solution for Project Euler Problem 389, Platonic Dice, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The experiment is a chain of random sums built from the five Platonic dice. First roll a 4-sided die and call the result \(T\). Then roll \(T\) fair 6-sided dice and sum them to obtain \(C\). Next roll \(C\) fair 8-sided dice to obtain \(O\), then \(O\) fair 12-sided dice to obtain \(D\), and finally \(D\) fair 20-sided dice to obtain \(I\). The task is to compute \(\operatorname{Var}(I)\). Mathematical Approach The local solution files do not enumerate huge distributions. Instead, they propagate only the mean and the variance from one layer to the next. That works because every stage has the same structure: a random number of independent fair dice are rolled and summed. Step 1: Moments of One Fair Die For a fair \(s\)-sided die \(X_s\) with values \(1,2,\dots,s\), the code uses $$\mathbb{E}[X_s]=\frac{s+1}{2},\qquad \operatorname{Var}(X_s)=\frac{s^2-1}{12}.$$ These are exactly the values returned by uniform_die_stats in C++, Python, and Java. For the first die, $$T \sim \text{Unif}\{1,2,3,4\},\qquad \mathbb{E}[T]=\frac{5}{2},\qquad \operatorname{Var}(T)=\frac{5}{4}.$$ Step 2: Random-Sum Identity Suppose \(N\) is a nonnegative integer-valued random variable, and conditional on \(N\) we roll \(N\) independent copies \(X_1,\dots,X_N\) of the same fair die, each with mean \(\mu\) and variance \(\sigma^2\)....

Detailed mathematical approach

Problem Summary

An unbiased \(4\)-sided die is thrown; its value is \(T\). Then \(T\) unbiased \(6\)-sided dice are thrown and their scores added together; the sum is \(C\). In the same way \(C\) unbiased \(8\)-sided dice give the sum \(O\), then \(O\) \(12\)-sided dice give \(D\), and finally \(D\) \(20\)-sided dice give \(I\). We must find the variance of \(I\), rounded to four decimal places. [1] The face counts are those of the five Platonic solids, and the letters are the initials of tetrahedron, cube, octahedron, dodecahedron and icosahedron. [2] The statement gives no checkpoint values, so the implementation must be measured against independent mathematics.

The whole solution in one sentence: the mean and the variance of each total can be computed from the mean and the variance of the previous total alone, so five short updates replace any attempt to describe the distribution of \(I\).

Mathematical Approach

1. The experiment written out

Write \(s_1=4\), \(s_2=6\), \(s_3=8\), \(s_4=12\), \(s_5=20\) for the face counts in the order they are used, put \(X_0=1\), and let \(X_k\) be the sum of \(X_{k-1}\) dice with \(s_k\) faces, so that \(X_1=T\), \(X_2=C\), \(X_3=O\), \(X_4=D\) and \(X_5=I\). The word unbiased fixes the model: inside a layer the dice are fair and mutually independent, and independent of all earlier layers, which decide only how many of them are added up.

Every die shows at least \(1\), so \(X_k\ge X_{k-1}\ge1\) and no layer ever asks for zero dice. The largest value of \(I\) is \(4\cdot6\cdot8\cdot12\cdot20=46080\) and the smallest is \(1\), while the number of throw sequences is about \(10^{3002}\). Even the distribution of \(I\), a list of \(46080\) probabilities, is more than the calculation needs: the recursion below carries only a mean and a variance from one layer to the next.

2. One fair die

The mean of \(X\) is the probability-weighted average \(\mathbb{E}[X]=\sum_x xp_x\), with \(p_x=\mathbb{P}(X=x)\), and the variance is the weighted average of the squared distances from it, \(\operatorname{Var}(X)=\sum_x(x-\mathbb{E}[X])^2p_x\), which expands to \(\mathbb{E}[X^2]-\mathbb{E}[X]^2\). A fair two-sided die shows \(1\) or \(2\): mean \(3/2\), both values \(1/2\) away, variance \(1/4\).

Now let \(X\) be one fair \(s\)-sided die's score, uniform on \(\{1,2,\dots,s\}\). [3] From \(\sum_{k=1}^{s}k=\tfrac{s(s+1)}2\) and \(\sum_{k=1}^{s}k^2=\tfrac{s(s+1)(2s+1)}6\),

$$\mu_s=\mathbb{E}[X]=\frac{s+1}{2},$$

$$\sigma_s^2=\mathbb{E}[X^2]-\mu_s^2=\frac{s^2-1}{12}.$$

For our five dice:

  • \(s=4\): \(\mu=5/2\), \(\sigma^2=5/4\);
  • \(s=6\): \(\mu=7/2\), \(\sigma^2=35/12\);
  • \(s=8\): \(\mu=9/2\), \(\sigma^2=21/4\);
  • \(s=12\): \(\mu=13/2\), \(\sigma^2=143/12\);
  • \(s=20\): \(\mu=21/2\), \(\sigma^2=133/4\).

3. Conditioning on what has already been thrown

Conditioning is the tool: once \(N\) is given, the layer is an ordinary sum of a fixed number of dice. Two laws turn conditional knowledge back into plain numbers. The law of total expectation [4] says that averaging the conditional averages gives the average,

$$\mathbb{E}[Y]=\mathbb{E}\bigl[\mathbb{E}[Y\mid N]\bigr],$$

and the law of total variance [5] splits the spread of \(Y\) into the average spread inside each group plus the spread of the group averages,

$$\operatorname{Var}(Y)=\mathbb{E}\bigl[\operatorname{Var}(Y\mid N)\bigr]+\operatorname{Var}\bigl(\mathbb{E}[Y\mid N]\bigr).$$

The second law follows from the first. Write \(A=\mathbb{E}[Y\mid N]\), itself random because \(N\) is. Applying \(\operatorname{Var}(Z)=\mathbb{E}[Z^2]-\mathbb{E}[Z]^2\) conditionally gives \(\mathbb{E}[Y^2\mid N]=\operatorname{Var}(Y\mid N)+A^2\). Take expectations and subtract \(\mathbb{E}[Y]^2=\mathbb{E}[A]^2\): the leftover terms are \(\mathbb{E}[\operatorname{Var}(Y\mid N)]\) and \(\operatorname{Var}(A)\). Both parts are needed, and dropping the second is the classic mistake here.

4. The variance of a sum of randomly many dice

Let \(N\) be a random count, let \(X_1,X_2,\dots\) be independent \(s\)-sided dice, independent of \(N\), and let \(S=X_1+\dots+X_N\). Given \(N=n\) it is a sum of \(n\) independent dice, so the conditional mean is \(n\mu_s\) and, variances of independent terms adding up, the conditional variance is \(n\sigma_s^2\). Section 3 then gives

$$\mathbb{E}[S]=\mu_s\,\mathbb{E}[N],$$

$$\operatorname{Var}(S)=\sigma_s^2\,\mathbb{E}[N]+\mu_s^2\operatorname{Var}(N).$$

The second identity is the Blackwell–Girshick equation. [6] Nothing else about \(N\) enters: its mean and variance are enough, which is why the layers can be handled one after another.

A tiny example can be checked by hand. Take two-sided dice and let \(N\) be such a die as well. If \(N=1\) the sum is \(1\) or \(2\); if \(N=2\) it is \(2\), \(3\) or \(4\) with probabilities \(\tfrac14,\tfrac12,\tfrac14\). Weighting both cases by \(\tfrac12\) and writing \(p_j=\mathbb{P}(S=j)\),

$$(p_1,p_2,p_3,p_4)=(\tfrac14,\tfrac38,\tfrac14,\tfrac18),$$

so \(\mathbb{E}[S]=\tfrac94\), \(\mathbb{E}[S^2]=6\) and \(\operatorname{Var}(S)=6-\tfrac{81}{16}=\tfrac{15}{16}\). The formulas give \(\tfrac32\cdot\tfrac32=\tfrac94\) and \(\tfrac14\cdot\tfrac32+\tfrac94\cdot\tfrac14=\tfrac{15}{16}\), as they must.

5. The chain as a two-number update

Section 4 turns the problem into a recursion on pairs. If the current total has mean \(m\) and variance \(v\), then after a layer of \(s\)-sided dice the new total has

$$m'=\mu_s\,m,\qquad v'=\sigma_s^2\,m+\mu_s^2\,v.$$

We start at \(X_0=1\), which is not random: \(m=1\), \(v=0\). The first layer then reproduces the single \(4\)-sided die, \(\mathbb{E}[T]=5/2\) and \(\operatorname{Var}(T)=5/4\). The second layer, with \(\mu_6=7/2\) and \(\sigma_6^2=35/12\), gives

$$\mathbb{E}[C]=\frac72\cdot\frac52=\frac{35}{4},$$

$$\operatorname{Var}(C)=\frac{35}{12}\cdot\frac52+\frac{49}{4}\cdot\frac54=\frac{1085}{48}.$$

This layer is small enough to confirm by brute force. For each \(T=n\), enumerate the \(6^n\) ordered score sequences, each of probability \(1/(4\cdot6^n)\): there are \(6+36+216+1296=1554\) of them, with different weights for different \(n\). Summing their weighted scores and squared deviations reproduces \(35/4\) and \(1085/48\). Two probabilities are easy to check by hand: \(\mathbb{P}(C=1)=\tfrac14\cdot\tfrac16=\tfrac1{24}\), and the maximum \(C=24\) needs \(T=4\) and four sixes, so \(\mathbb{P}(C=24)=\tfrac14\cdot6^{-4}=\tfrac1{5184}\).

6. Why both terms are needed

It is tempting to replace the random count by its mean and treat \(C\) as a sum of \(2.5\) dice. That gives \(\sigma_6^2\mathbb{E}[T]=175/24\), the first term alone. The second term, \(\mu_6^2\operatorname{Var}(T)=245/16\), is the larger of the two, and the correct variance \(1085/48\) is exactly \(31/10\) times the naive value: most of the spread of \(C\) comes from not knowing how many dice will be thrown, not from the scores themselves. The effect compounds in later layers, where the count is itself random.

7. The five layers

Applying the update again, with \(\mu_8=9/2\) and \(\sigma_8^2=21/4\),

$$\mathbb{E}[O]=\frac92\cdot\frac{35}{4}=\frac{315}{8},$$

$$\operatorname{Var}(O)=\frac{21}{4}\cdot\frac{35}{4}+\frac{81}{4}\cdot\frac{1085}{48}=\frac{32235}{64}.$$

Two more rounds finish the computation, which the program performs. The means simply multiply up, since \(m'=\mu_s m\) has no other contribution: after \(k\) layers,

$$\mathbb{E}[X_k]=\prod_{j=1}^{k}\frac{s_j+1}{2}.$$

That product cannot depend on the order of the layers. The variance update mixes the two quantities and has no such symmetry at first sight — yet it is order-independent too, for the reason that follows.

8. A hidden identity

In every pair computed so far, \((5/2,\,5/4)\), \((35/4,\,1085/48)\) and \((315/8,\,32235/64)\), the variance is determined by the mean through one formula,

$$\operatorname{Var}(Y)=\frac{\mathbb{E}[Y]\bigl(\mathbb{E}[Y]-1\bigr)}{3}.$$

Check the first, \(\tfrac52\cdot\tfrac32/3=\tfrac54\), and the third, \(\tfrac{315}8\cdot\tfrac{307}8/3=\tfrac{32235}{64}\). A single fair die obeys the same rule, because \(s=2\mu_s-1\) gives \(s^2-1=4\mu_s(\mu_s-1)\) and therefore

$$\sigma_s^2=\frac{s^2-1}{12}=\frac{\mu_s(\mu_s-1)}{3}.$$

So the identity holds at the start of the chain, and every layer preserves it. Suppose the count \(N\) satisfies it, write \(m=\mathbb{E}[N]\) and \(\mu=\mu_s\), and substitute both special forms into the update of Section 5:

$$\begin{aligned}3\operatorname{Var}(S)&=\mu(\mu-1)m+\mu^2m(m-1)\\&=\mu m\bigl(\mu-1+\mu m-\mu\bigr)\\&=\mu m(\mu m-1).\end{aligned}$$

Since \(\mathbb{E}[S]=\mu m\), the right-hand side is \(\mathbb{E}[S](\mathbb{E}[S]-1)\), which is the identity for \(S\). By induction every total in the chain satisfies it, the last one included: each layer's variance is a function of its mean alone, an algebraic check on the recursion. An equivalent form, \(\mathbb{E}[Y(Y-1)]=\tfrac43\mathbb{E}[Y](\mathbb{E}[Y]-1)\), makes the next section short.

9. Consequences and a second viewpoint

Since the mean is the product of the single-die means and the variance is a function of the mean, both survive any reordering of the dice: throwing the \(20\)-sided dice first and the \(4\)-sided die last gives the same answer. The proof covers all \(120\) orders, and computing them confirms a single pair. For a large mean \(\operatorname{Var}(Y)\approx\mathbb{E}[Y]^2/3\), so the standard deviation divided by the mean approaches \(1/\sqrt3\approx0.577\) as the mean grows.

Generating functions explain both facts at once. The generating function of a distribution is the polynomial \(G_Y(z)=\sum_j\mathbb{P}(Y=j)z^j\); multiplying two of them adds independent scores, since powers add while probabilities multiply. A sum of \(n\) dice with \(s\) faces thus has \(g_s(z)^n\) for \(g_s(z)=(z+\dots+z^{s})/s\), and averaging over a random count gives the composition \(G_N(g_s(z))\), [7] so \(I\) has \(g_4\circ g_6\circ g_8\circ g_{12}\circ g_{20}\). Differentiating \(z^j\) once or twice at \(z=1\) leaves \(j\) or \(j(j-1)\), so \(G'(1)=\mathbb{E}[Y]\) and \(G''(1)=\mathbb{E}[Y(Y-1)]\). The chain rule then makes the mean of a composition the product of the single-die derivatives, hence order-free, and Section 8 says every factor satisfies \(g''(1)=\tfrac43g'(1)(g'(1)-1)\), which composition preserves.

10. Rounding to four places

The exact value is a rational number whose denominator is a power of two, so its decimal expansion terminates and it lies more than \(10^{-5}\) from the nearest halfway point of the fourth decimal. A chain computed entirely in double arithmetic differs from the exact value by less than \(10^{-8}\), and casting a long double result to double adds at most one unit in its last place, so that margin is never approached and the four printed decimals are correct.

How the Code Works

The C++ program [8] keeps the mean and variance in a Stats struct of two long double fields. uniform_die_stats(sides) returns the pair \(\bigl((s+1)/2,\,(s^2-1)/12\bigr)\) of Section 2, and transition_sum_of_random_dice is the update of Section 5: it takes the statistics of the number of dice and a face count, fetches those of one such die, and returns \(\mu_s m\) and \(\sigma_s^2m+\mu_s^2v\). solve_variance starts from the \(4\)-sided die and calls it four times, for \(6\), \(8\), \(12\) and \(20\) faces, returning the last pair's variance. That is the entire computation.

The rest of the file checks itself. brute_one_layer_uniform_count builds the distribution of Section 5 honestly: for each count \(n\) up to its max_count it convolves the single-die distribution \(n\) times, spreading every probability over the sides faces, and adds the result into total_dist with weight \(1/\)max_count. The mean \(m\) comes from one pass over it and the variance from a second, as \(\sum_x p(x)(x-m)^2\), the stable form avoiding the cancellation in \(\mathbb{E}[X^2]-\mathbb{E}[X]^2\). run_checkpoints compares uniform_die_stats(4) with \(2.5\) and \(1.25\), then the update's pair for \(C\) with the brute-force pair for \((4,6)\), through nearly_equal, a relative test with tolerance \(10^{-14}\). It exits with code \(2\) on failure and \(1\) on an unknown argument; parse_arguments reads --skip-checkpoints.

Output goes through std::fixed and std::setprecision(4) after a cast of the long double result to double, so the printed digits come from double precision even where long double is wider, a step Section 10 covers. Python and Java implement the same two formulas and nothing else: plain tuples in Python's transition_sum, with solve looping over \([6,8,12,20]\), and a small Stats class with uniformDieStats and transitionSum in Java, both in ordinary double arithmetic and without the C++ checkpoint. Formatting differs: Java's String.format follows the default locale and the replace after it turns ASCII commas into points, which covers decimal-comma locales but not those formatting with another separator or with non-ASCII digits, as Arabic and Persian ones do; an explicit Locale.ROOT would. Java also rounds halfway cases up where C++ and Python round to even, immaterial so far from a tie.

Complexity and Verification

The solver computes one die's statistics five times (an addition, a subtraction, a multiplication and two divisions) and performs four transitions (four multiplications and one addition): forty-five arithmetic operations in the source, over a constant number of values. Time and memory are \(O(1)\) here, \(O(k)\) time and \(O(1)\) memory for \(k\) layers, whatever the face counts; nothing depends on \(46080\). The checkpoint restarts its convolution for every count, so for a maximum count \(M\) and \(s\) faces it costs \(\Theta(M^3s^2)\) time and \(O(Ms)\) memory: \(420\) inner steps on vectors of at most \(25\) entries for \((4,6)\). All three programs finish in well under a tenth of a second.

The checks below were run for this article, in exact rational arithmetic where fractions appear.

  • The per-die formulas agree with a direct sum over the faces for \(s=2,3,4,5,6,8,10,12,20,100\), and the two-sided example of Section 4 with a hand computation.
  • \(\mathbb{E}[C]\) and \(\operatorname{Var}(C)\) come out the same from the recursion as from all \(1554\) weighted outcomes of the first two layers, as do the pairs for \(T\) and \(O\) from their exact distributions.
  • The strongest check avoids the recursion: convolving the layers step by step gives the distribution of \(I\) over its \(46080\) values, about \(10^9\) multiply–add steps, whose mean and variance match the recursion and Section 8 to a relative error of \(3\cdot10^{-15}\).
  • A simulation of \(200\,000\) complete runs gave a sample variance within \(0.15\%\) of the exact value, inside the \(0.29\%\) standard error of such an estimate.
  • Section 8 holds for \(400\) random chains of up to seven dice with \(1\) to \(30\) faces, and for all \(120\) orders of the five Platonic dice.
  • The C++ program under GCC 16 and Clang 22, the Python version and the Java version print the same four-decimal string, which is the exact rational value rounded to four places.

These extra checks are not part of the published solutions. One portability remark: the C++ file calls std::max without including algorithm; both tested standard libraries supply the declaration through the file's other headers, so it compiles without diagnostics, though the standard does not guarantee that.

Footnotes and References

The references below supply the background facts used above. The two-number recursion, the identity of Section 8 and the independence of the order of the dice are derived explicitly in this article.

  1. Project Euler 389 — Platonic Dice. The official statement throws one unbiased 4-sided die for T, then T 6-sided dice for the sum C, C 8-sided dice for O, O 12-sided dice for D and D 20-sided dice for I, and asks for the variance of I rounded to 4 decimal places. It gives no intermediate values.
  2. Wikipedia — Platonic solid. It lists the five solids with four, six, eight, twelve and twenty faces, states that there are only five such polyhedra, and notes that they are often used as dice because these shapes can be made fair.
  3. Wikipedia — Discrete uniform distribution. For the integers from a to b it gives the mean (a + b)/2 and the variance ((b − a + 1)² − 1)/12, which for a = 1, b = s are the values of Section 2.
  4. Wikipedia — Law of total expectation. The identity E[X] = E[E[X | Y]], also called Adam's law or the tower rule, is stated and proved there; Section 3 applies it to a layer's conditional means.
  5. Wikipedia — Law of total variance. It states Var(Y) = E[Var(Y | X)] + Var(E[Y | X]), also known as Eve's law, and proves it from Var(Y) = E[Y²] − E[Y]² with the law of total expectation, as sketched in Section 3.
  6. Wikipedia — Blackwell–Girshick equation. For Y = X₁ + … + X_N with the X_i independent, identically distributed, independent of the non-negative integer count N and with second moments, it states and proves Var(Y) = Var(N)E(X₁)² + E(N)Var(X₁), the identity of Section 4.
  7. Wikipedia — Probability-generating function. It gives the generating function of a sum of N independent, identically distributed variables as the composition G_N(G_X(z)), and the mean and variance as G′(1) and G″(1) + G′(1) − G′(1)².
  8. C++ — Euler389.cpp. The linked immutable revision contains the two formulas, solve_variance and the checkpoint built on brute_one_layer_uniform_count, and is the implementation record for this article; the additional checks above are not packaged in it.

Problem 389 source code

C++

#include <cmath>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>

namespace {

struct Options {
    bool run_checkpoints = true;
};

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return true;
}

struct Stats {
    long double mean{0.0L};
    long double variance{0.0L};
};

Stats uniform_die_stats(const int sides) {
    const long double s = static_cast<long double>(sides);
    return Stats{(s + 1.0L) / 2.0L, (s * s - 1.0L) / 12.0L};
}

Stats transition_sum_of_random_dice(const Stats& count_stats, const int sides) {
    const Stats per_die = uniform_die_stats(sides);
    const long double mean =
        per_die.mean * count_stats.mean;
    const long double variance =
        per_die.variance * count_stats.mean +
        per_die.mean * per_die.mean * count_stats.variance;
    return Stats{mean, variance};
}

Stats brute_one_layer_uniform_count(const int max_count, const int sides) {
    // N is uniform on {1..max_count}; then S is sum of N fair `sides`-sided dice.
    const long double p_count = 1.0L / static_cast<long double>(max_count);
    const int max_sum = max_count * sides;
    std::vector<long double> total_dist(static_cast<std::size_t>(max_sum + 1), 0.0L);

    for (int n = 1; n <= max_count; ++n) {
        std::vector<long double> dist(1, 1.0L);  // sum=0 before rolling.
        for (int roll = 0; roll < n; ++roll) {
            std::vector<long double> next(static_cast<std::size_t>(dist.size() + sides), 0.0L);
            for (std::size_t sum = 0; sum < dist.size(); ++sum) {
                const long double base = dist[sum] / static_cast<long double>(sides);
                for (int face = 1; face <= sides; ++face) {
                    next[sum + static_cast<std::size_t>(face)] += base;
                }
            }
            dist.swap(next);
        }

        for (std::size_t sum = 0; sum < dist.size(); ++sum) {
            total_dist[sum] += p_count * dist[sum];
        }
    }

    long double mean = 0.0L;
    for (std::size_t sum = 0; sum < total_dist.size(); ++sum) {
        mean += total_dist[sum] * static_cast<long double>(sum);
    }

    long double variance = 0.0L;
    for (std::size_t sum = 0; sum < total_dist.size(); ++sum) {
        const long double diff = static_cast<long double>(sum) - mean;
        variance += total_dist[sum] * diff * diff;
    }

    return Stats{mean, variance};
}

bool nearly_equal(const long double a, const long double b, const long double rel_tol = 1e-14L) {
    const long double scale = std::max(std::fabsl(a), std::fabsl(b));
    if (scale == 0.0L) {
        return true;
    }
    return std::fabsl(a - b) <= rel_tol * scale;
}

bool run_checkpoints() {
    const Stats t = uniform_die_stats(4);
    if (!nearly_equal(t.mean, 2.5L) || !nearly_equal(t.variance, 1.25L)) {
        std::cerr << "Checkpoint failed for T\n";
        return false;
    }

    const Stats formula_c = transition_sum_of_random_dice(t, 6);
    const Stats brute_c = brute_one_layer_uniform_count(4, 6);
    if (!nearly_equal(formula_c.mean, brute_c.mean) ||
        !nearly_equal(formula_c.variance, brute_c.variance)) {
        std::cerr << "Checkpoint failed for C (formula vs brute distribution)\n";
        return false;
    }

    return true;
}

long double solve_variance() {
    const Stats t = uniform_die_stats(4);
    const Stats c = transition_sum_of_random_dice(t, 6);
    const Stats o = transition_sum_of_random_dice(c, 8);
    const Stats d = transition_sum_of_random_dice(o, 12);
    const Stats i = transition_sum_of_random_dice(d, 20);
    return i.variance;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << std::fixed << std::setprecision(4) << static_cast<double>(solve_variance()) << '\n';
    return 0;
}

Python

def uniform_die_stats(sides):
    s = float(sides)
    return ((s + 1.0) / 2.0, (s * s - 1.0) / 12.0)

def transition_sum(count_mean, count_var, sides):
    per_mean, per_var = uniform_die_stats(sides)
    mean = per_mean * count_mean
    var = per_var * count_mean + per_mean * per_mean * count_var
    return (mean, var)

def solve():
    mean, var = uniform_die_stats(4)
    for sides in [6, 8, 12, 20]:
        mean, var = transition_sum(mean, var, sides)
    return "{:.4f}".format(var)

if __name__ == '__main__':
    print(solve())

Java

public class Euler389 {
    static class Stats {
        double mean;
        double variance;

        Stats(double m, double v) {
            mean = m;
            variance = v;
        }
    }

    private static Stats uniformDieStats(int sides) {
        double s = sides;
        return new Stats((s + 1.0) / 2.0, (s * s - 1.0) / 12.0);
    }

    private static Stats transitionSum(Stats countStats, int sides) {
        Stats perDie = uniformDieStats(sides);
        double mean = perDie.mean * countStats.mean;
        double var = perDie.variance * countStats.mean + perDie.mean * perDie.mean * countStats.variance;
        return new Stats(mean, var);
    }

    public static String solve() {
        Stats s = uniformDieStats(4);
        s = transitionSum(s, 6);
        s = transitionSum(s, 8);
        s = transitionSum(s, 12);
        s = transitionSum(s, 20);
        return String.format("%.4f", s.variance).replace(',', '.');
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}