First of all I'm not sure where this question should be posted. I'm asking if a statistics problem is NP-Complete and if not to solve it programmatically. I'm posting it here because the statistics problem is the center point.

I'm trying to find a better formula for solving a problem. The problem is: if I have 4d6 (4 ordinary 6 sided dice) and roll them all at once, remove a die with the lowest number (called "dropping"), then sum the remaining 3, what is the probability of each possible result? I know the answer is this:

Sum (Frequency): Probability
3   (1):         0.0007716049
4   (4):         0.0030864198
5   (10):        0.0077160494
6   (21):        0.0162037037
7   (38):        0.0293209877
8   (62):        0.0478395062
9   (91):        0.0702160494
10  (122):       0.0941358025
11  (148):       0.1141975309
12  (167):       0.1288580247
13  (172):       0.1327160494
14  (160):       0.1234567901
15  (131):       0.1010802469
16  (94):        0.0725308642
17  (54):        0.0416666667
18  (21):        0.0162037037

The average is 12.24 and standard deviation is 2.847.

I found the above answer by brute force and don't know how or if there is a formula for it. I suspect this problem is NP-Complete and therefore can only be solved by brute force. It might be possible to get all probabilities of 3d6 (3 normal 6 sided dice) then skew each of them upward. This would be faster than brute force because I have a fast formula when all dice are kept.

I programmed the formula for keeping all dice in college. I had asked my statistics professor about it and he found this page, which he then explained to me. There is a big performance difference between this formula and brute force: 50d6 took 20 seconds but 8d6 drop lowest crashes after 40 seconds (chrome runs out of memory).

Is this problem NP-Complete? If yes please provide a proof, if no please provide a non-brute force formula to solve it.

Note that I don't know much about NP-Complete so I might be thinking of NP, NP-Hard, or something else. The proof for NP-Completeness is useless to me the only reason why I ask for it is to prevent people from guessing. And please bare with me as it's been a long time since I worked on this: I don't remember statistics as well as I might need to solve this.

Ideally I'm looking for a more generic formula for X number of dice with Y sides when N of them are dropped but am starting with something much more simple.


I would also prefer the formula to output frequencies but it is acceptable to only output probabilities.

For those interested I have programmed whuber's answer in JavaScript on my GitHub (in this commit only the tests actually use the functions defined).

  • 1
    $\begingroup$ This is an interesting question. I think it should be on-topic here. Thank you for your consideration. $\endgroup$ Commented Dec 22, 2014 at 16:40
  • 1
    $\begingroup$ Although the setting is interesting, you haven't yet asked an answerable question: the idea of NP-completeness depends on having a class of problems, while you have described only one. Exactly how do you want it to generalize? Although you hint that the numbers of dice could vary, various additional options are possible and they might yield different answers: you could change the numbers of faces, the values on the faces, the numbers of dice, and the numbers of dropped dice, all in various ways with various relationships among them. $\endgroup$
    – whuber
    Commented Dec 22, 2014 at 16:49
  • 1
    $\begingroup$ @whuber She doesn't know any complexity theory but I think it's clear that she's asking after the family of problems generated by changing the number of dice. I also think I have an efficient algorithm for it. $\endgroup$
    – Andy Jones
    Commented Dec 22, 2014 at 17:04
  • 2
    $\begingroup$ @Andy I do see at the end she is asking for "a more generic formula for X number of dice with Y sides when N of them are dropped". $\endgroup$
    – whuber
    Commented Dec 22, 2014 at 17:34
  • $\begingroup$ @whuber Hah! Apparently not as clear as I thought then. Sorry, my bad. $\endgroup$
    – Andy Jones
    Commented Dec 22, 2014 at 18:09

4 Answers 4



Let there be $n=4$ dice each giving equal chances to the outcomes $1, 2, \ldots, d=6$. Let $K$ be the minimum of the values when all $n$ dice are independently thrown.

Consider the distribution of the sum of all $n$ values conditional on $K$. Let $X$ be this sum. The generating function for the number of ways to form any given value of $X$, given that the minimum is at least $k$, is

