# Custom dice roll probability formula

I want to make an application that gives you the probability of reaching a certain threshold with a number of rolls of special dice. Just to clarify, I am not asking for help on programming, I need a formula to program from.

The variables are

• d: the number of dice, all with the sides [1, 2, 2, 3, 3, 4]

• x: the number the roll needs to reach

I need to work out the probability of rolling d dice (of the special type noted above) and getting a total equal to or greater than x

• I'd try simulation. Make a virtual die and roll $d$ of them $10000$ times. Then see how many times you get to that total.The alternative is just solving a generic equation which I'm sure you can do. – stuart stevenson Jan 1 '18 at 23:52

Let $p(t,d)$ be the probability that when you roll $d$ dice, the sum of their numbers is exactly $t.$

Then $p(t,d) = 0$ whenever $t < 0.$ Also $p(t,0) = 0$ when $t > 0,$ but $p(0,0) = 1.$

For $d \geq 0,$ $$p(t,d+1) = \frac16 p(t-1,d) + \frac13 p(t-2,d) + \frac13 p(t-3,d) + \frac16 p(t-4,d).$$

You can use that formula to find any value of $p(t,d)$ recursively. An efficient way to do this is to compute $p(t,1)$ for every relevant value of $t,$ then $p(t,2)$ for every relevant $t,$ and so forth until you have computed the necessary values of $p(t,d)$ in order to find the probability of rolling at least $x$; that probability is $$\sum_{t\geq x} p(t,d).$$

• Note that for $t \in \{d,d+1,\cdots, 4d\}$, $p(t,d)$ is the coefficient of $x^t$ in the polynomial $\phi^d(x)$ below. Exactly the same recursion appears, unsurprisingly, it is essentially a convolution. Do note that $p(t,d)$ is nonzero iff for $t \in \{d,d+1,\cdots, 4d\}$. – copper.hat Jan 3 '18 at 5:15

For small values of $d$ you could use the following generating function approach:

Let $\phi(x) = {x+ 2 x^2 + 2 x^3 + x^4 \over 6}$, then the probability that $d$ rolls gives exactly $m$ is the coefficient of $x^m$ in $\phi^d$, where $m=d,d+1,...,4d$.

To illustrate: If $d=7$ then $\phi^d(x) = \sum_{k=0}^{28} a_k x^k$ where $a={1 \over 6^7}(0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 14 , 98 , 455 , 1568 , 4256 , 9429 , 17446 , 27370 , 36771 , 42560 , 42560 , 36771 , 27370 , 17446 , 9429 , 4256 , 1568 , 455 , 98 , 14 , 1)$.

For example, the probability of getting a sum of exactly $26$ is the coefficient of $x^{26}$ which is ${98 \over 6^7}$.

The probability of the sum being $\ge x$ will be the sum of the coefficients of $x^k$ where $k \ge x$.

There are various ways of computing the polynomial coefficients depending on what you are trying to achieve.

# poly of phi*6, constant term last
f = [1 2 2 1 0];
# 7 rolls
d = 7;
# compute polynomial of (phi*6)^d
p = 1;
for i = 1:d
p = conv(p, f);
endfor
# scale result so the coefficients are probabilities
# constant term last
p = p/6^d;


If not, you could using something like Neville's algorithm to compute the interpolating polynomial of the points $\phi^d(0),\cdots, \phi^d(4d)$.

Another way would be to compute the coefficients by multiplying the polynomials $\phi^2 = \phi \cdot \phi$, $\phi^{k+1} = \phi^k \cdot \phi$. This is essentially the convolution of coefficients of the polynomials being multiplied.

Aside: My original Octave/Matlab code looked something like

f = [1 2 2 1 0];
r = roots(f);
d = 7;
v = poly(repmat(r',[1, d]));
v = v/6^d;


however, for $d \in \{5,6,7\}$, the poly call resulted in some small imaginary noise. The conv call avoids this entirely.

• How small? d will always be less than seven. – TheOnlyMrCat Jan 2 '18 at 8:46
• The call to poly produces an error with both Octave and MATLAB because the matrix argument is not square. – Fabio Somenzi Jan 2 '18 at 21:54
• @FabioSomenzi: My apologies, I will fix momentarily. – copper.hat Jan 2 '18 at 22:02
• @FabioSomenzi: Thanks for catching that, I gave up on my previous approach as it failed for $d=7$. There were a number of errors in my previous code, I should know not to wing it at this stage! – copper.hat Jan 2 '18 at 23:05
• Glad to be of some help! – Fabio Somenzi Jan 2 '18 at 23:13