Problem of dice

The initial purpose of this post was give a proper proof of a problem posted on Twitter by @jamestanton (it is hard within the 140 char limit), but the post was later extended to cover some other problems.

The Twitter-problem
Show that a four-sided and a nine-sided dice cannot be used to simulate the probability distribution of the product of outcomes when using two six-sided dice. In the original wording:

Is there a 4-sided die & a 9-sided die that together roll each of the products 1,2,3,…,30,36 w the same prob as two ordinary 6-sided dice?

We make an argument by contradiction, considering only the possible outcomes without taking the actual probabilities into account.

Obviously, to reach the same outcomes as for two normal dice \{1,2,3,4,5,6\} \times \{1,2,3,4,5,6\}, we need both dice to have the identity \{1\} (otherwise, we will not be able to reach 1 \cdot 1 = 1). So, \{1,*,*,*\} \times \{1,*,*,*,*,*,*,*,*\}.

Now, consider the prime 5. It must be on both dice, or we would have \mathbb{P}(5^2\cdot b)>0, b>1. So, \{1,5,*,*\} \times \{1,5,*,*,*,*,*,*,*\}. Also, since 5 appears on both dice, no dice can contain some product of the primes \{2,3,5\} and their powers (e.g 2^2 \cdot 3) that does not exist on the original dice, because then impossible products could be reached.

Hence, 6 must be on both dice, giving \{1,5,6,*\} \times \{1,5,6,*,*,*,*,*,*\}. There are 6 sides left on the larger die but we have more even products, so 2 must also be on each die. \{1,5,6,2\} \times \{1,5,6,2,*,*,*,*,*\}. Now, there is no space left for 3 on the smaller die. This means that 3^2 must be one the larger die, but then \mathbb{P}(3^2\cdot 5)>0, which is a contradiction.

(@FlashDiaz gave a shorter proof)

Project Euler 205

Peter has nine four-sided (pyramidal) dice, each with faces numbered 1, 2, 3, 4.
Colin has six six-sided (cubic) dice, each with faces numbered 1, 2, 3, 4, 5, 6.

Peter and Colin roll their dice and compare totals: the highest total wins. The result is a draw if the totals are equal.

What is the probability that Pyramidal Pete beats Cubic Colin? Give your answer rounded to seven decimal places in the form 0.abcdefg.

The probability functions of the nine four-sided dice and the six six-sided dice are given by the generating functions \frac{1}{4^9} \cdot (x^1+x^2+x^3+x^4)^9 and \frac{1}{6^6} \cdot (y^1+y^2+y^3+y^4+y^5+y^6)^6, respectively.

Let X_1,...,X_9 be i.i.d random variables taking values in the range [1,4] and let Y_1,...,Y_6 taking values in the range [1,6]. We want to determine the probability that \rho = \mathbb{P}(X_1+...+X_9 > Y_1+...+Y_6).

The distributions can be computed as

def rec_compute_dist(sides, nbr, side_sum):
        global dist
        if nbr == 1:
            for i in range(1, sides+1):
                dist[side_sum+i] += 1
            for i in range(1, sides+1):
                rec_compute_dist(sides, nbr-1, side_sum+i)

dist = [0]*37
dist_49 = dist
dist = [0]*37
dist_66 = dist

To determine \rho, we may express it as

\begin{array}{rl} \rho = & \sum_{t=6}^{36} \mathbb{P}(X_1+...+X_9 > t| Y_1+...+Y_6 = t)\cdot \mathbb{P}(Y_1+...+Y_6 = t) \\\\ = & \sum_{t=6}^{36} \mathbb{P}(X_1+...+X_9 > t)\cdot \mathbb{P}(Y_1+...+Y_6 = t) \end{array} .

Computing the sum using the following code,

probability = 0
for i in range(6,36+1):
    for j in range(i+1,36+1):
        probability += dist_66[i]*dist_49[j]

print 1.0 * probability/(6**6 * 4**9)

we obtain the answer. Great 🙂

Project Euler 240

There are 1111 ways in which five six-sided dice (sides numbered 1 to 6) can be rolled so that the top three sum to 15. Some examples are:

\begin{array}{rcl} D_1,D_2,D_3,D_4,D_5 &=& 4,3,6,3,5\\ D_1,D_2,D_3,D_4,D_5 &=& 4,3,3,5,6\\ D_1,D_2,D_3,D_4,D_5 &=& 3,3,3,6,6\\ D_1,D_2,D_3,D_4,D_5 &=& 6,6,3,3,3 \end{array}

In how many ways can twenty twelve-sided dice (sides numbered 1 to 12) be rolled so that the top ten sum to 70?

Let us first consider the simpler problem \left\{ d_1+..+d_{10}=70 \right\}. If we restrict the remaining ten dice to be less than or equal to the minimum value of the ten dice, we then can compute the cardinality. Let n_i denote the number of i‘s we got. Then, n_1 \cdot 1 + n_2 \cdot 2 + ... + n_{12} \cdot 12 = 70 where n_1  + n_2 + ... + n_{12}  = 10, n_i \geq 0.

All histograms of top-ten dice can be computed with

from copy import copy
d = [0] * 12
possible = []

def rec_compute(i, j, sum):
    global d
    if j == 0:
        if sum == 70:
    while i > 0:
        if sum + i <= 70:
            d[i - 1] += 1
            rec_compute(i, j - 1, sum + i)
            d[i - 1] -= 1
        i -= 1

rec_compute(12, 10, 0)

The code exhausts all solutions in 200ms. Call any solution H. For instance H = [0, 0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0]. The remaining dice can take any values in the range [1, j], where j is the left-most non-zero index (starting from 1). The number of configurations for this particular solution is then given by

20! \cdot \left((10+H_7)!H_6!H_5!H_4!H_3!H_2!H_1!\right)^{-1}, where \sum^7_1 H_i = 10.

Unfortunately, there is no good analytical way of computing this. So, the easiest way is to enumerate all possible H_i. Disregarding H_7, we compute all permutations of a given histogram in the same way (hence, we can make the search space a lot smaller) and the using the multiplicity to determine the exact number. All and all, the following code gives our answer:

def configurations(i, j, x, s, l):
    if sum(x) == s:
        # we can stop here as the sum cannot get smaller
        multiplicity = fact(l) / fact(l-len(x)) / \
                       reduce(lambda m, n: m * n, \
                       [fact(y) for y in \

        return fact(DICE) * multiplicity / \
               reduce(lambda m, n: m * n, \
               [fact(y) for y in x])

    if j == 0 or i == 0: return 0
    return configurations(i-1, j, x, s, l) + \
           configurations(i, j-1, x + [i], s, l)

S = 0
for H in possible_top_dice:
    min_index = next((i for i, \
                x in enumerate(H) if x), None)

    for j in range(0, REMAINING+1):
        u = reduce(lambda m, n: m * n, \
            [fact(y) for y in H])

        # mutually exclusive -- may add instead of union
        if j < REMAINING:
            q = configurations(REMAINING-j, min_index, \
                [], REMAINING-j, min_index) / u
            q = fact(DICE) / u
        H[min_index] += 1
        S += q
print S

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s