Calculating odds of rolling a given total with a specified number of dice

[Dave Britten]

I could swear I had a direct formula for this at one point, but I'm totally drawing a blank.

Suppose I've got n standard 6-sided dice, numbered 1 through 6. Is there a formula for calculating the frequency (i.e. number of different permutations) of rolling a given total? For example, what are the odds of rolling a total of 22 on 9 dice? (It's about 1.45% in this case.)

The brute-force approach is really simple on SQL Server (code below), and runs all ~10 million permutations in about 1-2 seconds on the server I've got access to, but my HP 48 probably won't be so quick.

Code:

WITH die AS (
    SELECT v FROM (VALUES (1), (2), (3), (4), (5), (6)) v(v)
),
freq AS (
    SELECT
        d1.v + d2.v + d3.v + d4.v + d5.v + d6.v + d7.v + d8.v + d9.v AS total,
        COUNT(*) AS freq
    FROM die d1
        CROSS JOIN die d2
        CROSS JOIN die d3
        CROSS JOIN die d4
        CROSS JOIN die d5
        CROSS JOIN die d6
        CROSS JOIN die d7
        CROSS JOIN die d8
        CROSS JOIN die d9
    GROUP BY d1.v + d2.v + d3.v + d4.v + d5.v + d6.v + d7.v + d8.v + d9.v
)
SELECT
    total,
    freq,
    CAST(freq AS float) / CAST(SUM(freq) OVER (PARTITION BY NULL) AS float) AS pct,
    SUM(freq) OVER (PARTITION BY NULL) AS perms
FROM freq
ORDER BY total

And for extra fun, I'm also wondering if such a formula could be extended to rolling a pool of non-standard, mixed dice, e.g. maybe three of them numbered 0, 0, 0, 1, 2, 3, two numbered 0, 0, 1, 2, 2, 3, and one with 0, 0, 1, 2, 3, 4. The brute-force approach can be adapted for this trivially.

The motivation behind this would be in case I'm playing a board game and want to quickly math out my odds of success before diving into a potentially pivotal move.

[Thomas Klemm]

