rand.c

/* See: https://mumble.net/~campbell/2014/04/28/uniform-random-float */
/* See: https://allendowney.com/research/rand/downey07randfloat.pdf */

double randf_3a() {
    double d;
    uint64_t u = rand_next();
    uint64_t exp = __builtin_ctzll(u | (1 << 12)) - 1022ull;
    u = (((u >> 11) + 1) >> 1) - (exp << 52);
    memcpy(&d, &u, sizeof(d));
    return d;
}

double randf_3b() {
    double d;
    uint64_t u = rand_next();
    uint64_t exp = __builtin_ctzll(u) - 11ull;
    if ((int64_t)exp >= 0) exp = __builtin_ctzll(rand_next());
    u = (((u >> 11) + 1) >> 1) - ((exp - 1011ull) << 52);
    memcpy(&d, &u, sizeof(d));
    return d;
}