Calculate probability of a fair dice roll (in non-exponential time)

Question

Variations on this are pretty common questions, but all my google-fu has left me stumped. I would like to calculate the odds of a fair dice roll, but I want to do it efficiently. There are lots of examples out there of how to do this, but all the algorithms I've found are too computationally expensive (exponential time) to work for large numbers of dice with many sides.

The Simple Problem: Calculate the odds of a roll of n, on x y sided dice.

The Simple Solution: Create the n-ary Cartesian product of the roll, sum each product, count the number of times the sum is the target, do a little division and voila.

Example of a Simple Solution in Go: https://play.golang.org/p/KNUS4YBQC0g

The Simple Solution works perfectly. I extended it to allow for cases like dropping the highest/lowest n faces, and the results hold up to spot testing.

But consider {Count: 20,Sides: 20,DropHighest: 0,DropLowest:0, Target: 200}.

If I evaluated that with the previous solution, my "table" would have 104 some odd septillion cells and will max the CPU pretty easily.

Is there a more efficient way to calculate probability for large numbers of dice with many sides? If so, can it account for more complex selection of "success" conditions like dropping some dice?

I'm convinced it's possible due to the existence of this beautiful website: https://anydice.com/program/969

EDIT:

The solution that worked best for me was David Eisenstat's answer, which I ported to go: https://play.golang.org/p/cpD51opQf5h

Answer 1

Here's some code that handles dropping low and high rolls. Sorry for switching to Python, but I needed easy bignums and a memoization library to keep my sanity. I think the complexity is something like O(count^3 sides^2 drop_highest).

The way this code works is to divide the space of possibilities for rolling count dice each with sides sides by how many dice are showing the maximum number (count_showing_max). There are binomial(count, count_showing_max) ways to achieve such a roll on uniquely labeled dice, hence the multiplier. Given count_showing_max, we can figure out how many maxed dice get dropped for being high and how many get dropped for being low (it happens) and add this sum to the outcomes for the remaining dice.

#!/usr/bin/env python3
import collections
import functools
import math


@functools.lru_cache(maxsize=None)
def binomial(n, k):
    return math.factorial(n) // (math.factorial(k) * math.factorial(n - k))


@functools.lru_cache(maxsize=None)
def outcomes(count, sides, drop_highest, drop_lowest):
    d = collections.Counter()
    if count == 0:
        d[0] = 1
    elif sides == 0:
        pass
    else:
        for count_showing_max in range(count + 1):  # 0..count
            d1 = outcomes(count - count_showing_max, sides - 1,
                          max(drop_highest - count_showing_max, 0),
                          drop_lowest)
            count_showing_max_not_dropped = max(
                min(count_showing_max - drop_highest,
                    count - drop_highest - drop_lowest), 0)
            sum_showing_max = count_showing_max_not_dropped * sides
            multiplier = binomial(count, count_showing_max)
            for k, v in d1.items():
                d[sum_showing_max + k] += multiplier * v
    return d


def main(*args):
    d = outcomes(*args)
    denominator = sum(d.values()) / 100
    for k, v in sorted(d.items()):
        print(k, v / denominator)


if __name__ == '__main__':
    main(5, 6, 2, 2)

Answer 2

You can compute the distribution of sums of x y-sided dice by multiplying out the following polynomial in variable Z:

(Z + Z^2 + ... + Z^x)^y / x^y.

For example, for two six-sided dice:

(Z + Z^2 + ... + Z^6)^2 / 6^2
  = (Z + Z^2 + ... + Z^6) * (Z + Z^2 + ... + Z^6) / 36
  = (Z^2 + 2Z^3 + 3Z^4 + 4Z^5 + 5Z^6 + 6Z^7 + 5Z^8 + 4Z^9 + 3Z^10 + 2Z^11 + Z^12) / 36,

so you can read off the probability of getting sum 6 as the coefficient of Z^6 (5/36).

For three "two-sided" dice:

(Z + Z^2)^3 / 2^3 = (Z + Z^2) * (Z + Z^2) * (Z + Z^2) / 8
                  = (Z^2 + 2Z^3 + Z^4) (Z + Z^2) / 8
                  = (Z^3 + 3Z^4 + 3Z^5 + Z^6) / 8,

so the probability of getting sum 4 is the coefficient of Z^4 (3/8).

You can use the school algorithm to get a polynomial algorithm for this problem. Lightly tested Go code:

package main

import "fmt"

func dieRolls(x, y int) map[int]float64 {
    d := map[int]float64{0: 1.0}
    for i := 0; i < x; i++ {
        d1 := make(map[int]float64)
        for j := 1; j <= y; j++ {
            for k, v := range d {
                d1[k+j] += v / float64(y)
            }
        }
        d = d1
    }
    return d
}

func main() {
    for k, v := range dieRolls(2, 6) {
        fmt.Printf("%d %g\n", k, v)
    }
    fmt.Printf("\n")
    for k, v := range dieRolls(3, 2) {
        fmt.Printf("%d %g\n", k, v)
    }
}

Run the code: https://play.golang.org/p/O9fsWy6RZKL