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 , we need both dice to have the identity (otherwise, we will not be able to reach ). So, .

Now, consider the prime . It must be on both dice, or we would have . So, . Also, since appears on both dice, no dice can contain some product of the primes and their powers (e.g ) that does not exist on the original dice, because then impossible products could be reached.

Hence, must be on both dice, giving . There are sides left on the larger die but we have more even products, so must also be on each die. . Now, there is no space left for on the smaller die. This means that must be one the larger die, but then , which is a contradiction.

(@FlashDiaz gave a shorter proof)

Project Euler 205

Peter has nine four-sided (pyramidal) dice, each with faces numbered .
Colin has six six-sided (cubic) dice, each with faces numbered .

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 .

The probability functions of the nine four-sided dice and the six six-sided dice are given by the generating functions and , respectively.

Let be i.i.d random variables taking values in the range and let taking values in the range . We want to determine the probability that .

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
        else:
            for i in range(1, sides+1):
                rec_compute_dist(sides, nbr-1, side_sum+i)

dist = [0]*37
rec_compute_dist(4,9,0)
dist_49 = dist
dist = [0]*37
rec_compute_dist(6,6,0)
dist_66 = dist

To determine , we may express it as

.

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 ways in which five six-sided dice (sides numbered to ) can be rolled so that the top three sum to . Some examples are:

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

Let us first consider the simpler problem . 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 denote the number of ‘s we got. Then,  where .

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:
            possible.append(copy(d))
        return
    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 . 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 , where is the left-most non-zero index (starting from ). The number of configurations for this particular solution is then given by

, where .

Unfortunately, there is no good analytical way of computing this. So, the easiest way is to enumerate all possible . Disregarding , 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 \
                       Counter(x).values()])

        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
        else:
            q = fact(DICE) / u
        H[min_index] += 1
        S += q
print S