Dice game probabilities

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>

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(10, 5, 7, 6));
   return 0;

}</lang>

Output:

0.5731440767829801
0.6427886287176260

D

version 1

<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!(10, 5, 7, 6));

}</lang>

Output:

0.5731440767829801
0.6427886287176262

version 2 (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, 11u), 5, iota(1u, 8u), 6));

}</lang>

Output:

0.5731440767829801
0.6427886287176262

ooRexx

Algorithm

<lang oorexx>Numeric Digits 30 Call test '9 4 6 6' Call test '5 10 6 7' Exit test: Parse Arg w1 s1 w2 s2 p1.=0 p2.=0 Call pp 1,w1,s1,p1.,p2. Call pp 2,w2,s2,p1.,p2. p2low.=0 Do x=w1 To w1*s1

 Do y=0 To x-1
   p2low.x+=p2.y
   End
 End

pwin1=0 Do x=w1 To w1*s1

 pwin1+=p1.x*p2low.x
 End

Say 'Player 1 has' w1 'dice with' s1 'sides each' Say 'Player 2 has' w2 'dice with' s2 'sides each' Say 'Probability for player 1 to win:' pwin1 Say Return

pp: Procedure /*---------------------------------------------------------------------

Compute and assign probabilities to get a sum x
when throwing w dice each having s sides (marked from 1 to s)
k=1 sets p1.*, k=2 sets p2.*
--------------------------------------------------------------------*/

Use Arg k,w,s,p1.,p2. str= cnt.=0 Do wi=1 To w

 str=str||'Do v'wi'=1 To' s';'
 End

str=str||'sum=' Do wi=1 To w-1

 str=str||'v'wi'+'
 End

str=str||'v'w';' str=str||'cnt.'sum'+=1;' Do wi=1 To w

 str=str||'End;'
 End

Interpret str psum=0 Do x=0 To w*s

 If k=1 Then
   p1.x=cnt.x/(s**w)
 Else
   p2.x=cnt.x/(s**w)
 psum+=p1.x
 End

Return</lang>

Output:

Player 1 has 9 dice with 4 sides each
Player 2 has 6 dice with 6 sides each
Probability for player 1 to win: 0.573144076782980082947530864198

Player 1 has 5 dice with 10 sides each
Player 2 has 6 dice with 7 sides each
Probability for player 1 to win: 0.642788628717626159168373721835

Algorithm using fraction arithmetic

<lang oorexx>Numeric Digits 30 Call test '9 4 6 6' Call test '5 10 6 7' Exit

test: Parse Arg w1 s1 w2 s2 p1.='0/1' p2.='0/1' Call pp 1,w1,s1,p1.,p2. Call pp 2,w2,s2,p1.,p2. p2low.='0/1' Do x=w1 To w1*s1

 Do y=0 To x-1
   p2low.x=fr_add(p2low.x,p2.y)
   End
 End

pwin1='0/1' Do x=w1 To w1*s1

 pwin1=fr_add(pwin1,fr_Mult(p1.x,p2low.x))
 End

Say 'Player 1 has' w1 'dice with' s1 'sides each' Say 'Player 2 has' w2 'dice with' s2 'sides each' Say 'Probability for player 1 to win:' pwin1 Parse Var pwin1 nom '/' denom Say ' ->' (nom/denom) Say Return

pp: Procedure /*---------------------------------------------------------------------

Compute and assign probabilities to get a sum x
when throwing w dice each having s sides (marked from 1 to s)
k=1 sets p1.*, k=2 sets p2.*
--------------------------------------------------------------------*/

Use Arg k,w,s,p1.,p2. str= cnt.=0 Do wi=1 To w

 str=str||'Do v'wi'=1 To' s';'
 End

str=str||'sum=' Do wi=1 To w-1

 str=str||'v'wi'+'
 End

str=str||'v'w';' str=str||'cnt.sum+=1;' Do wi=1 To w

 str=str||'End;'
 End

Interpret str

psum='0/1' Do x=0 To w*s

 If k=1 Then
   p1.x=cnt.x'/'||(s**w)
 Else
   p2.x=cnt.x'/'||(s**w)
 psum=fr_add(psum,p1.x)
 End

Return

fr_add: Procedure Parse Arg a,b parse Var a an '/' az parse Var b bn '/' bz rn=an*bz+bn*az rz=az*bz res=fr_cancel(rn','rz) Return res

fr_div: Procedure Parse Arg a,b parse Var a an '/' az parse Var b bn '/' bz rn=an*bz rz=az*bn res=fr_cancel(rn','rz) Return res

