Избегайте деления на ноль в 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? (В Фортране есть 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, так что вам не придется делать преобразование, чтобы удвоить и разделение снова и снова

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

u1 = (double) rand()/RAND_MAX; 

еще проще решение

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

Таким образом, вы не» t нужна петля больше

Однако вы должны сгенерировать 53 бита для двойной мантиссы, чтобы получить более правильную точность вместо (обычно) только 16 или 32 бит с rand(). Образец решения является 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 достаточен.

Тем не менее, я подозреваю, что ошибка «отключена на 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);