Efficient computation of homogeneous roll-and-keep dice pool probabilities

“Roll-and-keep” is a concept popularized by the Legend of the Five Rings RPG. From 1st to 4th Edition, L5R used the following system: the player rolls a pool of (exploding) d10s, keeping a number of the highest dice rolled. This is thus a generalization of keep-highest dice pools (= keep one) and success-counting dice pools (= keep all).

A similar concept can be applied to many schemes that have been proposed for D&D-style character generation. For example, roll 5d6, keeping the highest three dice; and/or rolling seven ability scores, keeping the highest six.

However, computing the probabilities for this type of dice pool is not as easy as it is for keep-highest or success-counting. In this article, we’ll look at systems where the player rolls a homogeneous pool of n dice and keeps a number k of the highest dice. Each die’s outcomes have a finite range r. Without loss of generality, we’ll say that the lowest roll on a die is 0.

Existing calculators

Popular existing dice probability calculators are not capable of efficiently computing the probabilities for this type of system in closed form. By “closed form” I mean “not Monte Carlo” ; while Monte Carlo methods are very flexible, they can take a long time to converge, especially compared to the precision of a direct calculation.

Here are some examples:

AnyDice

AnyDice appears to take exponential time in n. For example,

output [highest 5 of 10d[explode d10]]

times out. (n.b. At default, AnyDice explodes dice up to two times, not counting the base roll. In this example, the maximum outcome on each die is r = 30 or triple that of a single die.) Neither is exploding the only issue, though it exacerbates it;

output [highest 8 of 16d10]

times out as well, as does

output [highest 2 of 20d10]

Troll

The conference paper for the Troll dice roller and probability calculator says:

For example, the game “Legend of Five Rings” [18] uses a roll mechanism that can be described in Troll as

