The sum of the dice pieces
I am trying to calculate the probability of getting a specific sum of n s-sided bone results. I found the formula in this link (formula 10).
This is the code 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 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
negative probabilities !!! What's wrong with the code or the formula?
this is the 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)
source to share
Functions C
and fact
overflows. fact(13)
already overflows a signed 32-bit integer.
This definition can be used for C
, which completely eliminates the need for fact
:
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 double rather than int. Using double may not always be satisfactory, but your code converts the result to double anyway, so it looks great.
Here's my solution to the original problem. It completely eliminates computational power and combinations, resulting in a shorter, faster, and more numerically robust solution. You can run it like, for example, ./dice 3d12
to create the likelihood 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;
}
source to share
In a small version, it seems that the calculation of the coefficient assumes integer truncation and type substitution double
can lead to different results. uint64_t
provides sufficient storage for all reasonable requests, s
and n
. You can do something similar to the following, which also removes any dependency on pow(-1, k)
, replacing it with a simple bit comparison detailed in the comment above.
With that in mind, and using exact types for portability, you can do something similar to the following:
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <inttypes.h>
int32_t power (int32_t x, uint32_t y)
{
int32_t temp;
if (y == 0)
return 1;
temp = power (x, y / 2);
if ((y % 2) == 0)
return temp * temp;
else
return x * temp * temp;
}
uint64_t factorial (uint32_t v) {
if (v <= 1)
return 1;
uint64_t n = 1;
for (uint32_t i = 1; i <= v; i++)
n *= i;
return n;
}
int32_t coeff (uint32_t p, uint32_t n, uint32_t s)
{
uint32_t c = 0, pn_s = (p - n)/s;
for (uint32_t k = 0; k <= pn_s; k++)
{
int32_t neg1k = ((k & 1) ? -1 : 1),
vp = (factorial (n)/(factorial (k) * factorial (n - k))) *
(power (p - s * k - 1, n - 1))/(n - 1);
c += neg1k * vp;
}
return c;
}
double ppns (uint32_t p, uint32_t n, uint32_t s)
{
uint32_t s2n = 1;
double invs2n;
for (uint32_t i = 0; i < n; i++)
s2n *= s;
invs2n = (double)1.0 / s2n;
return invs2n * coeff (p, n, s);
}
int main (int argc, char **argv) {
uint32_t n, p, s;
n = argc > 1 ? (uint32_t)strtoul (argv[1], NULL, 10) : 2;
s = argc > 2 ? (uint32_t)strtoul (argv[2], NULL, 10) : 6;
for (p = n; p <= s * n; p++) {
double prob = ppns (p, n, s);
printf (" points: %3" PRIu32 ", probablility: %8.4lf\n", p, prob);
}
return 0;
}
Usage / output example
$ ./bin/dice_prob
points: 2, probablility: 0.0278
points: 3, probablility: 0.0556
points: 4, probablility: 0.0833
points: 5, probablility: 0.1111
points: 6, probablility: 0.1389
points: 7, probablility: 0.1667
points: 8, probablility: 0.1389
points: 9, probablility: 0.1111
points: 10, probablility: 0.0833
points: 11, probablility: 0.0556
points: 12, probablility: 0.0278
$ ./bin/dice_prob 3
points: 3, probablility: 0.0093
points: 4, probablility: 0.0185
points: 5, probablility: 0.0370
points: 6, probablility: 0.0556
points: 7, probablility: 0.0833
points: 8, probablility: 0.1111
points: 9, probablility: 0.1204
points: 10, probablility: 0.1250
points: 11, probablility: 0.1204
points: 12, probablility: 0.1065
points: 13, probablility: 0.0833
points: 14, probablility: 0.0509
points: 15, probablility: 0.0370
points: 16, probablility: 0.0185
points: 17, probablility: 0.0093
points: 18, probablility: 0.0000
source to share