Sum of dice rolls

Question

I am trying to compute the probability of getting a specific sum from n rolls of s-sided dice. I found the formula in this link (formula 10).

This is the code that I wrote in C:

# include <stdio.h>
# include <stdlib.h>
# include <math.h>

# define n      2   // number of dices
# define s      6   // number of sides of one dice

int fact(int x){
  int y = 1;
  if(x){
    for(int i = 1; i <= x; i++)
      y *= i;
  }
  return y;
}
int C(int x,int y){
  int z = fact(x)/(fact(y)*fact(x-y));
  return z;
}

int main(){
  int     p,k,kmax;
  double  proba;

  for(p = n; p <= s*n; p++){
    proba = 0.0;
    kmax = (p-n)/s;
    for(k = 0; k <= kmax; k++)
      proba += pow(-1.0,k)*C(n,k)*C(p-s*k-1,p-s*k-n);
    proba /= pow((float)s,n);
    printf("%5d %e\n",p,proba);
  }
}

Here are the results for:

Two 6-sided dice:

2    2.777778e-02
3    5.555556e-02
4    8.333333e-02
5    1.111111e-01
6    1.388889e-01
7    1.666667e-01
8    1.388889e-01
9    1.111111e-01
10   8.333333e-02
11   5.555556e-02
12   2.777778e-02

and for three 6-sided dice:

3    4.629630e-03
4    1.388889e-02
5    2.777778e-02
6    4.629630e-02
7    6.944444e-02
8    9.722222e-02
9    1.157407e-01
10   1.250000e-01
11   1.250000e-01
12   1.157407e-01
13   9.722222e-02
14  -1.805556e-01
15  -3.703704e-01
16  -4.768519e-01
17  -5.462963e-01
18  -6.203704e-01

Negatives probabilities? What's wrong in the code or in the formula?

This is Valgrind report

==9004== 
==9004== HEAP SUMMARY:
==9004==     in use at exit: 0 bytes in 0 blocks
==9004==   total heap usage: 1 allocs, 1 frees, 1,024 bytes allocated
==9004== 
==9004== All heap blocks were freed -- no leaks are possible
==9004== 
==9004== For counts of detected and suppressed errors, rerun with: -v
==9004== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 0 from 0)

Answer 1

Your C and fact functions experience overflow. fact(13) already overflows a signed 32 bit integer.

You can use this definition for C, which avoids the need for fact entirely:

double C(int n, int k) {
    double r = 1.0;
    for (int i = 0; i < k; i++) {
        r *= n - i;
        r /= i + 1;
    }
    return r;
}

This avoids large intermediate results, and accumulates the result in a double rather than an int. Using a double may not always be satisfactory, but your code is converting the result to a double anyhow, so it seems fine here.

Here's my solution to the original problem. It avoids computing powers and combinations entirely, resulting in a shorter, faster and more numerically stable solution. You can run it, for example like ./dice 3d12 to produce the probabilities of rolling 3 12-sided dice.

#include <stdio.h>
#include <stdlib.h>

void dice_sum_probabilities(int n, int d) {
    // Compute the polynomial [(x+x^2+...+x^d)/d]^n.
    // The coefficient of x^k in the result is the
    // probability of n dice with d sides summing to k.
    int N=d*n+1;
    // p is the coefficients of a degree N-1 polynomial.
    double p[N];
    for (int i=0; i<N; i++) p[i]=i==0;
    // After k iterations of the main loop, p represents
    // the polynomial [(x+x^2+...+x^d)/d]^k
    for (int i=0; i<n; i++) {
        // S is the rolling sum of the last d coefficients.
        double S = 0;
        // This loop iterates backwards through the array
        // setting p[j+d] to (p[j]+p[j+1]+...+p[j+d-1])/d.
        // To get the ends right, j ranges over a slightly
        // larger range than the array, and care is taken to
        // not write out-of-bounds, and to treat out-of-bounds
        // reads as 0.
        for (int j=N-1;j>=-d;j--) {
            if (j>=0) S += p[j];
            if (j+d < N) {
                S -= p[j+d];
                p[j+d] = S/d;
            }
        }
    }
    for (int i=n; i<N; i++) {
        printf("% 4d: %.08lf\n", i, p[i]);
    }
}

int main(int argc, char **argv) {
    int ndice, sides, ok;
    ok = argc==2 && sscanf(argv[1], "%dd%d", &ndice, &sides)==2;
    if (!ok || ndice<1 || ndice>1000 || sides<1 || sides>100) {
        fprintf(stderr, "Usage: %s <n>d<sides>. For example: 3d6\n", argv[0]);
        return 1;
    }
    dice_sum_probabilities(ndice, sides);
    return 0;
}