sum (largest M N#(sum accumulate x := d10 until x<10))

where M and N depend on the situation. With M = 3, N = 5 and the maximum number of iterations for accumulate set to 5, Troll takes nearly 500 seconds to calculate the result, and it gets much worse if any of the numbers increase.

To be fair, this paper was presented more than a decade ago, and hardware has advanced considerably since then. cpubenchmark.net shows a factor of about 5× in single-threaded performance and 33× in overall performance from the CPU used in the paper to my current desktop CPU. However, 500 seconds is a long time, and even this advancement in hardware is not nearly enough to overcome even modest increases in the number of dice due to the exponential time of the algorithm.

Troll does do better than AnyDice on the last of the three.

SnakeEyes

A Lua-language calculator that runs locally in the browser. (Hence why I don’t feel bad linking a futile L5R script in this case, as it’s not putting load on someone else’s server. I don’t recommend clicking on that if your computer is prone to hanging though.) Unfortunately, this too cannot complete the first two calculations in a reasonable amount of time, though it also does better than AnyDice on the third. The author also did come up with an optimized script for L5R in particular.

I believe the difference is that AnyDice enumerates all r^n possible rolls when doing roll-and-keep, whereas SnakeEyes and possibly Troll only keep track of the k highest dice, which has only O \left( r^k \right) possibilities (times some polynomial factor when looking at the net time). More on this later.

lynks.se

There does exist a website with L5R probabilities: lynks.se. However, if you look at that webpage’s source code, you’ll see that the probabilities have been precomputed and stored in static arrays. Furthermore, there is a small but noticeable amount of noise in the data, indicating that it was computed via Monte Carlo rather than a closed-form method.

pyPBE

Another Monte Carlo application, this time for D&D-style character generation. Web version.

I don’t mean this as harsh criticism of any of these resources. They are very accessible and are a great service to the RPG community. It also may not be simple to integrate semi-specialized algorithms into a generalist online dice calculator.

But, like the archetypical adventurer, we always want more.

Convolution

How can we compute these probabilities more efficiently? Let’s look at a simpler problem first: finding the probability mass function (PMF; “Normal” data on AnyDice) of a sum of dice without dropping any of them.

The brute-force method is to enumerate all r^n possible rolls for the entire pool and sum the dice for each, which would take O \left( nr^n \right) time, which is exponential.

A much faster method to to do the addition via iterated convolutions. Instead of the brute force method, we could iteratively take the PMF of the first i dice, which can take on a range O \left( nr \right) of possible outcomes, and convolve it with the PMF of a single additional die to get the PMF of the first i+1 dice. I would guess that most of the popular dice calculators do this, whether their creators conceptualized it as a convolution or not. The time to compute the PMF of the sum of n dice is then the sum of the time taken by the n convolutions, each of which takes O \left( nr \cdot r \right) time, for a net time of

O \left(n^2 r^2 \right)

which is nicely polynomial compared to the exponential time of the brute-force algorithm.1

Unfortunately, the roll-and-keep system makes things more difficult. If we keep only the PMF of the sum from iteration to iteration, we have no way of correctly dropping the new lowest die as we consider each die in turn. On the other hand, if we track the probability distribution of all k dice kept so far, this requires \Theta \left( r^k \right) possible sets of kept dice to be processed during each of \Theta \left( n \right) iterations. Without even counting the time it takes to process each set, this is \Omega \left( nr^k\right) time, which is exponential again. (Though this is still better than the brute-force method, since if k is constant this is polynomial.) It’s possible to do better by keeping the kept dice in sorted order, though this is still exponential.

Cutting up the probability space

Instead, what we’re going to do is cut up the probability space into a polynomial number of disjoint pieces, where the probability of landing in each is the product of a number of factors which we can compute using iterated convolutions just like the sum of dice above. We’re going to parameterize the pieces using the following parameters:

  • The lowest single roll \ell among the kept dice.
  • The number of kept dice s \leq k - 1 that rolled strictly greater than \ell.
  • The sum S of those s dice.

The total result of each piece is then S + \ell \cdot \left(k - s\right).

The probability of landing in one of these pieces is equal to the product of the following:

  • The lo factor: The probability that, out of n - s dice that were less than or equal to the lowest kept die, all of them will roll \ell or less, and at least k - s of them roll exactly equal to \ell (in order to be consistent with \ell being in fact the lowest kept die).
  • The sum factor: The probability that, out of s dice, all of them will roll strictly greater than \ell and have sum exactly S.
  • The binomial factor: The number of ways of assigning n dice to the above two categories.

The lo factor

To compute the chance that at least i_= = k - s dice will roll exactly equal to \ell, we can compute the chance that exactly i_= dice will roll exactly equal to \ell and then take a cumulative sum in reverse order. We’ll compute a full array of these values. Given indexes i_\leq, \ell, i_=, the array will tell us the chance that i_\leq dice will all roll less than or equal to \ell, with i_= of those dice rolling exactly equal to \ell.

Consider each vector corresponding to i_\leq, \ell . The elements of the vector are indexed by i_=. The vector for 1, \ell just consists of two elements \left[ P_<\left( \ell \right), P_=\left( \ell \right) \right], where the elements are the chance of a single die rolling strictly less than \ell and exactly equal to \ell respectively. To compute the vectors for i_\leq+1, \ell, we can just convolve the vector for i_\leq, \ell with the vector for a single die 1, \ell.

i_= can range up to n, and \ell has r possible values. Computing this array for all possible values of i_\leq, \ell, i_= therefore requires O \left( nr \right) convolutions and cumulative sums, each of which takes O \left(n \right) time, for a net time of O \left( n^2 r \right) for this step.

Here’s a pictorial representation for \ell = 3 on d6s:

⚂ × 0⚂ × 1 ⚂ × 2 ⚂ × 3
1 die⚀, ⚁ = 2/6⚂ = 1/6
2 dice⚀⚀, ⚀⚁, ⚁⚀, ⚁⚁ = 4/36⚀⚂, ⚁⚂, ⚂⚀, ⚂⚁ = 4/36⚂⚂ = 1/36
3 dice⚀⚀⚀, ⚀⚀⚁, ⚀⚁⚀, ⚀⚁⚁, ⚁⚀⚀, ⚁⚀⚁, ⚁⚁⚀, ⚁⚁⚁ = 8/216⚀⚀⚂, ⚀⚁⚂, ⚀⚂⚀, ⚀⚂⚁, ⚁⚀⚂, ⚁⚁⚂, ⚁⚂⚀, ⚁⚂⚁, ⚂⚀⚀, ⚂⚀⚁, ⚂⚁⚀, ⚂⚁⚁ = 12/216⚀⚂⚂, ⚁⚂⚂, ⚂⚀⚂, ⚂⚁⚂, ⚂⚂⚀, ⚂⚂⚁ = 6/216 ⚂⚂⚂ = 1/216

Each row/vector is the convolution of the previous vector and the first vector. Note how instead of counting all the individual possibilities that fall into each cell, we can convolve the previous vector and the first vector.

The sum factor

Again, we’ll compute an array, this time indexed by s, \ell, S. The vector for s, \ell with s = 1 is just the PMF for a single die with all outcomes with index less than or equal to \ell removed. Just like the lo factor, we can compute the vector for s+1, \ell by convolving the vector s, \ell with the vector 1, \ell.

s can be up to k-1, and \ell once again has r possible values. Therefore, this requires O \left( kr \right) convolutions, each of which takes time O \left( kr \cdot r \right), for a net time of O \left( k^2r^3 \right).

The binomial factor

This is just the binomial coefficient {n \choose s}, which can be computed for all s \leq k - 1 in O \left( n^2 \right) time, which is less than that of the the lo factor.

The outer loop

Now we just have to loop, multiply, and accumulate the three factors over all values of s, \ell, S, which is O \left( k \cdot r \cdot nr \right) = O \left( knr^2 \right) iterations. The net time is thus

O \left( n^2 r + knr^2 + k^2 r^3 \right) = O \left( r \left( n + kr \right)^2 \right)

On my desktop (Intel i5-8400), my NumPy-based implementation can compute the PMF of any of the examples above in under 5 milliseconds.

Similar strategies can be used to compute the PMF of:

  • The sum of the lowest k dice.
  • The kth lowest or highest single die.
  • The sum of the k_ath to the k_bth lowest dice, i.e. dropping dice from both the top and the bottom, with the binomial factor becoming a multinomial factor. Though this is getting sketchy from both a computational and from a game design perspective.

Source code

GitHub. At time of writing the API is far from stable, so I wouldn’t advise building anything serious around this just yet.

1 It’s possible to do even better asymptotically using a fast Fourier transform. However, according to scipy this is estimated to only be faster than the quadratic direct algorithm past about a range of 500, which is greater than pretty much any tabletop RPG will make use of. There’s also the option of doing binary splits but a) this only improves the time by a constant factor, and b) our final algorithm needs the intermediate PMFs anyways.

It is conceivable will lose less to floating-point precision, but I haven’t worked this out yet.

Leave a comment