Dice game probabilities

Dice game probabilities is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Two players have a set of dice each. The first player has nine dice with four faces each, with numbers one to four. The second player has six normal dice with six faces each, each face has the usual numbers from one to six.

They roll their dice and sum the totals of the faces. The player with the highest total wins (it's a draw if the totals are the same). What's the probability of the first player beating the second player?

Later the two players use a different set of dice each. Now the first player has five dice with ten faces each, and the second player has six dice with seven faces each. Now what's the probability of the first player beating the second player?

This task was adapted from the Project Euler Problem n.205: https://projecteuler.net/problem=205

C

<lang c>#include <stdio.h>

1. include <stdint.h>

typedef uint32_t uint; typedef uint64_t ulong;

ulong ipow(const uint x, const uint y) {

   ulong result = 1;
   for (uint i = 1; i <= y; i++)
       result *= x;
   return result;

}

uint min(const uint x, const uint y) {

   return (x < y) ? x : y;

}

void throw_die(const uint n_sides, const uint n_dice, const uint s, uint counts[]) {

   if (n_dice == 0) {
       counts[s]++;
       return;
   }

   for (uint i = 1; i < n_sides + 1; i++)
       throw_die(n_sides, n_dice - 1, s + i, counts);

}

double beating_probability(const uint n_sides1, const uint n_dice1,

                          const uint n_sides2, const uint n_dice2) {
   const uint len1 = (n_sides1 + 1) * n_dice1;
   uint C1[len1];
   for (uint i = 0; i < len1; i++)
       C1[i] = 0;
   throw_die(n_sides1, n_dice1, 0, C1);

   const uint len2 = (n_sides2 + 1) * n_dice2;
   uint C2[len2];
   for (uint j = 0; j < len2; j++)
       C2[j] = 0;
   throw_die(n_sides2, n_dice2, 0, C2);

   const double p12 = (double)(ipow(n_sides1, n_dice1) * ipow(n_sides2, n_dice2));

   double tot = 0;
   for (uint i = 0; i < len1; i++)
       for (uint j = 0; j < min(i, len2); j++)
           tot += (double)C1[i] * C2[j] / p12;
   return tot;

}

int main() {

   printf("%1.16f\n", beating_probability(4, 9, 6, 6));
   printf("%1.16f\n", beating_probability(5, 10, 6, 7));
   return 0;

}</lang>

Output:

0.5731440767829801
0.7826269636740740

D

<lang d>import std.stdio, std.range, std.algorithm;

void throwDie(in uint nSides, in uint nDice, in uint s, uint[] counts) pure nothrow @safe @nogc {

   if (nDice == 0) {
       counts[s]++;
       return;
   }

   foreach (immutable i; 1 .. nSides + 1)
       throwDie(nSides, nDice - 1, s + i, counts);

}

real beatingProbability(uint nSides1, uint nDice1,

                       uint nSides2, uint nDice2)()

pure nothrow @safe /*@nogc*/ {

   uint[(nSides1 + 1) * nDice1] C1;
   throwDie(nSides1, nDice1, 0, C1);

   uint[(nSides2 + 1) * nDice2] C2;
   throwDie(nSides2, nDice2, 0, C2);

   immutable p12 = real((ulong(nSides1) ^^ nDice1) *
                        (ulong(nSides2) ^^ nDice2));

   return cartesianProduct(C1[].enumerate, C2[].enumerate)
          .filter!(p => p[0][0] > p[1][0])
          .map!(p => real(p[0][1]) * p[1][1] / p12)
          .sum;

}

void main() @safe {

   writefln("%1.16f", beatingProbability!(4, 9, 6, 6));
   writefln("%1.16f", beatingProbability!(5, 10, 6, 7));

}</lang>

Output:

0.5731440767829801
0.7826269636740741

Faster Alternative Version

Translation of: Python

<lang d>import std.stdio, std.range, std.algorithm;

ulong[] combos(R)(R sides, in uint n) pure nothrow @safe if (isForwardRange!R) {

   if (sides.empty)
       return null;
   if (!n)
       return [1];
   auto ret = new typeof(return)(reduce!max(sides[0], sides[1 .. $]) * n + 1);
   foreach (immutable i, immutable v; enumerate(combos(sides, n - 1))) {
       if (!v)
           continue;
       foreach (immutable s; sides)
           ret[i + s] += v;
   }
   return ret;

}

real winning(R)(R sides1, in uint n1, R sides2, in uint n2) pure nothrow @safe if (isForwardRange!R) {

   static void accumulate(T)(T[] arr) pure nothrow @safe @nogc {
       foreach (immutable i; 1 .. arr.length)
           arr[i] += arr[i - 1];
   }

   immutable p1 = combos(sides1, n1);
   auto p2 = combos(sides2, n2);
   immutable s = p1.sum * p2.sum;
   accumulate(p2);
   ulong win = 0; // 'win' is 1 beating 2.
   foreach (immutable i, immutable x1; p1.dropOne.enumerate)
       win += x1 * p2[min(i, $ - 1)];
   return win / real(s);

}

void main() @safe {

   writefln("%1.16f", winning(iota(1u, 5u),  9, iota(1u, 7u), 6));
   writefln("%1.16f", winning(iota(1u, 6u), 10, iota(1u, 7u), 7));

}</lang>

Output:

0.5731440767829801
0.7826269636740741

Python

<lang python>from itertools import product

def gen_dict(n_faces, n_dice):

   counts = [0] * ((n_faces + 1) * n_dice)
   for t in product(range(1, n_faces + 1), repeat=n_dice):
       counts[sum(t)] += 1
   return counts, n_faces ** n_dice

def beating_probability(n_sides1, n_dice1, n_sides2, n_dice2):

   c1, p1 = gen_dict(n_sides1, n_dice1)
   c2, p2 = gen_dict(n_sides2, n_dice2)
   p12 = float(p1 * p2)

   return sum(p[1] * q[1] / p12
              for p, q in product(enumerate(c1), enumerate(c2))
              if p[0] > q[0])

print beating_probability(4, 9, 6, 6) print beating_probability(5, 10, 6, 7)</lang>

Output:

0.573144076783
0.782626963674

To handle larger number of dice (and faster in general): <lang python>from __future__ import print_function from __future__ import division

def combos(sides, n):

   if not n: return [1]
   ret = [0] * (max(sides)*n + 1)
   for i,v in enumerate(combos(sides, n - 1)):
       if not v: continue
       for s in sides: ret[i + s] += v
   return ret

def winning(sides1, n1, sides2, n2):

   p1, p2 = combos(sides1, n1), combos(sides2, n2)
   win,loss,tie = 0,0,0 # 'win' is 1 beating 2
   for i,x1 in enumerate(p1):
       # using accumulated sum on p2 could save some time
       win += x1*sum(p2[:i])
       tie += x1*sum(p2[i:i+1])
       loss+= x1*sum(p2[i+1:])
   s = sum(p1)*sum(p2)
   return win/s, tie/s, loss/s

print(winning(range(1,5), 9, range(1,7), 6)) print(winning(range(1,6), 10, range(1,7), 7)) # this seem hardly fair

1. mountains of dice test case
2. print(winning((1, 2, 3, 5, 9), 700, (1, 2, 3, 4, 5, 6), 800))</lang>

Output:

(0.5731440767829801, 0.070766169838454, 0.3560897533785659)
(0.782626963674074, 0.043287775394238684, 0.17408526093168725)

If we further restrict die faces to be 1 to n instead of arbitrary values, the combo generation can be made much faster: <lang python>from __future__ import division, print_function from itertools import accumulate # Python3 only

def combos(sides, n):

   ret = [1] + [0]*(n + 1)*sides # extra length for negative indices
   for p in range(1, n + 1):
       rolling_sum = 0
       for i in range(p*sides, p - 1, -1):
           rolling_sum += ret[i - sides] - ret[i]
           ret[i] = rolling_sum
       ret[p - 1] = 0
   return ret

def winning(d1, n1, d2, n2):

   c1, c2 = combos(d1, n1), combos(d2, n2)
   ac = list(accumulate(c2 + [0]*(len(c1) - len(c2))))

   return sum(v*a for  v,a in zip(c1[1:], ac)) / (ac[-1]*sum(c1))

print(winning(4, 9, 6, 6)) print(winning(5, 10, 6, 7))

1. print(winning(6, 700, 8, 540))</lang>

Output:

0.5731440767829801
0.782626963674074