It is a very old exercise to compute the probability to get a given sum s in a given number m of dice throws. E.g., for two throws the combinations, which yield the sum of 7 are

The total number of combinations is 1/36. This yields chance of 1/6 to get the sum of 7.

We can simulate that in Euler. For this, we take a matrix of nx2 dice throws, sum the rows, make that a row vector, and count the multiplicities of 2 to 12.

>n=1000000; D=sum(randint(n,2,6))'; c=getmultiplicities(2:12,D); c/n

[0.02805, 0.055958, 0.083095, 0.111181, 0.139039, 0.166, 0.139126, 0.111384, 0.083124, 0.055471, 0.027572]

The 6th value belongs to the sum 7, and it should be close to 1/6.

>(c/n)[6], 1/6

0.166 0.166666666667

We get the following experimental fractions for each sum from our sample.

>columnsplot(c/n,lab=2:12):

Our first approach to get the exact frequencies of sums in m dice throws involves a function, which creates all possible dice throws and calls another function f$ for each. This is the same approach as in the permutation utilities of Euler. See

Our function is recursive.

>function forAllRek (f$,n,x,k) ... ## create all combinations of 1 to n in ## x[1:k]. For each combination, call f$(x). if k>1 then loop 1 to n x[k]=#; forAllRek(f$,n,x,k-1); end else loop 1 to n x[1]=#; f$(x); end; endif; endfunction

We call this starting with k=m.

>function forAll (f$,n,m) ... x=zeros(1,m); forAllRek(f$,n,x,m); endfunction

Let us test by simply printing the combination.

>function pr (x) ... x, endfunction

>forAll("pr",6,2)

[1, 1] [2, 1] [3, 1] [4, 1] [5, 1] [6, 1] [1, 2] [2, 2] [3, 2] [4, 2] [5, 2] [6, 2] [1, 3] [2, 3] [3, 3] [4, 3] [5, 3] [6, 3] [1, 4] [2, 4] [3, 4] [4, 4] [5, 4] [6, 4] [1, 5] [2, 5] [3, 5] [4, 5] [5, 5] [6, 5] [1, 6] [2, 6] [3, 6] [4, 6] [5, 6] [6, 6]

That worked.

Now let us increase a vector v[s] by 1 for each combination x, which has the sum s. We use a global vector v$$ for this.

For simplicity, we take the dice values as 0 to n-1. So we use x-1 instead of x.

>function sumup (x) ... s=sum(x-1)+1; v$$[s]=v$$[s]+1; endfunction

The following function initializes v$$ and does the job.

>function getsums (n,m) ... v$$=zeros(1,m*(n-1)+1); forAll("sumup",n,m); return v$$; endfunction

Here are the number of times, each sum from 0 to 10 appeared. We know that the 6th values must be 6.

>getsums(6,2)

[1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1]

We can plot this, carefully labelling from 2 to 12.

>k=2; columnsplot(getsums(6,k),lab=k:6*k):

Here is the result for three dice throws.

>k=3; columnsplot(getsums(6,k),lab=k:6*k):

For 6 dice throws, we can get numbers from 6 to 36.

>k=6; columnsplot(getsums(6,k),lab=(0:5*k)+k,width=0.1):

The result looks normally distributed. In fact, the central limit theorem tells us, that a normalized sum of equal random variables converges to the 0-1-normal Gauss distribution.

Here, we have the following numbers.

>z=1:6

[1, 2, 3, 4, 5, 6]

The mean value is 3.5.

>sum(z)/6

3.5

The variance of one dice throw is 35/12.

>fraction sum((z-3.5)^2)/6

35/12

Now, we can plot the distribution and add the correct normal distribution with mean k*3.5, and variance k*35/12.

>k=6; v=getsums(6,k); ... plot2d((k:6*k)-0.5,v/sum(v),>bar); ... plot2d("qnormal(x,k*3.5,sqrt(k*35/12))",>add,color=red):

Of course, the recursive function is not the most speedy way to get the frequencies. However, it requires a bit of thinking to get a faster algorithm.

The idea is the following.

Assume, we know the frequencies for m-1 dice throws. Then a sum of s in m dice throws can occur, if the previous m-1 dice throws had sums s-5 to s, and correspondingly the current dice throw was 0 to 5. So

which is a nice iterative algorithm. Euler lets us take the sum of 6 adjacent values in a vector by folding with [1,1,1,1,1,1]. We need to extend the vector with 0 to the left and right before the folding, which is correct, since those numbers are impossible to get in m-1 dice throws.

>function getv (n,m=6) ... v=ones(1,m); for k=2 to n; v=fold(zeros(1,m-1)|v|zeros(1,m-1),ones(1,m)); end; return v; endfunction

The result is the same, but many times faster.

>getv(6)

[1, 6, 21, 56, 126, 252, 456, 756, 1161, 1666, 2247, 2856, 3431, 3906, 4221, 4332, 4221, 3906, 3431, 2856, 2247, 1666, 1161, 756, 456, 252, 126, 56, 21, 6, 1]