Avoid division by zero in C when registering the logarithm with respect to a random number
I am currently using C to generate gaussian noise. In one step, I need to take an evenly distributed number log u1 = (double) rand() / RAND_MAX
. Since it u1
can be zero, there is a risk of 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 do you choose the epsilon value? (There is in Fortran TINY
, but I don't know what to do in C).
Thank!
The complete code is attached:
#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);
}
}
source to share
You just need to make sure the result is rand()
not 0, so you don't have to do double conversion and division over and over.
int r = rand();
while (r == 0)
r = rand();
u1 = (double) rand() / RAND_MAX;
An even simpler solution is
u1 = (double)((unsigned)rand() + 1U)/((unsigned)RAND_MAX + 1U);
This way you don't need a loop anymore
However, you have to generate 53 bits for the double mantissa to get more correct precision instead of (usually) only 16 or 32 bits s rand()
. The sample is JavaRandom().nextDouble()
public double nextDouble() {
return (((long)next(26) << 27) + next(27))
/ (double)(1L << 53);
}
source to share
As @nhahtdh pointed out, the OP's current code is sufficient.
However, I suspect the error is "off by 1".
The OP's code will generate a range: 0.0 (exclusive) to 1.0 (inclusive). However, for this task, I would expect 0.0 (inclusive) to 1.0 (exclusive) with 0.0 removed.
If your code wants the range 0.0 (exclusive) to be 1.0 (excluding), consider
u1 = ((double) rand() + 0.5) / (RAND_MAX + 1.0);
source to share