(06-15-2015 03:25 PM)Dave Britten Wrote:  For example, what are the odds of rolling a total of 22 on 9 dice? (It's about 1.45% in this case.)

It's the coefficient of \(x^{22}\) in:

\(
(\frac{x^6}{6}+\frac{x^5}{6}+\frac{x^4}{6}+\frac{x^3}{6}+\frac{x^2}{6}+\frac{x}{​6})^9
\)

The exact value is:

\(\frac{16211}{1119744}\)

Cheers
Thomas

[Dave Britten]

(06-15-2015 06:59 PM)Thomas Klemm Wrote:  

(06-15-2015 03:25 PM)Dave Britten Wrote:  For example, what are the odds of rolling a total of 22 on 9 dice? (It's about 1.45% in this case.)

It's the coefficient of \(x^{22}\) in:

\(
(\frac{x^6}{6}+\frac{x^5}{6}+\frac{x^4}{6}+\frac{x^3}{6}+\frac{x^2}{6}+\frac{x}{​6})^9
\)

The exact value is:

\(\frac{16211}{1119744}\)

Cheers
Thomas

That makes sense. What's the practical way to do a multinomial expansion like that, aside from feeding it through a CAS and letting it sum like terms? Seems like you'd have to find all permutations of 9 integers between 1 and 6 that add up to, e.g., 22, then you're back to square one.

[Thomas Klemm]

(06-15-2015 08:14 PM)Dave Britten Wrote:  What's the practical way to do a multinomial expansion like that, aside from feeding it through a CAS and letting it sum like terms?

You could use matrix multiplication to multiply two polynomials:

\(\begin{bmatrix}
1\\
1\\
1\\
1\\
1\\
1
\end{bmatrix}\cdot\begin{bmatrix}
1 & 1 & 1 & 1 & 1 & 1
\end{bmatrix}=\begin{bmatrix}
1 & 1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 & 1 \\
\end{bmatrix}\)

At least the HP-42S allows to grow the matrix at the bottom:

\(\begin{bmatrix}
1 & 1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 & 1 \\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 \\
\end{bmatrix}\)

Then you can transpose it:

\(\begin{bmatrix}
1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
\end{bmatrix}\)

Now when you shrink it just by the size of the 2nd matrix we get this:

\(\begin{bmatrix}
1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 \\
0 & 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 \\
\end{bmatrix}\)

To add up all the columns we can use again matrix multiplication:

\(\begin{bmatrix}
1 & 1 & 1 & 1 & 1 & 1
\end{bmatrix}\cdot\begin{bmatrix}
1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 \\
0 & 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 \\
\end{bmatrix}=\begin{bmatrix}
1 & 2 & 3 & 4 & 5 & 6 & 5 & 4 & 3 & 2 & 1
\end{bmatrix}\)

Now you can repeat this to calculate \(\begin{bmatrix}
1 & 1 & 1 & 1 & 1 & 1
\end{bmatrix}^4\) and \(\begin{bmatrix}
1 & 1 & 1 & 1 & 1 & 1
\end{bmatrix}^8\).

A final multiplication by \(\begin{bmatrix}
1 & 1 & 1 & 1 & 1 & 1
\end{bmatrix}\) will give you the desired list of coefficients:

CoefficientList[(1 + x + x^2 + x^3 + x^4 + x^5)^9, x]

Cheers
Thomas

[Thomas Klemm]

It's much easier to use lists:

Code:

\<< DUP
    1 5 START
        { 0 } SWAP +
        SWAP 0 +
        OVER ADD SWAP
    NEXT
    DROP
\>>

Start with this list:

{ 1 1 1 1 1 1 }

Repeat the program above 8 times.
To get the coefficient of \(x^{22}\) use:

14 GET

That's because 14 = 22 - 9 + 1 and the list is symmetric.

You should get: 145,899
Divide this by 6⁹ = 10,077,696 to get the desired probability.

HTH
Thomas

[Dave Britten]

(06-15-2015 11:52 PM)Thomas Klemm Wrote:  It's much easier to use lists:

Code:

\<< DUP
    1 5 START
        { 0 } SWAP +
        SWAP 0 +
        OVER ADD SWAP
    NEXT
    DROP
\>>

Start with this list:

{ 1 1 1 1 1 1 }

Repeat the program above 8 times.
To get the coefficient of \(x^{22}\) use:

14 GET

That's because 14 = 22 - 9 + 1 and the list is symmetric.

HTH
Thomas

Thanks Thomas, that looks like a good trick for dealing with pools of standard six-sided dice. Not sure if that can be applied to multiplication of arbitrary polynomials, i.e. pools of non-uniform, non-standard dice with any number of sides. I may end up having to use the matrix approach there and just multiply the polynomials, collect terms, etc.

[Thomas Klemm]

Another approach could be to use:

\(x^5+x^4+x^3+x^2+x+1=\frac{x^6-1}{x-1}\) if \(x\ne 0\).

Now we can use the binomial formulae to calculate the coefficients of \(x^k\) for both expressions \((x^6-1)^9\) and \((x-1)^9\) and just have to implement polynomial division.

Cheers
Thomas

[Thomas Klemm]

(06-15-2015 03:25 PM)Dave Britten Wrote:  And for extra fun, I'm also wondering if such a formula could be extended to rolling a pool of non-standard, mixed dice, e.g. maybe three of them numbered 0, 0, 0, 1, 2, 3, two numbered 0, 0, 1, 2, 2, 3, and one with 0, 0, 1, 2, 3, 4.

CoefficientList[(3+x+x^2+x^3)^3 (2+x+2 x^2+x^3)^2 (2+x+x^2+x^3+x^4), x]

In this case it appears that 7 is the most probable total.

Cheers
Thomas

[Dave Britten]

Looks like there's a ton of options for multinomial expansion if you're raising a single multinomial to a power. I'll probably have to take the direct polynomial multiplication approach, though, since I want to handle arbitrary numbers of dice which may be non-uniform and non-standard. That should be a hell of a lot more efficient than summing up every permutation, though.

[Thomas Klemm]

(06-16-2015 12:31 PM)Dave Britten Wrote:  I'll probably have to take the direct polynomial multiplication approach, though, since I want to handle arbitrary numbers of dice which may be non-uniform and non-standard.

You might consider using FFT as it is a built-in function.
Polynomial Multiplication and Fast Fourier Transform

Cheers
Thomas

[Thomas Klemm]

Just tested it with HP-48GX.
You start with a 64-dimensional vector:

[0, 0, … , 0 , 0, 1, 1, 1, 1, 1, 1]

Run the discrete Fourier transform:

FFT

Transmogrify the result to a list of 64 complex numbers:

OBJ→
OBJ→
DROP
→LIST

{ (6,0) (5.5701807… }

Apply the 9th power (point-wise):

9 y^x

{ (10077696,0) (-8… }

Transmogrify back to a vector:

OBJ→
{ }
+
→ARRY

Compute the inverse Fourier transform:

IFFT

[ (.00000034,0) (.… ]

Drawback is that you have now an array of complex numbers that are nearly real.
Thus transmogrify back to a list and run:

ABS
0 RND

And finally you have the list of coefficients:

{ 0 0 0 0 0 0 0 0 0
0 1 9 45 165 495
1287 2994 6354
12465 22825 39303

Of course the list is a little longer than just that.

Cheers
Thomas

[Dave Britten]

(06-16-2015 02:10 PM)Thomas Klemm Wrote:  Just tested it with HP-48GX.
You start with a 64-dimensional vector:

[0, 0, … , 0 , 0, 1, 1, 1, 1, 1, 1]

Run the discrete Fourier transform:

FFT

Transmogrify the result to a list of 64 complex numbers:

OBJ→
OBJ→
DROP
→LIST

{ (6,0) (5.5701807… }

Apply the 9th power (point-wise):

9 y^x

{ (10077696,0) (-8… }

Transmogrify back to a vector:

OBJ→
{ }
+
→ARRY

Compute the inverse Fourier transform:

IFFT

[ (.00000034,0) (.… ]

Drawback is that you have now an array of complex numbers that are nearly real.
Thus transmogrify back to a list and run:

ABS
0 RND

And finally you have the list of coefficients:

{ 0 0 0 0 0 0 0 0 0
0 1 9 45 165 495
1287 2994 6354
12465 22825 39303

Of course the list is a little longer than just that.

Cheers
Thomas

Very interesting stuff. I have to admit, I haven't done anything with FFTs before (though I have a general idea of what they're used for), so this is a bit out of my depth.

However, a simple C# version on my desktop gives results for 9 standard 6-sided (1-6) dice nearly instantaneously. The code is below; I've avoided heavily optimizing the array use to keep the indexing simple and make sure it's working correctly first, so there's a lot of zero padding for the lower powers of x. Note that results aren't necessarily symmetrical for for "unusual" die pools, so it can't just compute up to the halfway mark.

I should be able to do a C version for my 200LX, though the input parsing and mallocing will be a bit more labor-intensive. Everything else is just simple looping and arithmetic. I have a feeling trying to handle this with nested DOLISTs on a 48 would be fairly slow, but I might give it a shot.

Code:

    class Program {
        static void Main(string[] args) {
            List<int[]> dice = new List<int[]>(); //Input dice, transformed to polynomial coefficients (e.g. 1, 1, 1, 1, 1, 1)
            int[] d = null; //Single input die, for parsing
            int dq = 0; //Number of dice
            int[] result = null; //Accumulated result
            int[] lastresult = null; //Previous result, used in repeated multiplication
            int[] current = null; //Current die being multiplied with lastresult
            string line; //User input
            int q = 0; //Quantity (multiplicity) of input die
            List<int> inputs; //Parsed user input, with values of die faces (e.g. 1, 2, 3, 4, 5, 6)

            double p; //Total permutations

            //Input parsing loop
            while (true) {
                Console.Write("Die {0}? ", dq + 1);
                line = Console.ReadLine();
                if (String.IsNullOrWhiteSpace(line)) {
                    break;
                }
                Console.Write("Die {0} quantity? [1] ", dq + 1);
                if (!Int32.TryParse(Console.ReadLine(), out q)) {
                    q = 1;
                }
                inputs = line.Split(new char[]{' ', ',', ';'}, StringSplitOptions.RemoveEmptyEntries).Select(x => Int32.Parse(x.Trim())).ToList<int>();
                d = new int[inputs.Max() + 1]; //There's a zero element that remains 0, just to keep the indexing straight-forward. If I do a C version, it'll probably store the array length.
                for (int i = 0; i < inputs.Count; i++) {
                    d[inputs[i]]++;
                }
                for (int i = 0; i < q; i++) {
                    dq++;
                    dice.Add(d);
                }
            };

            //Calculation loop(s)
            lastresult = dice[0];
            p = (double)lastresult.Length - 1.0D; //Don't count the zero element.
            for (int i = 1; i < dq; i++) { //Loop through the dice.
                current = dice[i];
                p *= (double)current.Length - 1.0D; //Don't count the zero element.
                result = new int[lastresult.Length + current.Length - 1];
                for (int x = 0; x < lastresult.Length; x++) {
                    if (lastresult[x] == 0) {
                        continue;
                    }
                    for (int y = 0; y < current.Length; y++) {
                        result[x + y] += lastresult[x] * current[y];
                    }
                }
                lastresult = result;
            }

            //Result display loop
            bool skipzero = true;
            Console.WriteLine("Permutations: {0:n0}", p);
            Console.WriteLine("Total\tFreq\tProb");
            for (int i = 0; i < lastresult.Length; i++) {
                if (lastresult[i] == 0 && skipzero) {
                    continue;
                }
                skipzero = false;
                Console.WriteLine("{0}\t{1}\t{2:p4}", i, lastresult[i], (double)lastresult[i] / p);
            }
        }
    }

Usage:
Run the program, enter the values of the faces of each die, separated with spaces, commas, or semicolons. For example, "1 2 3 4 5 6", or "0 0 0 1 2 3". Enter the quantity after each given die, or press enter for the default of 1. When finished, enter a blank line when prompted for the next die faces.

Sample output:

Code:

Die 1? 0 0 0 1 2 3
Die 1 quantity? [1] 3
Die 4? 0 0 1 1 2 3
Die 4 quantity? [1] 2
Die 6? 0 0 1 2 3 4
Die 6 quantity? [1]
Die 7?
Permutations: 46,656
Total   Freq    Prob
0       216     0.4630 %
1       756     1.6204 %
2       1584    3.3951 %
3       2816    6.0357 %
4       4166    8.9292 %
5       5293    11.3447 %
6       6038    12.9415 %
7       6163    13.2094 %
8       5662    12.1356 %
9       4737    10.1530 %
10      3602    7.7203 %
11      2475    5.3048 %
12      1546    3.3136 %
13      871     1.8669 %
14      434     0.9302 %
15      193     0.4137 %
16      74      0.1586 %
17      23      0.0493 %
18      6       0.0129 %
19      1       0.0021 %

[Thomas Klemm]

(06-16-2015 03:00 PM)Dave Britten Wrote:  Very interesting stuff. I have to admit, I haven't done anything with FFTs before (though I have a general idea of what they're used for), so this is a bit out of my depth.

This is the relevant Convolution Theorem:

It's similar to how we avoid the multiplication with: \(\log(a\cdot b) = \log(a) + \log(b)\).

Code:

%%HP: T(3)A(D)F(.);
DIR
  A
[ 3 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
  B
[ 2 1 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
  C
[ 2 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
  TO
    \<< FFT 
        OBJ\-> OBJ\-> DROP
        \->LIST
    \>>
  FRO
    \<< OBJ\-> { } + \->ARRY
        IFFT
        OBJ\-> OBJ\-> DROP
        \->LIST
        ABS
        0 RND
    \>>
END

With these two programs you can now run:

A TO 3 y^x
B TO x² ×
C TO ×
FRO

{ 216 540 1314 2393
3618 4880 5774 6112
5878 5116 4022 2886
1870 1080 562 256
98 32 8 1 0 0 0 0 0
0 0 0 0 0 0 0 }

HTH
Thomas

[Dave Britten]

I've been playing around with Pythonista, a rather nice Python IDE, on my iPad. Here's a Python version of the program, which I imagine will work on plenty more than just Pythonista:

Code:

dice = [] #Input dice as polynomial coefficients    
dice_quantity = 0 #Total number of dice

while True:
    rawdie = [] #List of face values
    qty = 1 #Quantity of a single die spec

    inp = raw_input("Die " + str(dice_quantity + 1) + "? ")
    if inp == "":
        break
    for face in inp.split(" "):
        rawdie.append(int(face))
    qty = raw_input("Quantity? [1] ")

    if qty == "":
        qty = 1
    else:
        qty = int(qty)

    die = [0] * (max(rawdie) + 1)
    for i in rawdie:
        die[i] += 1

    for i in range(qty):
        dice.append(die)

    print(str(rawdie) + " x " + str(qty))
    dice_quantity = dice_quantity + qty

lastresult = dice[0]

for current in dice[1:]:
    result = [0] * (len(lastresult) + len(current) - 1)
    for x in range(len(lastresult)):
        if lastresult[x] == 0:
            continue
        for y in range(len(current)):
            result[x + y] += lastresult[x] * current[y]
    lastresult = result

skipzero = True
p = 0.0
ft = 0.0
for i in lastresult:
    p += i

print "Permutations: {:.0f}".format(p)
print "{:<6} {:<12} {:>10} {:>10} {:>10}".format("Total","Freq","Prob", "p>=", "p<=")
for i in range(len(lastresult)):
    ft += lastresult[i]
    if lastresult[i] == 0 and skipzero:
        continue
    skipzero = False
    print "{:<6} {:<12} {:>10.3%} {:>10.3%} {:>10.3%}".format(i, lastresult[i], lastresult[i] / p, (p - ft + lastresult[i]) / p, ft / p)