fr_mult: Procedure Parse Arg a,b parse Var a an '/' az parse Var b bn '/' bz rn=an*bn rz=az*bz res=fr_cancel(rn','rz) Return res

fr_cancel: Procedure Parse Arg n ',' z k=ggt(n,z) Return n%k'/'z%k

ggt: Procedure /**********************************************************************

ggt (gcd) Greatest common Divisor
Recursive procedure as shown in PL/I
- - - - /

Parse Arg a,b if b = 0 then return abs(a) return ggt(b,a//b)</lang>

Output:

Player 1 has 9 dice with 4 sides each
Player 2 has 6 dice with 6 sides each
Probability for player 1 to win: 48679795/84934656
                              -> 0.573144076782980082947530864198

Player 1 has 5 dice with 10 sides each
Player 2 has 6 dice with 7 sides each
Probability for player 1 to win: 3781171969/5882450000
                              -> 0.642788628717626159168373721834

Test

Result from 10 million tries. <lang oorexx>oid='diet.xxx'; Call sysFileDelete oid Call test '9 4 6 6' Call test '5 10 6 7' Exit test: Parse Arg n1 s1 n2 s2 Call o 'Player 1:' n1 'dice with' s1 'sides each' Call o 'Player 2:' n2 'dice with' s2 'sides each' cnt1.=0 cnt2.=0 win.=0 nn=10000000 Call time 'R' Do i=1 To nn

 sum1=sum(n1 s1) ; cnt1.sum1+=1
 sum2=sum(n2 s2) ; cnt2.sum2+=1
 Select
   When sum1>sum2 Then win.1+=1
   When sum1<sum2 Then win.2+=1
   Otherwise           win.0+=1
   End
 End

Call o win.1/nn 'player 1 wins' Call o win.2/nn 'player 2 wins' Call o win.0/nn 'draws' /* Do i=min(n1,n2) To max(n1*s1,n2*s2)

 Call o right(i,2) format(cnt1.i,7) format(cnt2.i,7)
 End

/

Call o time('E') 'seconds elapsed' Return

sum: Parse Arg n s sum=0 Do k=1 To n

 sum+=rand(s)
 End

Return sum

rand: Parse Arg n

Return random(n-1)+1

o: Say arg(1) Return lineout(oid,arg(1))</lang>

Output:

Player 1: 9 dice with 4 sides each
Player 2: 6 dice with 6 sides each
0.5730344 player 1 wins
0.3562464 player 2 wins
0.0707192 draws
186.794000 seconds elapsed
Player 1: 5 dice with 10 sides each
Player 2: 6 dice with 7 sides each
0.6425906 player 1 wins
0.312976 player 2 wins
0.0444334 draws
149.784000 seconds elapsed

Perl 6

<lang perl6>sub likelihoods ($roll) {

   my ($dice, $faces) = $roll.comb(/\d+/);
   my @counts;
   @counts[$_]++ for [X+] |((1..$faces,) xx $dice);
   return [@counts[]:p], $faces ** $dice;

}

sub beating-probability ([$roll1, $roll2]) {

   my (@c1, $p1) := likelihoods $roll1;
   my (@c2, $p2) := likelihoods $roll2;
   my $p12 = $p1 * $p2;

   [+] gather for @c1 X @c2 -> $p, $q {

take $p.value * $q.value / $p12 if $p.key > $q.key;

}

We're using standard DnD notation for dice rolls here.

say .gist, "\t", .perl given beating-probability < 9d4 6d6 >; say .gist, "\t", .perl given beating-probability < 5d10 6d7 >;</lang>

Output:

0.573144077	<48679795/84934656>
0.64278862872	<3781171969/5882450000>

Note that all calculations are in integer and rational arithmetic, so the results in fractional notation are exact.

PL/I

version 1

<lang pli>*process source attributes xref;

dicegame: Proc Options(main);
Call test(9, 4,6,6);
Call test(5,10,6,7);
test: Proc(w1,s1,w2,s2);
Dcl (w1,s1,w2,s2,x,y) Bin Fixed(31);
Dcl p1(100)    Dec Float(18) Init((100)0);
Dcl p2(100)    Dec Float(18) Init((100)0);
Dcl p2low(100) Dec Float(18) Init((100)0);
Call pp(w1,s1,p1);
Call pp(w2,s2,p2);
Do x=w1 To w1*s1;
  Do y=0 To x-1;
    p2low(x)+=p2(y);
    End;
  End;
pwin1=0;
Do x=w1 To w1*s1;
 pwin1+=p1(x)*p2low(x);
 End;
Put Edit('Player 1 has ',w1,' dice with ',s1,' sides each')
        (Skip,3(a,f(2)));
Put Edit('Player 2 has ',w2,' dice with ',s2,' sides each')
        (Skip,3(a,f(2)));
Put Edit('Probability for player 1 to win: ',pwin1)(Skip,a,f(20,18));
Put Edit()(Skip,a);
End;
pp: Proc(w,s,p);
/*--------------------------------------------------------------------
* Compute and assign probabilities to get a sum x
* when throwing w dice each having s sides (marked from 1 to s)
*-------------------------------------------------------------------*/
Dcl (w,s)    Bin Fixed(31);
Dcl p(100)   Dec Float(18);
Dcl cnt(100) Bin Fixed(31);
Dcl (a(12),e(12),v(12),sum,i,n) Bin Fixed(31);
a=0;
e=0;
Do i=1 To w;
  a(i)=1;
  e(i)=s;
  End;
