Multiset enumeration

For the most part, AnyDice operates by enumerating all possible multisets that could come out of a pool of dice. For \$n\$ dice of \$y\$ sides each, the number of possible multisets is

$$\binom{n + y - 1}{n}$$

An example of how multiset enumeration could be implemented is given in the documentation for Python's itertools.combinations_with_replacement():

def combinations_with_replacement(iterable, r):
    # combinations_with_replacement('ABC', 2) --> AA AB AC BB BC CC
    pool = tuple(iterable)
    n = len(pool)
    if not n and r:
        return
    indices = [0] * r
    yield tuple(pool[i] for i in indices)
    while True:
        for i in reversed(range(r)):
            if indices[i] != n - 1:
                break
        else:
            return
        indices[i:] = [indices[i] + 1] * (r - i)
        yield tuple(pool[i] for i in indices)

Each multiset should then be weighted according to the multinomial distribution.

Multiset enumeration is easy to use since we can see the entire sequence of sorted rolls, and it's fast enough for most cases on the tabletop. However, it can struggle with large numbers of many-sided or exploding dice such as classic Legend of the Five Rings. The 5-second limit tops out at around eight to nine d20s, or between 4.5 and 7 million possibilities. Opposed pools are also challenging since we need to consider the product of possibilities of each pool; for a RISK-like mechanic, AnyDice times out on 4d6 vs. 4d6.

Order statistics

Wikipedia has an explicit formula (permalink) for extracting a single ranked die from a pool, which computes the same result as AnyDice's @ operator with a single index. This is formally known as an order statistic. There are also formulas for joint order statistics (i.e. selecting multiple indexes at a time), but these get complicated quickly -- and once more than a few dice are kept, not particularly computationally efficient either.

Fast Fourier transform

That the sum of a set of dice can be computed via convolution is a classical result in discrete statistics. Convolution (sometimes also framed in terms of generating functions) allows us to use the fantastically efficient fast Fourier transform (FFT). With some more work, FFTs can also be used the case where only some of the dice are kept. Approaches of this type have been proposed on StackExchange, with examples including:

• Formula for dropping dice (non-brute force)
• Generating Function for sum of N dice [or other multinomial distribution] where lowest N values are "dropped" or removed

Dynamic programming

However, my preferred approach is dynamic programming (sometimes also framed in terms of recurrence relations). Examples include:

• Roll and Keep in Anydice
• How to calculate the probabilities for eliminative dice pools (dice cancelling mechanic) in Neon City Overdrive?

The short of it is that we compute the answer to the problem to "roll \$n\$ d20, keep \$x\$" by considering how many dice \$k = 0 \ldots n\$ might possibly roll a 20, and then use the solutions to "roll \$n - k\$ d19, keep \$x - k\$" -- namely, the remaining pool after the \$k\$ dice that rolled a 20 are removed -- to complete the solution. In turn, the solutions for d19s use the solutions for d18s, and so forth until we reach d1s and all of the remaining dice must roll 1. While not as fast as FFTs, the dynamic programming approach is much more flexible, covering general "single-pass" functions over the order statistics: in addition to summing the dice, it can also be used to produce efficient solutions for matching sets, straights, RISK-like mechanics, and so forth.

If you are interested in learning more, you can check out my Icepool Python package or read my paper on the subject.