Calculate the exact result of a complex throw of two D30

Well, this has been bothering me for several years. If you sucked statistics and advanced math at school, turn around now. Too late.

Good. Take a deep breath. Here are the rules. Take two thirty-sided dice (yes, they exist ) and scan them at the same time.

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

If you write some code (see below), roll these bones several million times, and you will calculate how often you get each number as the final result, you get a curve that is quite flat to the left of 1, about 45 degrees from 1 to 60 and flat above 60. The probability of rolling 30.5 or higher is more than 50%, for rolling it is better 18 - 80%, and for rolling it is better than 0 - 97%.

Now the question is: is it possible to write a program to calculate the exact value of f (x), that is, the probability of overturning a certain value?

Reference Information. For our role-playing game "Jungle of Stars", we were looking for a way to keep random events under control. The rules above guarantee a much more stable result for what you are trying :)

For geeks around, 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