random_normal.c

浏览该文件的文档.

#include <math.h>
#include "random_internal.h"

 
#define ZIG_N 256
#define ZIG_R 3.6541528853610088
 
static const double zig_xtab[ZIG_N + 1] = {
#include "ziggurat_table.inc"
};
static inline double zig_pdf(double x) {
    return exp(-0.5 * x * x);
}
 

static double zig_sample_tail(uint64_t* state) {
    double x, y, u1, u2;
    for (;;) {
        do {
            u1 = u64_to_double01(xoshiro256ss_next(state));
        } while (u1 == 0.0);
        do {
            u2 = u64_to_double01(xoshiro256ss_next(state));
        } while (u2 == 0.0);
 
        x = -log(u1) / ZIG_R;
        y = -log(u2);
 
        if (2.0 * y >= x * x) {
            return x + ZIG_R;
        }
    }
}

double lmmc_rng_standard_normal(uint64_t* state) {
    uint64_t u, u2;
    int i, sign;
    double x;
 
 
    for (;;) {
        u = xoshiro256ss_next(state);
        i = (int)(u & 0xFF);  
        sign = (u & 0x100) ? -1 : 1;  

        u2 = xoshiro256ss_next(state);
        x = u64_to_double01(u2) * zig_xtab[i];
 
        if (x < zig_xtab[i + 1]) {
            return sign * x;
        }

        if (i == 0) {
            if (x < zig_xtab[1]) {
                return sign * x;
            }
            double tail = zig_sample_tail(state);
            return sign * tail;
        }

        {
            double f_x = zig_pdf(x);
            double f_outer = zig_pdf(zig_xtab[i]);     
            double f_inner = zig_pdf(zig_xtab[i + 1]); 
            double u_wedge = u64_to_double01(xoshiro256ss_next(state));
 
            if (u_wedge * (f_inner - f_outer) < (f_x - f_outer)) {
                return sign * x;
            }
        }
    }
}
 
lmmc_status_t lmmc_rng_normal(
    lmmc_rng_t* rng,
    lmmc_real_t mean,
    lmmc_real_t stddev,
    lmmc_real_t* out_value)
{
    if (rng == NULL || out_value == NULL) {
        return LMMC_STATUS_INVALID_ARGUMENT;
    }
    if (!isfinite(mean) || !isfinite(stddev) || stddev <= 0.0) {
        return LMMC_STATUS_INVALID_ARGUMENT;
    }
 
    *out_value = mean + stddev * lmmc_rng_standard_normal(rng->state);
    return LMMC_STATUS_OK;
}