$$f_{(n,d,k)}(x) = x^k+x^{k+1} + \cdots + x^d = x^k\frac{1-x^{d-k+1}}{1-x}.\tag{1}$$

Since the dice are independent, the generating function for the number of ways to form values of $X$ where all $n$ dice show values of $k$ or greater is

$$f_{(n,d,k)}(x)^n = x^{kn}\left(\frac{1-x^{d-k+1}}{1-x}\right)^n.\tag{2}$$

This generating function includes terms for the events where $K$ exceeds $k$, so we need to subtract them off. Therefore the generating function for the number of ways to form values of $X$, given $K=k$, is

$$f_{(n,d,k)}(x)^n - f_{(n,d,k+1)}(x)^n.\tag{3}$$

Noting that the sum of the $n-1$ highest values is the sum of all values minus the smallest, equal to $X-K$. The generating function therefore needs to be divided by $k$. It becomes a probability generating function upon multiplying by the common chance of any combination of dice, $(1/d)^n$:

$$d^{-n}\sum_{k=1}^dx^{-k}\left(f_{(n,d,k)}(x)^n - f_{(n,d,k+1)}(x)^n\right).\tag{4}$$

Since all the polynomial products and powers can be computed in $O(n\log n)$ operations (they are convolutions and therefore can be carried out with the discrete Fast Fourier Transform), the total computational effort is $O(k\,n\log n)$. In particular, it is a polynomial time algorithm.


Let's work through the example in the question with $n=4$ and $d=6$.

Formula $(1)$ for the PGF of $X$ conditional on $K\ge k$ gives

$$\eqalign{ f_{(4,6,1)}(x) &= x+x^2+x^3+x^4+x^5+x^6 \\ f_{(4,6,2)}(x) &= x^2+x^3+x^4+x^5+x^6 \\ \ldots \\ f_{(4,6,5)}(x) &= x^5+x^6 \\ f_{(4,6,6)}(x) &= x^6 \\ f_{(4,6,7)}(x) &= 0. }$$

Raising them to the $n=4$ power as in formula $(2)$ produces

$$\eqalign{ f_{(4,6,1)}(x)^4 &= x^4 + 4x^5 + 10 x^6 + \cdots + 4x^{23} + x^{24} \\ f_{(4,6,2)}(x)^4 &= x^8 + 4x^9 + 10x^{10}+ \cdots + 4x^{23} + x^{24} \\ \ldots \\ f_{(4,6,5)}(x)^4 &=x^{20} + 4 x^{21} + 6 x^{22} + 4x^{23} +x^{24}\\ f_{(4,6,6)}(x)^4 &= x^{24}\\ f_{(4,6,7)}(x)^4 &= 0 }$$

Their successive differences in formula $(3)$ are

$$\eqalign{ f_{(4,6,1)}(x)^4 - f_{(4,6,2)}(x)^4 &= x^4 + 4x^5 + 10 x^6 + \cdots + 12 x^{18} + 4x^{19} \\ f_{(4,6,2)}(x)^4 - f_{(4,6,3)}(x)^4 &= x^8 + 4x^9 + 10x^{10} + \cdots + 4 x^{20} \\ \ldots \\ f_{(4,6,5)}(x)^4 - f_{(4,6,6)}(x)^4 &=x^{20} + 4 x^{21} + 6 x^{22} + 4x^{23} \\ f_{(4,6,6)}(x)^4 - f_{(4,6,7)}(x)^4 &= x^{24}. }$$

The resulting sum in formula $(4)$ is

$$6^{-4}\left(x^3 + 4x^4 + 10x^5 + 21x^6 + 38x^7 + 62x^8 + 91x^9 + 122x^{10} + 148x^{11} + \\167x^{12} + 172x^{13} + 160x^{14} + 131x^{15} + 94x^{16} + 54x^{17} + 21x^{18}\right).$$

For example, the chance that the top three dice sum to $14$ is the coefficient of $x^{14}$, equal to

