Dice game probabilities: Difference between revisions
Walterpachl (talk | contribs) m (→{{header|Rexx}}: corrected source to ooRexx) |
Walterpachl (talk | contribs) (→{{header|ooRexx}}: add algorithm) |
||
Line 155: | Line 155: | ||
=={{header|ooRexx}}== |
=={{header|ooRexx}}== |
||
===Algorithm=== |
|||
<lang oorexx>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+=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 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'+=1;' |
|||
Do wi=1 To w |
|||
str=str||'End;' |
|||
End |
|||
/* Say str */ |
|||
Interpret str |
|||
psum=0 |
|||
Do x=0 To w*s |
|||
p.x=cnt.x/(s**w) |
|||
psum+=p.x |
|||
/*Say format(x,2) format(cnt.x,3) format(p.x,8,3)*/ |
|||
End |
|||
/*Say psum*/ |
|||
res='' |
|||
Do x=w To s*w |
|||
res=res p.x |
|||
End |
|||
Return res</lang> |
|||
{{out}} |
|||
<pre>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.573144077 |
|||
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.642788632</pre> |
|||
===Test=== |
|||
Result from 10 million tries. |
Result from 10 million tries. |
||
<lang oorexx>oid='diet.xxx'; Call sysFileDelete oid |
<lang oorexx>oid='diet.xxx'; Call sysFileDelete oid |
Revision as of 23:25, 18 January 2015
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
<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.6427886287176260
Faster Alternative Version
<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.6427886287176260
ooRexx
Algorithm
<lang oorexx>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+=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 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'+=1;' Do wi=1 To w
str=str||'End;' End
/* Say str */ Interpret str psum=0 Do x=0 To w*s
p.x=cnt.x/(s**w) psum+=p.x /*Say format(x,2) format(cnt.x,3) format(p.x,8,3)*/ End
/*Say psum*/ res= Do x=w To s*w
res=res p.x End
Return res</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.573144077 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.642788632
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
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
Rexx
(actually just copied)
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: 9 dice with 4 sides each Player 2: 6 dice with 6 sides each 0.57406 player 1 wins 0.35606 player 2 wins 0.06988 draws 1.190000 seconds elapsed Player 1: 5 dice with 10 sides each Player 2: 6 dice with 7 sides each 0.64066 player 1 wins 0.31463 player 2 wins 0.04471 draws 0.952000 seconds elapsed