mardi 30 décembre 2014

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

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).


Thanks!


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




Aucun commentaire:

Enregistrer un commentaire