掷骰子的总和

Question

我正在尝试计算从 n 卷 s 面骰子中获得特定总和的概率。我在这个link（公式10）中找到了公式。

这是我用 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);
  }
}

以下是结果：

两个 6 面骰子：

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

对于三个 6 面骰子：

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

负概率？代码或公式有什么问题？

这是 Valgrind 报告

==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

您的 C 和 fact 函数出现溢出。 fact(13) 已经溢出有符号的 32 位整数。

您可以对C 使用此定义，这样就完全不需要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;
}

这避免了大的中间结果，并将结果累积在 double 而不是 int 中。使用双精度可能并不总是令人满意，但您的代码无论如何都会将结果转换为双精度，所以这里看起来不错。

这是我对原始问题的解决方案。它完全避免了计算能力和组合，从而产生更短、更快且数值更稳定的解决方案。你可以运行它，例如 ./dice 3d12 来产生掷出 3 个 12 面骰子的概率。

#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;
}

Answer 2

在细微的变化中，系数计算似乎打算截断整数，并且替换 double 类型可能会产生不同的结果。 uint64_t 为所有合理的s 和n 查询提供足够的存储空间。您可以执行类似于以下的操作，这也消除了对 pow(-1, k) 的任何依赖，将其替换为上面评论中详述的简单位比较。

考虑到这一点，并利用精确类型实现可移植性，您可以执行类似于以下的操作：

#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;
}

使用/输出示例

$ ./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