I am trying to compute the probability of getting a specific sum of n s-sided dice outcomes. 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);
}
}
and here are the results for : two 6-sided dices :
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 dices :
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)
Your
C
andfact
functions experience overflow.fact(13)
already overflows a signed 32 bit integer.You can use this definition for
C
, which avoids the need forfact
entirely: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.In a slight variation, it appears the coefficient calculation intends integer truncation and substituting a
double
type may produce varied results.uint64_t
provides sufficient storage for all reasonables
andn
inquiries. You could do something similar to the following, which also eliminates any dependency onpow(-1, k)
replacing the same with a simple bit-comparison detailed in the comment above.With that in mind, and making use of exact types for portability, you could do something similar to the following:
Example Use/Output