Dice probability problem

[Dave Mandelin]

That looks reasonable. The operation you are implementing is known as
'convolution' and is equivalent to multiplying polynomials. It would be
a little more general if you had the input 'die' be a sequence of the
count for each outcome, so d6 would be [1]*6 (or [0]+[1]*6 if you
prefer). That would allow allow you to represent unfair dice and also
to add not just a die to a distribution, but to add any two
distributions, so you can play tricks like computing 16d6 as
(d6)*2*2*2*2. (The code above is a convolution that restricts the
second distribution 'die' to have only 0 and 1 coefficients.) The
general convolution can be implemented much like what you have, except
that you need another multiplication (to account for the fact that the
coefficient is not always 0 or 1). My not particularly efficient
implementation:

def vsum(seq1, seq2):
    return [ a + b for a, b in zip(seq1, seq2) ]

def vmul(s, seq):
    return [ s * a for a in seq ]

# Convolve 2 sequences
# equivalent to adding 2 probabililty distributions
def conv(seq1, seq2):
    n = (len(seq1) + len(seq2) -1)
    ans = [ 0 ] * n
    for i, v in enumerate(seq2):
        vec = [ 0 ] * i + vmul(v, seq1) + [ 0 ] * (n - i - len(seq1))
        ans = vsum(ans, vec)
    return ans

# Convolve a sequence n times with itself
# equivalent to multiplying distribution by n
def nconv(n, seq):
    ans = seq
    for i in range(n-1):
        ans = conv(ans, seq)
    return ans

print nconv(3, [ 1 ] * 6)
print nconv(3, [ 1.0/6 ] * 6)
print nconv(2, [ .5, .3, .2 ])

[1, 3, 6, 10, 15, 21, 25, 27, 27, 25, 21, 15, 10, 6, 3, 1]
[0.0046296296296296294, 0.013888888888888888, 0.027777777777777776,
0.046296296296296294, 0.069444444444444448, 0.097222222222222238,
0.11574074074074074, 0.125, 0.125, 0.11574074074074073,
0.097222222222222224, 0.069444444444444448, 0.046296296296296294,
0.027777777777777776, 0.013888888888888888, 0.0046296296296296294]
[0.25, 0.29999999999999999, 0.29000000000000004, 0.12,
0.040000000000000008]