Calculate exact result of complex throw of two D30

Question

Okay, this bugged me for several years, now. If you sucked in statistics and higher math at school, turn away, now. Too late.

Okay. Take a deep breath. Here are the rules. Take two thirty sided dice (yes, they do exist) and roll them simultaneously.

• Add the two numbers
• If both dice show <= 5 or >= 26, throw again and add the result to what you have
• If one is <= 5 and the other >= 26, throw again and subtract the result from what you have
• Repeat until either is > 5 and < 26!

If you write some code (see below), roll those dice a few million times and you count how often you receive each number as the final result, you get a curve that is pretty flat left of 1, around 45° degrees between 1 and 60 and flat above 60. The chance to roll 30.5 or better is greater than 50%, to roll better than 18 is 80% and to roll better than 0 is 97%.

Now the question: Is it possible to write a program to calculate the exact value f(x), i.e. the probability to roll a certain value?

Background: For our role playing game "Jungle of Stars" we looked for a way to keep random events in check. The rules above guarantee a much more stable outcome for something you try :)

For the geeks around, the code in Python:

import random
import sys

def OW60 ():
    """Do an open throw with a "60" sided dice"""
    val = 0
    sign = 1

    while 1:
        r1 = random.randint (1, 30)
        r2 = random.randint (1, 30)

        #print r1,r2
        val = val + sign * (r1 + r2)
        islow = 0
        ishigh = 0
        if r1 <= 5:
            islow += 1
        elif r1 >= 26:
            ishigh += 1
        if r2 <= 5:
            islow += 1
        elif r2 >= 26:
            ishigh += 1

        if islow == 2 or ishigh == 2:
            sign = 1
        elif islow == 1 and ishigh == 1:
            sign = -1
        else:
            break

        #print sign

    #print val
    return val

result = [0] * 2000
N = 100000
for i in range(N):
    r = OW60()
    x = r+1000
    if x < 0:
        print "Too low:",r
    if i % 1000 == 0:
        sys.stderr.write('%d\n' % i)
    result[x] += 1

i = 0
while result[i] == 0:
    i += 1

j = len(result) - 1
while result[j] == 0:
    j -= 1

pSum = 0
# Lower Probability: The probability to throw this or less
# Higher Probability: The probability to throw this or higher
print "Result;Absolut Count;Probability;Lower Probability;Rel. Lower Probability;Higher Probability;Rel. Higher Probability;"
while i <= j:
    pSum += result[i]
    print '%d;%d;%.10f;%d;%.10f;%d;%.10f' % (i-1000, result[i], (float(result[i])/N), pSum, (float(pSum)/N), N-pSum, (float(N-pSum)/N))
    i += 1

Answer 1

Well, let's see. The second throw (which will sometimes be added or subtracted to the first roll) has a nice easily predictable bell curve around 31. The first roll, of course, is the problem.

For the first roll, we have 900 possible combinations.

• 50 combinations result in adding the second roll.
• 25 combinations result in subtracting the second roll.
• Leaving 825 combinations which match the bell curve of the second roll.

The subtracting set (pre-subtraction) will form a bell curve in the range (27..35). The lower half of the adding set will form a bell curve in the range (2..10), while the upper half will form a bell curve in the range (52...60)

My probablity is a bit rusty, so I can't figure the exact values for you, but it should be clear that these lead to predictable values.

Answer 2

Compound unbounded probability is... non-trivial. I was going to tackle the problem the same way as James Curran, but then I saw from your source code that there could be a third set of rolls, and a fourth, and so on. The problem is solvable, but far beyond most die rolling simulators.

Is there any particular reason that you need a random range from -Inf to +Inf with such a complex curve around 1-60? Why is the bell curve of 2D30 not acceptable? If you explain your requirements, it is likely someone could provide a simpler and more bounded algorithm.