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

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