Euler Problem 205

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.

  1. C = 0, Die 1:
    • INT(-66/6^0) = INT(-66/1) = -66
    • –66 = 66
    • 66 Mod 6 = 0, Return 6
  2. C = 1, Die 2:
    • INT(-66/6^1) = INT(-66/6) = -11
    • –11 = 11
    • 11 Mod 6 = 5, Return 5
  3. C = 2, Die 3:
    • INT(-66/6^2) = INT(-66/36) = INT(-1.83333) = -2
    • –2 = 2
    • 2 Mod 6 = 2, Return 2
  4. 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.

Sub Problem_205()
   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 P = 1 to 4^9
   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

Posted in Uncategorized

24 thoughts on “Euler Problem 205

  1. 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.

  2. Recursion is the ideal way to calculate the probability distributions for the different sets of dice.

    Function dice(n As Long, f As Long, Optional s As Variant) As Variant
      ‘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.

  3. 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

  4. 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.

  5. 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

  6. 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:

    Sub MeaCulpa()
        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

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

    CR = 6 ^ 6  ‘initially ttl # Colin possible rolls

    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
  8. 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

  9. 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.

  10. 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

  11. 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.

  12. 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

  13. 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

  14. 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

  15. 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

  16. 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.

  17. 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.

Leave a Reply

Your email address will not be published. Required fields are marked *