n=0;
cnt=0;
Do v(1)=a(1) To e(1);
  Do v(2)=a(2) To e(2);
    Do v(3)=a(3) To e(3);
      Do v(4)=a(4) To e(4);
        Do v(5)=a(5) To e(5);
          Do v(6)=a(6) To e(6);
            Do v(7)=a(7) To e(7);
              Do v(8)=a(8) To e(8);
                Do v(9)=a(9) To e(9);
                  Do v(10)=a(10) To e(10);
                    sum=0;
                    Do i=1 To 10;
                      sum=sum+v(i);
                      End;
                    cnt(sum)+=1;
                    n+=1;
                    End;
                  End;
                End;
              End;
            End;
          End;
        End;
      End;
    End;
  End;
Do k=w To w*s;
  p(k)=divide(cnt(k),n,18,16);
  End;
End;
End;</lang>

Output:

Player 1 has  9 dice with  4 sides each
Player 2 has  6 dice with  6 sides each
Probability for player 1 to win: 0.573013663291931152

Player 1 has  5 dice with 10 sides each
Player 2 has  6 dice with  7 sides each
Probability for player 1 to win: 0.642703175544738770

version 2 using fraction arithmetic

<lang pli>*process source attributes xref;

dgf: Proc Options(main);
Call test(9, 4,6,6);
Call test(5,10,6,7);
test: Proc(w1,s1,w2,s2);
Dcl (w1,s1,w2,s2,x,y) Dec Float(18);
Dcl 1 p1(100),
     2 nom      Dec Float(18) Init((100)0),
     2 denom    Dec Float(18) Init((100)1);
Dcl 1 p2(100),
     2 nom      Dec Float(18) Init((100)0),
     2 denom    Dec Float(18) Init((100)1);
Dcl 1 p2low(100),
     2 nom      Dec Float(18) Init((100)0),
     2 denom    Dec Float(18) Init((100)1);
Dcl 1 pwin1,
     2 nom      Dec Float(18) Init(0),
     2 denom    Dec Float(18) Init(1);
Dcl 1 prod Like pwin1;
Call pp(w1,s1,p1);
Call pp(w2,s2,p2);
Do x=w1 To w1*s1;
  Do y=0 To x-1;
    Call fr_add(p2low(x),p2(y),p2low(x));
    End;
  End;
Do x=w1 To w1*s1;
 Call fr_mult(p1(x),p2low(x),prod);
 Call fr_add(pwin1,prod,pwin1);
 End;
Put Edit('Player 1 has ',w1,' dice with ',s1,' sides each')
        (Skip,3(a,f(2)));
Put Edit('Player 2 has ',w2,' dice with ',s2,' sides each')
        (Skip,3(a,f(2)));
Put Edit('Probability for player 1 to win: ',
         str(pwin1.nom),'/',str(pwin1.denom))(Skip,4(a));
Put Edit('                              -> ',
         pwin1.nom/pwin1.denom)(Skip,a,f(20,18));
Put Edit()(Skip,a);
End;

pp: Proc(w,s,p);
/*--------------------------------------------------------------------
* Compute and assign probabilities to get a sum x
* when throwing w dice each having s sides (marked from 1 to s)
*-------------------------------------------------------------------*/
Dcl (w,s)    Dec Float(18);
Dcl 1 p(100),
     2 nom   Dec Float(18),
     2 denom Dec Float(18);
Dcl cnt(100) Dec Float(18);
Dcl (a(12),e(12),v(12),sum,i,n) Dec Float(18);
a=0;
e=0;
Do i=1 To w;
  a(i)=1;
  e(i)=s;
  End;