$$6^{-4}\times 160 = 10/81 = 0.123\,456\,790\,123\,456\,\ldots.$$

It is in perfect agreement with the probabilities quoted in the question.

By the way, the mean (as calculated from this result) is $15869/1296 \approx 12.244598765\ldots$ and the standard deviation is $\sqrt{13\,612\,487/1\,679\,616}\approx 2.8468444$.

A similar (unoptimized) calculation for $n=400$ dice instead of $n=4$ took less than a half a second, supporting the contention that this is not a computationally demanding algorithm. Here is a plot of the main part of the distribution:


Since the minimum $K$ is highly likely to equal $1$ and the sum $X$ will be extremely close to having a Normal$(400\times 7/2, 400\times 35/12)$ distribution (whose mean is $1400$ and standard deviation is approximately $34.1565$), the mean must be extremely close to $1400-1=1399$ and the standard deviation extremely close to $34.16$. This nicely describes the plot, indicating it is likely correct. In fact, the exact calculation gives a mean of around $2.13\times 10^{-32}$ greater than $1399$ and a standard deviation around $1.24\times 10^{-31}$ less than $\sqrt{400\times 35/12}$.

  • 1
    $\begingroup$ Your answer is fast and is correct so I've marked it as the answer. Also in an edit I said it would also be nice to have frequencies if possible. For that you don't need to edit your answer since I can see that the 6^-4 multiplier is used to convert from frequency to probability. $\endgroup$
    – SkySpiral7
    Commented Dec 5, 2016 at 1:00
  • 1
    $\begingroup$ This is a really good answer. What modifications would have to be made to take only the highest 2 dice (or, the highest $h$ dice in general), instead of the highest 3? I assume it would modify equation 4 but I can't work it out or find a reference. $\endgroup$ Commented Sep 4, 2020 at 5:35

Edit: @SkySpiral has had trouble getting the below formula to work. I currently don't have time to work out what the issue is, so if you're reading this it's best to proceed under the assumption it's incorrect.

I'm not sure about the general problem with varying numbers of dice, sides, and drops, but I think I can see an efficient algorithm for the drop-1 case. The qualifier is that I'm not completely sure that it's correct, but right now I can't see any flaws.

Let's start by not dropping any dice. Suppose $X_n$ represents the $n$th die, and suppose $Y_n$ represents the sum of $n$ dice. Then

$$p(Y_n = a) = \sum_k p(Y_{n-1} = a - k)p(X_n=k)$$

Now suppose $Z_n$ is the sum of $n$ dice when one die is dropped. Then

$$p(Z_n = a) = p(\text{$n$th die is the smallest})p(Y_{n-1} = a) + \\ p(\text{$n$th die is not the smallest})\sum_k p(Z_{n-1} = a - k)p(X_n=k)$$

If we define $M_n$ to be distribution of the minimum of $n$ dies, then

$$p(Z_n = a) = p(X_n \leq M_{n-1})p(Y_{n-1} = a | X_n \leq M_{n-1}) + \\ p(X_n > M_{n-1})\sum_k p(Z_{n-1} = a - k)p(X_n=k | X_n > M_{n-1})$$

and we can calculate $M_n$ using

$$p(M_n = a) = p(X_n \leq M_{n-1})p(X_n = a |X_n \leq M_{n-1}) + p(X_n > M_{n-1})p(M_{n-1} = a|X_n > M_{n-1})$$

Anyway, together this all suggests a dynamic programming algorithm based on $Y_n, Z_n$ and $M_n$. Should be quadratic in $n$.

edit: A comment has been raised on how to calculate $p(X_n \leq M_{n-1})$. Since $X_n, M_{n-1}$ can each only take on one of six values, we can just sum over all possibilities:

$$p(X_n \leq M_{n-1}) = \sum_{a,b} p(X_n = a, M_{n-1} = b, a \leq b)$$

