Avoid division by zero in C when taking log with respect to a random number

Question

I am currently using C to generate a Gaussian noise. In one step, I need to take the log of a uniformly distributed number u1 = (double) rand() / RAND_MAX. Because u1 may be zero, there is a risk in doing log(u1). So, I need to check. Should I use

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

Or, should I use

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

where epsilon is a small number? If the latter is preferred, how should I choose the value of epsilon? (In Fortran there is TINY, but I do not know what to do in C).

Attached is the complete code:

#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

As nhahtdh says, zero is the only problem number, so it's the only one to throw out. Comparing floating point numbers with == and != is usually not preferable. In this case, however, it's exactly what you want: anything that isn't exactly the floating point (double)(0.) passes.

Answer 2

You just need to make sure that the result of rand() isn't 0 so that you won't have to do the conversion to double and division again and again

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

u1 = (double) rand() / RAND_MAX;

An even simpler solution is

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

This way you don't need the loop anymore

However you should generate 53 bits for the double mantissa to get more proper precision and distribution instead of (typically) only 15 or 31 bits with rand(). A sample solution is Java Random().nextDouble()

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

[In early versions of Java, the result was incorrectly calculated as:

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

This might seem to be equivalent, if not better, but in fact it introduced a large nonuniformity because of the bias in the rounding of floating-point numbers: it was three times as likely that the low-order bit of the significand would be 0 than that it would be 1! This nonuniformity probably doesn't matter much in practice, but we strive for perfection.]