n=0;
cnt=0;
Do v(1)=a(1) To e(1);
  Do v(2)=a(2) To e(2);
    Do v(3)=a(3) To e(3);
      Do v(4)=a(4) To e(4);
        Do v(5)=a(5) To e(5);
          Do v(6)=a(6) To e(6);
            Do v(7)=a(7) To e(7);
              Do v(8)=a(8) To e(8);
                Do v(9)=a(9) To e(9);
                  Do v(10)=a(10) To e(10);
                    sum=0;
                    Do i=1 To 10;
                      sum=sum+v(i);
                      End;
                    cnt(sum)+=1;
                    n+=1;
                    End;
                  End;
                End;
              End;
            End;
          End;
        End;
      End;
    End;
  End;
Do k=w To w*s;
  p(k).nom=cnt(k);
  p(k).denom=n;
  End;
End;

fr_add: Proc(a,b,res);
Dcl 1 a,
     2 nom   Dec Float(18),
     2 denom Dec Float(18);
Dcl 1 b Like a;
Dcl res like a;
/* Put Edit('fr_add',a.nom,a.denom,b.nom,b.denom)(Skip,a,4(f(15))); */
res.nom=a.nom*b.denom+b.nom*a.denom;
res.denom=a.denom*b.denom;
Call fr_cancel(res,res);
End;

fr_mult: Proc(a,b,res);
Dcl 1 a,
     2 nom   Dec Float(18),
     2 denom Dec Float(18);
Dcl 1 b Like a;
Dcl res like a;
/* Put Edit('fr_mult',a.nom,a.denom,b.nom,b.denom)(Skip,a,4(f(15)));*/
res.nom=a.nom*b.nom;
res.denom=a.denom*b.denom;
Call fr_cancel(res,res);
End;

fr_cancel: Proc(a,res);
Dcl 1 a,
     2 nom   Dec Float(18),
     2 denom Dec Float(18);
Dcl k Dec Float(18);
Dcl 1 res like a;
/* Put Edit('fr_cancel',a.nom,a.denom)(Skip,a,4(f(15)));            */
k=ggt(a.nom,a.denom);
res=a/k;
End;

ggt: Proc(a,b) Recursive Returns(Dec Float(18));
/**********************************************************************
* ggt (gcd) Greatest common Divisor
* Recursive Proc(a,b)) as shown in PL/I
**********************************************************************/
Dcl (a,b) Dec Float(18);
if b = 0 then return (abs(a));
return (ggt(b,mod(a,b)));
End;

str: Proc(x) Returns(Char(20) Var);
Dcl x Dec Float(18);
Dcl res Char(20) Var;
Put String(res) Edit(x)(f(20));
Return (trim(res));
End;

End;</lang>

Output:

Player 1 has  9 dice with  4 sides each
Player 2 has  6 dice with  6 sides each
Probability for player 1 to win: 48679795/84934656
                              -> 0.573144076782980083

Player 1 has  5 dice with 10 sides each
Player 2 has  6 dice with  7 sides each
Probability for player 1 to win: 3781171969/5882450000
                              -> 0.642788628717626159

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(10, 5, 7, 6)</lang>

Output:

0.573144076783
0.642788628718

To handle larger number of dice (and faster in general): <lang python>from __future__ import print_function, 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,11), 5, range(1,8), 6)) # this seem hardly fair

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

Output:

(0.5731440767829801, 0.070766169838454, 0.3560897533785659)
(0.6427886287176262, 0.044496030310499875, 0.312715340971874)

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))

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

Output:

0.5731440767829801
0.6427886287176262

Racket

<lang racket>#lang racket

