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 6^{2} different outcomes. Rolling 6 dice will have 6^{6} outcomes, or 46,656 different rolls.

Peter has 4^{9} 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 6^{6}), 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 6^{0} is 1, and 6^{1} 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 66^{th} 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 66^{th} 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 (6^{6})*(4^{9}),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 6^{ths} 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/10thsof a second. Plenty fast enough. It was the brute force way I discarded than took 6 and a half minutes.…mrtMichael,

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.

…mrtSgn 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 ;-)

…mrtfzz,

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.

…mrtWell 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 ;-)

…mrtwondering where the smiley gouge isMichael,

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.