Euler Problem 205 asks:
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
As a quick review, if we roll 2 of Colin’s dice, we expect 62 different outcomes. Rolling 6 dice will have 66 outcomes, or 46,656 different rolls.
Peter has 49 different outcomes, or 262,144 different rolls. Peter’s least roll (nine 1’s) will best the one way Colin can roll a 6 (six 1’s), the six ways he can roll a 7 (five 1’s and a 2 six times) and the twenty-one ways he can roll an 8 (a 3 and five 1’s six times, or two 2’s and four 1’s fifteen times). Peter’s meager 9 wins over 28 of Colin’s possible rolls. Peter’s 10, which he can roll 9 ways, bests 84 of Colin’s rolls.
VBA does not have a CEILING function, and I needed one for this problem. We could use Application.Worksheetfunction.Ceiling, but there is a quicker way execution-wise by a factor of 5. The INT function always rounds down. When the argument to INT is negative, INT rounds down or away from zero. INT(-3.14159) is -4, and -INT(-3.14159) is 4, rounding pi() up! Very useful when you need more area in your circles. It works this way in both the VBA and the spreadsheet implementations.
Easier than developing the usage for Problem 205, I’ll show it and explain how it works. The code we want to use for Colin is “-INT(-N/6^C) Mod 6” for C from zero to five, where N is the number of the roll (1 to 66), and when Mod 6 = zero, substitute 6. In a spreadsheet, this would be =IF(MOD(-INT(-N/6^C),6), MOD(-INT(-N/6^C),6), 6)
Remembering that 60 is 1, and 61 is 6, this is how the first four of Colin’s dice (C=0,1,2,3) look on Roll 66, N = 66, -N = -66.
- C = 0, Die 1:
- INT(-66/6^0) = INT(-66/1) = -66
- 66 = 66
- 66 Mod 6 = 0, Return 6
- C = 1, Die 2:
- INT(-66/6^1) = INT(-66/6) = -11
- 11 = 11
- 11 Mod 6 = 5, Return 5
- C = 2, Die 3:
- INT(-66/6^2) = INT(-66/36) = INT(-1.83333) = -2
- 2 = 2
- 2 Mod 6 = 2, Return 2
- C = 3, Die 4:
- INT(-66/6^3) = INT(-66/216) = INT(-0.30555) = -1
- 1 = 1
- 1 Mod 6 = 1, Return 1
Dice 5 and 6 (with C of 4 and 5) also return 1. Colin’s 66th roll is {6,5,2,1,1,1}. We do the same thing for Peter, where the code is “-INT(-N/4^P) Mod 4” for P from zero to eight, returning 4 when Mod 4 is zero. Peter’s 66th roll is {2,1,1,2,1,1,1,1,1}, summing 11. Peter gets 11 forty-five ways, on which he beats the 210 of Colin’s rolls (but not Colin’s #66) that sum 10 or below.
If we aggregate the number of times Colin sums from 6 to 36 in his 46,656 possible rolls, and the number of times Peter gets a particular sum from 9 to 36 in his 262,144 different rolls, we can then loop through Peter’s aggregation and see how many of Colin’s rolls lose to that number. If we then multiply that discovery by the number of ways Peter achives that aggregation, keep a grand sum of winners, and then divide by the product of (66)*(49),we will have our percentage of Peter’s winning. Format the answer to 7 decimals to the right. Format() will take care of the necessary rounding.
This is the code that does this. It runs in about 6ths of a second.
Dim N As Long, TEMP As Long, Sum As Long
Dim Answer As Double, T As Single, Count As Double
Dim PP(1 To 36) As Long, P As Long ‘Pyramidal Pete
Dim CC(1 To 36) As Long, C As Long ‘Cubic Colin
Dim LosersToPete As Long
T = Timer
For N = 1 To 6 ^ 6 ‘Cubic Colin
Sum = 0
For C = 0 To 5
TEMP = -Int(-N / 6 ^ C) Mod 6
If TEMP = 0 Then TEMP = 6
Sum = Sum + TEMP
Next C
CC(Sum) = CC(Sum) + 1 ‘Incrementing Colin’s ways this value can happen
Next N
For N = 1 To 4 ^ 9 ‘Pyramidal Pete
Sum = 0
For P = 0 To 8
TEMP = -Int(-N / 4 ^ P) Mod 4
If TEMP = 0 Then TEMP = 4
Sum = Sum + TEMP
Next P
PP(Sum) = PP(Sum) + 1 ‘Incrementing Pete’s ways this value can happen
Next N
For P = 9 To 36 ‘ Pete’s rolls
LosersToPete = 0
For C = 6 To P – 1
LosersToPete = LosersToPete + CC(C) ‘ Num Colin’s rolls (all losses) below Pete’s roll
Next C
Count = Count + (LosersToPete * PP(P))
‘Incrementing the winning Count by the # of ways Colin’s roll can lose to Pete
Next P
Answer = Count / (CDbl(4 ^ 9) * CDbl(6 ^ 6))
Debug.Print Format(Answer, “0.0000000”);” Time:”; Timer – T; Count
End Sub
If, instead of -INT, I use Application.Worksheetfunction.Ceiling as:
- TEMP = Application.WorksheetFunction.Ceiling(N / 6 ^ C, 1) Mod 6 and
- TEMP = Application.WorksheetFunction.Ceiling(N / 4 ^ P, 1) Mod 4
the runtime is 3.5 seconds! Using ROUNDUP() is even slower. The really wrong way to do this problem is to match each of Peter’s rolls with each of Colin’s, or something like this, where larger PP() and CC() now hold each roll and not the occurrances of each sum.
For C = 1 to 6^6
If PP(P) > CC(C) then Count = Count + 1
Next C
Next P
That’s 12,230,590,464 loops. Been there, did that. Takes 6 and a half minutes. No tee shirt.
…mrt
I have always found worksheet functions run appallingly slowly in VBA (unless you need to do something like sort a worksheet of course). MAX and MIN run about 50x slower than the equivalent VBA code.
Recursion is the ideal way to calculate the probability distributions for the different sets of dice.
‘n = number of dice, f = number of faces (assume points go from 1 to f)
‘optional s parameter returns a single probability value when given
‘otherwise returns the probability distribution as a 2D array
Dim j As Long, k As Long, lo As Long, hi As Long, u As Long
Dim rv() As Double, prv As Variant
If n LT 1 Or f LT 1 Then
dice = CVErr(xlErrNum)
ElseIf IsMissing(s) Then
‘return full probability table as an array
ElseIf Not IsNumeric(s) Then
dice = CVErr(xlErrValue)
ElseIf CDbl(s) NE CLng(s) Then
dice = CVErr(xlErrValue)
ElseIf CLng(s) LT n Or n * f LT CLng(s) Then
dice = CVErr(xlErrNA)
End If
If IsError(dice) Then Exit Function
If n = 1 Or f = 1 Then ‘handle trivial cases quickly
If IsMissing(s) Then
ReDim rv(0 To f – 1, 0 To 1)
For k = 0 To f – 1
rv(k, 0) = n + k
rv(k, 1) = CDbl(f) ^ -1
Next k
dice = rv
Else
dice = CDbl(f) ^ -1
End If
Exit Function
Else
prv = dice(n – 1, f)
lo = 0
hi = 0
u = n * (f – 1)
ReDim rv(0 To u, 0 To 1)
For k = 0 To Int(u / 2) ‘use symmetry
rv(k, 0) = n + k
rv(k, 1) = 0
For j = lo To hi
rv(k, 1) = rv(k, 1) + prv(j, 1)
Next j
rv(k, 1) = rv(k, 1) / CDbl(f)
rv(u – k, 0) = n + u – k
rv(u – k, 1) = rv(k, 1)
hi = hi + 1
If hi – lo = f Then lo = lo + 1
Next k
If IsMissing(s) Then dice = rv Else dice = rv(s – n, 2)
End If
End Function
This particular problem can be solved using the following array formula calling this function as a udf.
=MMULT(TRANSPOSE(INDEX(dice(9,4),0,2)),
MMULT((INDEX(dice(9,4),0,1)>TRANSPOSE(INDEX(dice(6,6),0,1))),
INDEX(dice(6,6),0,2)))
or using a VBA procedure like
Sub test()
'shamelessly using the fact that Peter and Colin
'have the same maximum possible score of 36
Dim p As Variant, c As Variant, pu As Long, cu As Long
Dim s As Double, t As Double, k As Double, dt As Date
Dim n As Long
dt = Timer
For n = 1 To 10000
p = dice(9, 4)
pu = UBound(p, 1)
c = dice(6, 6)
cu = UBound(c, 1)
t = 1 - c(cu, 1)
s = p(pu, 1) * t
Do While Sgn(pu) = 1
pu = pu - 1
cu = cu - 1
t = t - c(cu, 1)
s = s + p(pu, 1) * t
Loop
Next n
Debug.Print s; CDbl(Timer - dt)
End Sub
The ten thousand iterations run time averages about 2 seconds, implying each iteration's run time is .2 milliseconds.
fzz beats me to it! I thought that computing the distribution of rolls of a sum of dice might be more efficient, and wrote a recursion-based solution in C#, which runs in under a second:
http://clear-lines.com/blog/post/Euler-Problem-205.aspx
That being said, I wonder if the recursive approach is not more sensitive to rounding errors…
In any case, thanks for pointing this problem out, it was really fun!
Mathias
Mathias – if you believe the recursive approach is more sensitive to rounding errors, make the distribution one of counts rather than probabilities, i.e., skip dividing by CDbl(f) in my function until finishing the initial pass (would need a static boolean variable to detect this). No rounding errors with integers, minimal/unavoidable rounding errors dividing partial counts by total count.
Hi Mathias –
Just to correct your record at your website, the full solution above runs in 6/10ths of a second. Plenty fast enough. It was the brute force way I discarded than took 6 and a half minutes.
…mrt
Michael,
I bow my head in shame, and have set the record straight! I confess – I read your post a bit too fast. So here it is:
Dim i As Integer
For i = 1 To 100
MsgBox (“Next Time, I Will Read Posts Carefully Before Making Comments.”)
Next
End Sub
BTW I like your usage of Modulus, it’s pretty clever.
Mathias
Michael – your brute force approach run time is only a bit less than the run time for 10,000 iterations of a recursive approach.
The point is that your approach uses 6 * 6^6 + 9 * 4^9 + (36 – 9 + 1) * (36 – 6 + 1) = 2,640,100 iterations in different loops while recursion takes 4 + (9 – 1) * 4 / 2 + 6 + (6 – 1) * 6 / 2 + 36 – 9 + 1 = 69 iterations and 9 – 1 + 6 – 1 = 13 recursive function calls.
The two biggest time wasters in brute force are in calculating every entry in CC() and PP() arrays rather than using symmetry (for 2 6-sided dice the number of rolls producing 2 is the same as the number of rolls producing 12 and similarly for 3 and 11, 4 and 10, 5 and 9, 6 and 8), and effectively calculating the crossproduct of CC() and PP() rather than using something like
For P = 36 To 9 Step -1 ‘Peter’s SCORES in reverse order
CR = CR – CC(P) ‘now ttl # Colin possible rolls w/scores LT P
Count = Count + CR * PP(P)
Next P
Sorry. Make the 10,000 just 1,000.
Hi fzz –
But it’s a small brute…the big brute got sent home ;-)
What is Sgn(pu) ?
Mathias –
No problem, really. Thank you for the comments, and the compliment. Welcome to the Euler conversations.
…mrt
Sgn is a built-in VBA function, so a complete description is available in online help.
It’s not the unnecessarily large number of iterations themselves that cause the performance hit, it’s the statement
TEMP = -Int(-N / faces ^ C) Mod faces ‘where faces is actually a numeric constant, 6 or 9
that kills performance. When working with programming languages that don’t provide common subexpression elimination, such as VBA, you need to eliminate them yourself. Which means you should have used a Double-type variable to hold the value of faces ^ C, and used that variable rather than recalculating it in every iteration.
A couple million unnecessary exponentiation operations tends to affect execution time.
It seems I’m a bit late to this thread. My code runs in 0.06 seconds, significantly slower than fzz’s, so I won’t post it. Very impressive, fzz. I’m stepping through your dice function, trying very hard to understand it. I have a feeling that this approach may be the key to solving several other problems, notably #249. I’m not very good at dynamic programming, though I’m trying. Any help you could give with the logic of the function would be greatly appreciated. Thanks. By the way fzz, do you have a user name / preferred language on PE?
-Josh
Since I have no idea what PE is, I guess the answer is no.
My dice function is based on Lisp code a saw long ago. With one die, the score is discretely uniformly distributed with probability equal to 1 / number of faces. For 2 dice, the lowest score is 2 and the highest is 2 * number of faces, and with n+1 dice there are (faces – 1) more possible scores than with just n dice, with the lowest possible score going from n to n+1 and the highest possible score going from n*faces to (n+1)*faces. So the number of possible scores from rolling n dice with f faces is n*(f-1)+1, and the mean score is n*f/2.
Let S(n,f) represents the random variable of the score from rolling n dice with f faces with sequential scores 1 to f. Then S(n+1,f) = S(n,f) + S(1,f), so Pr[S(n+1,f) = s] = sum(Pr[S(n,f) = s-k] * Pr[S(1,f) = k], k = 1..f). That establishes the recursive relationship.
This is easier to see with formulas. Enter the following.
B1:0
C1:=B1+1
Fill C1 right into D1:H1.
B2:=SUM(B3:B65536)
Fill B2 right into C2:H2.
A9:0
A10:=A9+1
Fill A10 down into A11:A46.
B9:1
B10:0
C9:=IF(VAR(SIGN($A9-C$1+1),SIGN(C$1*6-$A9+1),1)=0,SUM(B3:B8)/6,0)
Fill C9 right into D9:H9.
Fill C9:H9 down into C10:H10.
Fill B10:H10 down into B11:H46.
[The VAR(SIGN(..),SIGN(..),1) but is equivalent to C1 less than or equal to A9 and A9 less than or equal to C1*6, but without the hassle of this blog fubarring angle brackets.]
The resulting table in A9:H46 shows the score probabilities for 0 to 6 dice in columns B through H, respectively, for the scores shown in column A. The sums in row 2 confirm that individual score probabilities sum to 1. The lo and hi variables in my dice function eliminate the need to store unnecessary blank or zero entries in the probability table array for n-1 dice when calculating probabilities for n dice.
I believe “PE” stands for Project Euler, the source for these problems:
http://projecteuler.net/
Screwed up the generalized mean. It’s n*(f+1)/2.
BTW Michael, you seem to get dangerously close from finishing the problems, so I thought I would point out that The Daily WTF has a new section, “Programming Praxis”, which could give you some extra miles to code, if they keep it going ;)
http://thedailywtf.com/Series/Programming_Praxis.aspx
fzz –
Thank you. So obvious in the telling, but I never would have thought of it. I’ll make the changes and see what the differences are.
Mathias –
No, not close. I just went over the 100-problem mark a few weeks ago. That qualified me as a novice, by the way the site assigns grades. Euler seems to be adding a problem a week, and my solution rate is falling off. Trends are going the wrong way ;-)
…mrt
fzz,
You’ve always been so quick and helpful with the Project Euler problems, even posting solutions, that I just assumed you were a member. My bad.
Thanks a lot for your explanation above. It’s a huge help.
-Josh
fzz –
Your exponential advice, considering we are talking 10ths of a second, was dramatic.
Old Time: 0.578125
New Time: 0.28125
With this as a prelude:
Dim Power6(0 To 5) As Double
Dim Power4(0 To 8) As Double
Power6(0) = 1
Power6(1) = 6
For N = 2 To 5
Power6(N) = 6 ^ N
Next N
Power4(0) = 1
Power4(1) = 4
For N = 2 To 8
Power4(N) = 4 ^ N
Next N
The 2 affected lines now look like this:
TEMP = -Int(-N / Power6(C)) Mod 6
TEMP = -Int(-N / Power4(P)) Mod 4
Thanks again.
…mrt
Well I just learned something: an 8 next to a close paren makes this 8)
Josh –
fzz should join, only he has to do the problems from the high end down to give us a chance ;-)
…mrt wondering where the smiley gouge is
Michael,
I agree – fzz should join. How about it fzz, got some extra time on your hands?
-Josh
For a Excel solution:
Project Euler – Problem 205
http://www.tushar-mehta.com/misc_tutorials/project_euler/euler205.html
Tushar does this by generating all possible permutations?! One of the goals of Project Euler is solving the problems efficiently. Generating all dice permutations is completely unnecessary. Brute force taken to to the extreme.
Answering my own question, the Smiley gouge is here:
http://codex.wordpress.org/Using_Smilies
…mrt :roll:
How about instead using the Troll probability calculator?
There you simply say count (sum 9d4) > (sum 6d6) and press a button!
Google “Troll dice roller” and you’ll find its webpage easily.