对随机数取对数时，避免在 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

正如 nhahtdh 所说，零是唯一的问题编号，所以它是唯一可以丢弃的编号。将浮点数与== 和!= 进行比较通常是不可取的。然而，在这种情况下，这正是您想要的：任何不完全是浮点 (double)(0.) 的东西都会通过。

Answer 2

您只需要确保rand() 的结果不为0，这样您就不必一次又一次地进行双精度和除法的转换

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

[在Java的早期版本中，结果被错误地计算为：

return (((long)next(27) << 27) + next(27))
    / (double)(1L << 54);

这似乎是等效的，如果不是更好的话，但实际上它引入了很大的不均匀性，因为浮点数的舍入偏差：它是有效数字的低位的三倍将是 0，而不是 1！这种不均匀性在实践中可能无关紧要，但我们力求完美。]

Answer 3

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

但我怀疑存在“关闭 1 个错误”。

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

代码是否应该考虑 0.0（不包含）到 1.0（不包含）的范围

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