How I generate the distribution for 3d6 (and others)

Introduction

How does one generate the distribution for rolling s-sided dice r times? If you think about it, this is equivalent to finding the coefficients of the expanded expression:

(1+x+x²+ ... +x^s)^r (1)

Straightforward, eh? Well, anyway, this next section demonstrates an efficient algorithm for computing these coefficients. (I'm open to comments. Anybody have a better way?)

Note that the obvious solution to this, counting in base 's', is extremely inefficient (i.e. - O(S^N)), and as a result, completely unreasonable.

How to calculate those damned coefficients

Some identities demonstrate that

(1) == [(1-x^s+1) / (1-x)]^r (2)
(2) == (1-x^s+1)^r * [1 / (1-x)]^r (3)

Now say that

A == (1-x^s+1)^r (4)
B == [1 / (1-x)]^r (5)

(So that (2) == A * B)

In the next section, we define C(n, r) == n! / (r! * (n-r)!). This is the standard 'n choose r' thing we keep hearing about.

Some more identities tell us the following:

A == 1 - C(r, 1)*x^s+1 + C(r, 2)*x^2(s+1) + ... + (-1)^k * C(r, k) * x^k(s+1) + ... + (-1)^r * C(r, r) * x^r(s+1)
B == 1 + C(1+r-1, 1)*x + ... + C(k+r-1, k)*x^k + ...

But now we know that (2) == A * B, and more to the point, (1) == A * B! This is what we need. Finding the coefficients in A * B turns out to be really easy, and cheap.

To find the coefficient of x^v in A * B, just choose from A all those coefficients for powers of xⁱ in A such that i < v. Then, for each of these coefficients, choose the coefficient of x^j in B such that i+j == v. Then just multiply the two coefficients together. Sum all of these products, and you have the coefficient of x^v!

I believe that this is O(N²), although I could be convinced otherwise.

Comments? Criticisms? Confusion? Mail me, and I'll try to help.