(define probs# (make-hash))

(define (NdD n d)

 (hash-ref!
  probs# (cons n d)
  (λ ()
    (cond
      [(= n 0) ; every chance of nothing!
       (hash 0 1)]
      [else
       (for*/fold ((hsh (hash))) (((i p) (in-hash (NdD (sub1 n) d))) (r (in-range 1 (+ d 1))))
         (hash-update hsh (+ r i) (λ (p+) (+ p+ (/ p d))) 0))]))))

(define (game-probs N1 D1 N2 D2)

 (define P1 (NdD N1 D1))
 (define P2 (NdD N2 D2))  
 (define-values (W D L)
   (for*/fold ((win 0) (draw 0) (lose 0)) (((r1 p1) (in-hash P1)) ((r2 p2) (in-hash P2)))
     (define p (* p1 p2))
     (cond
       [(< r1 r2) (values win draw (+ lose p))]
       [(= r1 r2) (values win (+ draw p) lose)]
       [(> r1 r2) (values (+ win p) draw lose)])))
 
 (printf "P(P1 win): ~a~%" (real->decimal-string W 6))
 (printf "P(draw):   ~a~%" (real->decimal-string D 6))
 (printf "P(P2 win): ~a~%" (real->decimal-string L 6))  
 (list W D L))

(printf "GAME 1 (9D4 vs 6D6)~%") (game-probs 9 4 6 6) (newline)

(printf "GAME 2 (5D10 vs 6D7) [what is a D7?]~%") (game-probs 5 10 6 7)</lang>

Output:

GAME 1 (9D4 vs 6D6)
P(P1 win): 0.573144
P(draw):   0.070766
P(P2 win): 0.356090
(48679795/84934656 144252007/2038431744 725864657/2038431744)

GAME 2 (5D10 vs 6D7) [what is a D7?]
P(P1 win): 0.642789
P(draw):   0.044496
P(P2 win): 0.312715
(3781171969/5882450000 523491347/11764900000 735812943/2352980000)

REXX

Algorithm

Translation of: ooRexx

(adapted for Classic Rexx) <lang rexx>Numeric Digits 30 Call test '9 4 6 6' Call test '5 10 6 7' Exit test: Parse Arg w1 s1 w2 s2 plist1=pp(w1,s1) p1.=0 Do x=w1 To w1*s1

 Parse Var plist1 p1.x plist1
 End

plist2=pp(w2,s2) p2.=0 Do x=w2 To w2*s2

 Parse Var plist2 p2.x plist2
 End

p2low.=0 Do x=w1 To w1*s1

 Do y=0 To x-1
   p2low.x=p2low.x+p2.y
   End
 End

pwin1=0 Do x=w1 To w1*s1

 pwin1=pwin1+p1.x*p2low.x
 End

Say 'Player 1 has' w1 'dice with' s1 'sides each' Say 'Player 2 has' w2 'dice with' s2 'sides each' Say 'Probability for player 1 to win:' pwin1 Say Return

pp: Procedure /*---------------------------------------------------------------------

Compute and return the probabilities to get a sum x
when throwing w dice each having s sides (marked from 1 to s)
--------------------------------------------------------------------*/

Parse Arg w,s str= cnt.=0 Do wi=1 To w

 str=str||'Do v'wi'=1 To' s';'
 End

str=str||'sum=' Do wi=1 To w-1

 str=str||'v'wi'+'
 End

str=str||'v'w';' str=str||'cnt.'sum'=cnt.'sum'+1;' Do wi=1 To w

 str=str||'End;'
 End

Interpret str psum=0 Do x=0 To w*s

 p.x=cnt.x/(s**w)
 psum=psum+p.x
 End

res= Do x=w To s*w

 res=res p.x
 End

Return res
Result from 100.000 tries. <lang rexx>oid='diet.xxx'; Call sysFileDelete oid Call test '9 4 6 6' Call test '5 10 6 7' Exit test: Parse Arg n1 s1 n2 s2 Call o 'Player 1:' n1 'dice with' s1 'sides each' Call o 'Player 2:' n2 'dice with' s2 'sides each' cnt1.=0 cnt2.=0 win.=0 nn=100000 Call time 'R' Do i=1 To nn

 sum1=sum(n1 s1) ; cnt1.sum1+=1
 sum2=sum(n2 s2) ; cnt2.sum2+=1
 Select
   When sum1>sum2 Then win.1+=1
   When sum1<sum2 Then win.2+=1
   Otherwise           win.0+=1
   End
 End

Call o win.1/nn 'player 1 wins' Call o win.2/nn 'player 2 wins' Call o win.0/nn 'draws' /* Do i=min(n1,n2) To max(n1*s1,n2*s2)

 Call o right(i,2) format(cnt1.i,7) format(cnt2.i,7)
 End

/

Call o time('E') 'seconds elapsed' Return

sum: Parse Arg n s sum=0 Do k=1 To n

 sum+=rand(s)
 End

Return sum

rand: Parse Arg n

Return random(n-1)+1

o: Say arg(1) Return lineout(oid,arg(1))</lang>

Output:

Player 1 has 9 dice with 4 sides each
Player 2 has 6 dice with 6 sides each
Probability for player 1 to win: 0.573144076782980082947530864198

Player 1 has 5 dice with 10 sides each
Player 2 has 6 dice with 7 sides each
Probability for player 1 to win: 0.642788628717626159168373721835