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
else:
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:
@cache
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:
- 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
else:
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},
volume={18},
number={1},
pages={258-265},
year={2022},
month={Oct.},
eventdate={2022-10-24/2022-10-28},
venue={Pomona, California},
url={https://ojs.aaai.org/index.php/AIIDE/article/view/21971},
doi={10.1609/aiide.v18i1.21971}
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.