Similarly, $p(X_n = k | X_n > M_{n-1})$ can be calculated by applying Bayes rule then summing over the possible values of $X_n, M_{n-1}$.

  • 1
    $\begingroup$ +1 This looks correct and you said that's it's quadratic. But it's been a few years since I took statistics (I'm primarily a programmer). So I'd like to fully understand this before marking it as the answer. Also I see you have p(nth is the smallest die) does this include if nth is tied with the smallest? Such as rolling all 3s. $\endgroup$
    – SkySpiral7
    Commented Dec 26, 2014 at 16:21
  • $\begingroup$ Good catch. If the $n$th die rolled is the same as the current minimum, we can regard that die as the one to be dropped. In which case the distribution is $Y_{n-1}$. I've swapped some $(<)$s for $(\leq)$s to reflect this. $\endgroup$
    – Andy Jones
    Commented Dec 26, 2014 at 18:58
  • $\begingroup$ Thank you. If I understand this correctly I think your formulas are the answer. However I don't know how to calculate p(X(n) > M(n-1)) (or the negation of it) or p(X(n)=k|X(n) > M(n-1)) so I can't use this answer yet. I'll mark this as the answer but I'd like more information. Can you edit your answer to explain these or should I post it as another question? $\endgroup$
    – SkySpiral7
    Commented Dec 27, 2014 at 2:56
  • $\begingroup$ Edited my answer. $\endgroup$
    – Andy Jones
    Commented Dec 27, 2014 at 10:42
  • 1
    $\begingroup$ Sorry I know it's been a year and a half but I've finally gotten around to implementing this formula into code. However the p(Z(n)=a) formula appears incorrect. Suppose 2 dice with 2 sides (drop lowest), what are the chances of the result being 1? The chance of X(n) being the smallest or tied is 3/4 and p(Y(n-1)=1) is 1/2 so that Z(n) returns at least 3/8 even though the correct answer is 1/4. The Z formula looks correct to me and I don't know how to fix it. So if it's not too much to ask: what do you think? $\endgroup$
    – SkySpiral7
    Commented Jul 10, 2016 at 20:36

I have a reasonably efficient algorithm for this that, on testing, seems to match results of pure brute force while relying less heavily on enumerating all possibilities. It's actually more generalized than the above problem of 4d6, drop 1.

