当相对于随机数进行记录时，避免在C中除以零

Question

我目前正在使用C来生成高斯噪声。在一个步骤中，我需要记录均匀分布的数字u1 = (double) rand() / RAND_MAX。因为u1可能为零，所以有做log(u1)的风险。所以，我需要检查一下。我应该用吗？

do {
    u1 = ((double) rand() / RAND_MAX);
} while (u1 == 0.);

或者，我应该使用

do {
    u1 = ((double) rand() / RAND_MAX);
} while (u1 < epsilon);

哪里epsilon是一个小数字？如果后者是首选，我该如何选择epsilon的值？ （在Fortran中有TINY，但我不知道在C中做什么）。

附件是完整的代码：

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

double gaussian_noise(double mean, double std)
{
    static int have_spare = 0;
    static double u1, u2, z1, z2;
    if(have_spare)
    {
        have_spare = 0;
        z2 = sqrt(-2. * log(u1)) * sin(2. * M_PI * u2);
        return mean + std * z2;
    }
    have_spare = 1;
    do {
    u1 = ((double) rand() / RAND_MAX);
    } while (u1 == 0.);
    u2 = ((double) rand() / RAND_MAX);
    z1 = sqrt(-2. * log(u1)) * cos(2. * M_PI * u2);
    return mean + std * z1;
}

void main()
{
    const double mean = 0., std = 1.;
    double noise;
    int i;
    for(i=0; i<100000; i++)
    {
        noise = gaussian_noise(mean, std);
        printf("%lf\t", noise);
    }
}

Answer 1

你只需要确保rand()的结果不是0，这样你就不必一次又一次地转换为double和division

int r = rand();
while (r == 0)
    r = rand();

u1 = (double) rand() / RAND_MAX;

更简单的解决方案是

u1 = (double)(rand() + 1U)/(RAND_MAX + 1U);

这样你就不再需要循环了

但是，您应该为双尾数生成53位以获得更合适的精度，而不是（通常）使用rand()仅15或31位。一个示例解决方案是Java Random().nextDouble()

public double nextDouble() {
    return (((long)next(26) << 27) + next(27))
         / (double)(1L << 53);
}

Answer 2

正如nhahtdh所说，零是唯一的问题数字，因此它是唯一被抛出的问题。将浮点数与==和!=进行比较通常不是优选的。然而，在这种情况下，它正是你想要的：任何不完全是浮点(double)(0.)的东西。

Answer 3

如@nhahtdh所述，OP的当前代码就足够了。

但我怀疑是“一分错误”。

OP的代码将生成范围：0.0（不包括）到1.0（包括）。然而，对于这个任务，我期望0.0（包括）到1.0（不包括），0.0被删除。

代码是否希望范围0.0（不包括）为1.0（独占）考虑

u1 = ((double) rand() + 0.5) / (RAND_MAX + 1.0);