Some notation first: Let $X_NdY$ indicate that you are rolling $X$ dice with $Y$ faces (integer values $1$ to $Y$), and considering only the highest $N$ dice rolled. The output is a sequence of dice values, e.g. $4_3d6$ yields $3, 4, 5$ if you rolled $1, 3, 4, 5$ on the four dice. (Note that I'm calling it a "sequence," but the order is not important here, particularly since all we care about in the end is the sum of the sequence.)

The probability $P(X_NdY = S)$ (or more specifically, $P(4_3d6 = S)$) is a simplified version of the original problem, where we are only considering a specific set of dice, and not all possible sets that add up to a given sum.

Suppose $S$ has $k$ distinct values, $s_0, s_1, ..., s_k$, such that $s_i > s_{i+1}$, and each $s_i$ has a count of $c_i$. For example, if $S = 3, 4, 4, 5$, then $(s_0,c_0) = (5,1)$, $(s_1,c_1) = (4,2)$, and $(s_2,c_2) = (3,1)$.

You can calculate $P(X_NdY = S)$ in the following way:

$$ P(X_NdY = S) = \frac{ \left( \prod_{i=0}^{k-1} {X - \sum_{h=0}^{i-1} c_h \choose c_i} \right) \left( \sum_{j=0}^{X-N} { c_k+X-N \choose c_k+X-N-j} (s_k-1)^j \right)}{ Y^X }$$

That's pretty messy, I know.

The product expression $\prod_{i=0}^{k-1}$ is iterating through all but the lowest of the values in $S$, and calculating all the ways those values may be distributed among the dice. For $s_0$, that's just $X \choose c_i$, but for $s_1$, we have to remove the $c_0$ dice that have already been set aside for $s_0$, and likewise for $s_i$ you must remove $\sum_{h=0}^{i-1}c_h$.

The sum expression $\sum_{j=0}^{X-N}$ is iterating through all the possibilities of how many of the dropped dice were equal to $s_k$, since that affects the possible combinations for the un-dropped dice with $s_k$ as their value.

By example, let's consider $P[4_3d6=(5,4,4)]$:

$$ (s_1, c_1) = (5, 1) $$ $$ (s_2, c_2) = (4, 2) $$

So using the formula above:

$$ P[4_3d6=(5,4,4)] \\ = \frac{ {4 \choose 1} \left( {3 \choose 3} \cdot 3^0 + {3 \choose 2} \cdot 3^1 \right) }{ 6^4 } \\ = \frac{5}{162} = 0.0\overline{308641975}$$

The formula breaks down on a domain issue when $s_k=1$ and $j=0$ in the summation, leading to a first term of $0^0$, which is indeterminate and needs to be treated as $1$. In such a case, a summation is not actually necessary at all, and can be omitted, since all the dropped dice will also have a value of $s_k = 1$.

Now here's where I do need to rely on some brute force. The original problem was to calculate the probability of the sum being some value, and $X_NdY$ represents the individual dice left after dropping. This means you must add up the probabilities for all possible sequences $S$ (ignoring ordering) whose sum is the given value. Perhaps there is a formula to calculate this across all such values of $S$ at once, but I haven't even tried broaching that yet.

I've implemented this in Python first, and the above is an attempt to express it mathematically. My Python algorithm is accurate and reasonably efficient. There are some optimizations that could be made for the case of calculating the entire distribution of $\sum X_NdY$, and maybe I'll do that later.

  • $\begingroup$ As a programmer it might be easier for me to understand your Python code (although I've never used Python so it might be the same). Posting the code here is off topic but you could post a link to github etc. $\endgroup$
    – SkySpiral7
    Commented Oct 31, 2016 at 1:47
  • 1
    $\begingroup$ Your answer may be correct and it seems to reduce the complexity from O(Y^X) to O((Y+X-1)!/(X!*(Y-1)!)) but it still isn't as efficient as whuber's answer of O(c*X*log(X)). Thanks for your answer though +1. $\endgroup$
    – SkySpiral7
    Commented Nov 29, 2016 at 4:21

whuber's fast Fourier transform is fantastically efficient. However, I'm going to take another angle at the problem. Let's go back to the beginning:

Why might one think this problem is NP-complete?

If the result of the roll is allowed to depend on the entire sequence of individual dice, then each of the $d^n$ possible rolls must be considered, which results in an exponential number of evaluations. If we only consider the multiset of dice rolls (or equivalently, the possible sorted sequences that could result from the roll), then this reduces to $\binom{n + d - 1}{d - 1}$, which is still exponential if both the number of dice $n$ and die size $d$ are both allowed to increase linearly. If $d$ is fixed, then this becomes polynomial---though of order $d - 1$, which can be quite large if the dice have many faces, explode, or are themselves the result of other computations, e.g. a "die" with faces 3 to 18 representing a D&D ability score.

Of course, as whuber's answer shows, the problem is well within P. How can this be? We are not evaluating an arbitrary function over the dice rolls; rather, the action of dropping the lowest and taking the sum has structure that we can exploit using convolution and FFTs.

This raises the question: is there some more general type of structure that we can exploit in some way?

"Single-pass" functions over order statistics

Consider this procedure for evaulating the roll of a dice pool. I first encountered the seeds of this idea in this answer by Ilmari Karonen, and similar ideas were also used in this paper.

  • First, I tell you how many ones you rolled.
  • Then, I tell you how many twos you rolled.
  • Then, I tell you how many threes you rolled.
  • Then, I tell you how many fours you rolled.
  • Then, I tell you how many fives you rolled.
  • Then, I tell you how many sixes you rolled.

You're allowed to remember things between each step, typically a "running total" of some sort, ideally without too many distinct possibilities.

Or in terms of Python code with a transition function next_state:

state = None
state = next_state(state, 1, num_ones)
state = next_state(state, 2, num_twos)
state = next_state(state, 3, num_threes)
state = next_state(state, 4, num_fours)
state = next_state(state, 5, num_fives)
state = next_state(state, 6, num_sixes)

For example, if we wanted to find the sum of the dice, we could use the transition function

def next_state(state, outcome, count):
    if state is None:
        return outcome * count
        return state + outcome * count

What about dropping dice? We can augment the state with a "count-list", and whenever we decide that $k$ dice rolled the current outcome, we pop $k$ elements off the list and use them to determine how many of those dice should contribute. Thus, this formulation also covers SkySpiral7's follow-up question where two dice are dropped, or if an arbitrary number of dice are dropped (granted, not in the form of a generating function). In fact, we can count each sorted position multiple times, or even a negative number of times.

The algorithm

The key to the algorithm's efficiency is dynamic programming. To find the probability distribution for a dice pool of $\ell$ outcomes per die and $n$ dice, the algorithm recursively uses memoized solutions for dice pools of $\ell-1$ outcomes per die and $0 \ldots n$ dice. For example, to compute the solution to a pool of 4d6, we use memoized solutions to pools of 0d5, 1d5, 2d5, 3d5, and 4d5.

Specifically, in each call we pop the greatest outcome from the die, and then for $k = 0 \ldots n$:

  • Compute how many ways there are for $k$ out of $n$ dice to roll the current outcome. This is just the binomial coefficient $\binom{n}{k}$. These coefficients can be efficiently computed and memoized using Pascal's triangle.
  • Recursively compute the solution for a pool with the current outcome removed from the die, and $n - k$ dice in the pool.
  • For each state in the recursive distribution, apply the transition function, giving the current outcome and $k$ as the other arguments. Then add the resulting state to the output distribution with weight equal to its recursive weight times the binomial coefficient.

If the last remaining outcome is popped, all $n$ dice must roll that outcome, leading to the base case of a die with an empty set of outcomes and 0 dice in the pool. This is considered to produce a distribution consisting of just the state None (Python's null value) with weight 1.

Here is sample Python code for the basic algorithm:

def solve(die, n):
    outcome, die = die.pop()
    if len(die) == 0:
        state = next_state(None, outcome, n)
        return {state : 1}
    result = defaultdict(int)
    for k in range(n + 1):
        tail = solve(die, n - k)
        for state, weight in tail.items():
            state = next_state(state, outcome, k)
            weight *= comb(n, k)
            result[state] += weight
    return result

Here is an example call graph:

Example call graph.

  • Vertexes are unique calls. Since the algorithm is memoized, the state distribution at each vertex will be computed exactly once.
  • Edges are all calls. Each edge has weight equal to a binomial coefficient.
  • Each path from a vertex to the sink (base case) corresponds one-to-one with a possible sorted sequence of dice rolls of the starting vertex's dice pool. The product of the weights of the edges on the path is the weight of rolling that sorted sequence. This is equivalent to the decomposition of a multinomial coefficient as the product of binomial coefficients.

While not as fast as the FFT, this dynamic programming approach is much more flexible while remaining reasonably efficient over a variety of dice mechanics: the number of unique calls is $O\left(n^2d\right)$, times some hopefully small number of states and time needed to transition each state and accumulate weights. By changing the transition function, we can do things like finding the size and outcome of the largest matching set:

def next_state(state, outcome, count):
    if state is None:
        return count, outcome
        return max(state, (count, outcome))

Or other things such as looking for straights, evaluating RISK-like mechanics, and more.

Further reading

If you would like to know more, you can read my paper on the subject:

    title={Icepool: Efficient Computation of Dice Pool Probabilities},
    author={Albert Julius Liu},
    booktitle={Eighteenth AAAI Conference on Artificial Intelligence and Interactive Digital Entertainment},
    venue={Pomona, California},

Or try out my Icepool implementation, which includes several demo web apps and JupyterLite notebooks. For the particular problem of drop-lowest, you can try my Cortex Prime calculator, which will compute the probability distribution of the sum of a (possibly mixed) pool of standard dice.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.