Program Listing for File cmb_random.h

↰ Return to documentation for file (include/cmb_random.h)

/*
 * Copyright (c) Asbjørn M. Bonvik 1994, 1995, 2025-26.
 *
 * The normal and exponential distributions below are based on code at
 *      https://github.com/cd-mcfarland/fast_prng
 *      Copyright (c) Chris D McFarland 2025.
 *      Used with permission by author.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#ifndef CIMBA_CMB_RANDOM_H
#define CIMBA_CMB_RANDOM_H

#include <inttypes.h>
#include <math.h>
#include <stdbool.h>
#include <stddef.h>

#include "cmb_assert.h"

#include "cmi_memregistry.h"

extern void cmb_random_initialize(uint64_t seed);

extern void cmb_random_terminate(void);

extern uint64_t cmb_random_hwseed(void);

extern uint64_t cmb_random_curseed(void);

extern uint64_t cmb_random_sfc64(void);

extern uint64_t cmb_random_fmix64(uint64_t seed, uint64_t nonce);

CMB_MAYBE_UNUSED
static inline double cmb_random(void)
{
    const uint64_t x = cmb_random_sfc64();
    const double r = (double)(x >> 11) * 0x1.0p-53;

    cmb_assert_debug((r >= 0.0) && (r < 1.0));
    return r;
}

CMB_MAYBE_UNUSED
static inline double cmb_random_uniform(const double min, const double max)
{
    cmb_assert_release(min < max);

    const double r = min + (max - min) * cmb_random();
    cmb_assert_debug((r >= min) && (r < max));

    return r;
}

extern double cmb_random_triangular(double min, double mode, double max);

extern const uint8_t cmi_random_nor_zig_max;
extern const double cmi_random_nor_zig_pdf_x[];
extern double cmi_random_nor_not_hot(int64_t i_cand_x);
CMB_MAYBE_UNUSED
static inline double cmb_random_std_normal(void)
{
    uint64_t bits = cmb_random_sfc64();
    const int64_t i_cand_x = *(int64_t *) &bits;
    const uint8_t idx = i_cand_x & 0xFF;

    return (idx <= cmi_random_nor_zig_max) ?
        cmi_random_nor_zig_pdf_x[idx] * (double) i_cand_x :
        cmi_random_nor_not_hot(i_cand_x);
}

CMB_MAYBE_UNUSED
static inline double cmb_random_normal(const double mu, const double sigma)
{
    cmb_assert_release(sigma > 0.0);

    return mu + sigma * cmb_random_std_normal();
}

CMB_MAYBE_UNUSED
static inline double cmb_random_lognormal(const double m, const double s)
{
    cmb_assert_release(s > 0.0);
    const double r = exp(cmb_random_normal(m, s));

    cmb_assert_debug(r >= 0.0);
    return r;
}

CMB_MAYBE_UNUSED
static inline double cmb_random_logistic(const double m, const double s)
{
    cmb_assert_release(s > 0.0);

    /* Loop to protect against the extremely unlikely case of drawing 0.0 */
    double x = 0.0;
    while (x == 0.0) {
        x = cmb_random();
    }

    return m + s * log(x / (1.0 - x));
}

CMB_MAYBE_UNUSED
static inline double cmb_random_cauchy(const double mode, const double scale)
{
    cmb_assert_release(scale > 0.0);

    const double x = cmb_random_std_normal();
    double y;
    while ((y = cmb_random_std_normal()) == 0.0) {}

    return mode + scale * x / y;
 }

extern const uint8_t cmi_random_exp_zig_max;
extern const double cmi_random_exp_zig_pdf_x[];
extern double cmi_random_exp_not_hot(int64_t i_cand_x);
CMB_MAYBE_UNUSED
static inline double cmb_random_std_exponential(void)
{
    uint64_t bits = cmb_random_sfc64();
    /* Clear the sign bit for 63 useful bits */
    const int64_t i_cand_x = (*(int64_t *) &bits) & INT64_MAX;
    const uint8_t idx = i_cand_x & 0xFF;
    const double r = (idx <= cmi_random_exp_zig_max) ?
                        cmi_random_exp_zig_pdf_x[idx] * (double) i_cand_x :
                        cmi_random_exp_not_hot(i_cand_x);

    cmb_assert_debug(r >= 0.0);
    return r;
}

CMB_MAYBE_UNUSED
static inline double cmb_random_exponential(const double mean)
{
    cmb_assert_release(mean > 0.0);

    const double r = mean * cmb_random_std_exponential();

    cmb_assert_debug(r >= 0.0);
    return r;
}

extern double cmb_random_erlang(const unsigned k, const double m);

CMB_MAYBE_UNUSED
static inline double cmb_random_hypoexponential(const unsigned n, const double *ma)
{
    cmb_assert_release(n > 0);
    cmb_assert_release(ma != NULL);

    double x = 0.0;
    for (unsigned i = 0; i < n; i++) {
        cmb_assert_release(ma[i] > 0.0);
        x += cmb_random_exponential(ma[i]);
    }

    cmb_assert_debug(x >= 0.0);
    return x;
}

extern double cmb_random_hyperexponential(uint64_t n,
                                          const double *ma,
                                          const double *pa);

extern double cmb_random_std_gamma(double shape);

CMB_MAYBE_UNUSED
static inline double cmb_random_gamma(const double shape, const double scale)
{
    cmb_assert_release(shape > 0.0);
    cmb_assert_release(scale > 0.0);

    const double r = scale * cmb_random_std_gamma(shape);

    cmb_assert_debug(r >= 0.0);
    return r;
}

CMB_MAYBE_UNUSED
static inline double cmb_random_std_beta(const double a, const double b)
{
    cmb_assert_release(a > 0.0);
    cmb_assert_release(b > 0.0);

    double x;
    double y;
    do {
        /* Guard against the exceedingly unlikely double zero outcome */
        x = cmb_random_std_gamma(a);
        y = cmb_random_std_gamma(b);
    } while ((x == 0.0) && (y == 0.0));

    const double r = x / (x + y);

    cmb_assert_debug((r >= 0.0) && (r <= 1.0));
    return r;
}

CMB_MAYBE_UNUSED
static inline double cmb_random_beta(const double a, const double b,
                                     const double min, const double max)
{
    cmb_assert_release(a > 0.0);
    cmb_assert_release(b > 0.0);
    cmb_assert_release(min < max);

    const double x = min + (max - min) * cmb_random_std_beta(a, b);

    cmb_assert_debug((x >= min) && (x <= max));
    return x;
}

extern double cmb_random_PERT_mod(double min,
                                  double mode,
                                  double max,
                                  double lambda);

CMB_MAYBE_UNUSED
static inline double cmb_random_PERT(const double min,
                                     const double mode,
                                     const double max)
{
    cmb_assert_release(min < mode);
    cmb_assert_release(mode < max);

    const double x = cmb_random_PERT_mod(min, mode, max, 4.0);

    cmb_assert_debug((x >= min) && (x <= max));
    return x;
}

CMB_MAYBE_UNUSED
static inline double cmb_random_weibull(const double shape, const double scale)
{
    cmb_assert_release(shape > 0.0);
    cmb_assert_release(scale > 0.0);

    const double u = cmb_random_exponential(1.0);
    const double x = scale * pow(u, 1.0 / shape);

    cmb_assert_debug(x >= 0.0);
    return x;
}

CMB_MAYBE_UNUSED
static inline double cmb_random_pareto(const double shape, const double mode)
{
    cmb_assert_release(shape > 0.0);
    cmb_assert_release(mode > 0.0);

    /* Generate uniform sample on (0.0, 1.0] */
    const double y = 1.0 - cmb_random();
    const double x = mode / pow(y, 1.0 / shape);

    cmb_assert_debug(x >= mode);
    return x;
}

CMB_MAYBE_UNUSED
static inline double cmb_random_chisquared(const double k)
{
    cmb_assert_release(k > 0.0);

    const double x = cmb_random_gamma(k / 2.0, 2.0);

    cmb_assert_debug(x >= 0.0);
    return x;
}

CMB_MAYBE_UNUSED
static inline double cmb_random_F_dist(const double a, const double b)
{
    cmb_assert_release(a > 0.0);
    cmb_assert_release(b > 0.0);

    const double x = cmb_random_chisquared(a) / a;
    double y;
    while ((y = cmb_random_chisquared(b) / b) == 0.0) {}

    const double r = x / y;

    cmb_assert_debug(r >= 0.0);
    return r;
}

CMB_MAYBE_UNUSED
static inline double cmb_random_std_t_dist(const double v)
{
    cmb_assert_release(v > 0.0);

    const double x = cmb_random_std_normal();
    double y;
    while ((y = cmb_random_chisquared(v)) == 0.0) {}

    const double r = x / sqrt( y / v);
    return r;
}

CMB_MAYBE_UNUSED
static inline double cmb_random_t_dist(const double m,
                                const double s,
                                const double v)
{
    cmb_assert_release(s > 0.0);
    cmb_assert_release(v > 0.0);

    const double x = m + s * cmb_random_std_t_dist(v);
    return x;
}

CMB_MAYBE_UNUSED
static inline double cmb_random_rayleigh(const double s)
{
    cmb_assert_release(s > 0.0);

    const double x = cmb_random_normal(0.0, s);
    const double y = cmb_random_normal(0.0, s);
    const double r = sqrt(x * x + y * y);

    cmb_assert_debug(r >= 0.0);
    return r;
}

extern int cmb_random_flip(void);

CMB_MAYBE_UNUSED
static inline unsigned int cmb_random_bernoulli(const double p)
{
    cmb_assert_release((p >= 0.0) && (p <= 1.0));

    return (cmb_random() < p) ? 1 : 0;
}

extern uint64_t cmb_random_geometric(double p);

extern uint64_t cmb_random_binomial(uint64_t n, double p);

extern uint64_t cmb_random_negative_binomial(uint64_t m, double p);

CMB_MAYBE_UNUSED
static inline uint64_t cmb_random_pascal(const uint64_t m, const double p)
{
    return cmb_random_negative_binomial(m, p);
}

extern uint64_t cmi_random_poisson_chopdown(double r);
extern uint64_t cmi_random_poisson_ptrd(double r);
extern uint64_t cmb_random_poisson(double r);

extern uint64_t cmb_random_discrete_uniform(uint64_t n);

CMB_MAYBE_UNUSED
static inline int64_t cmb_random_dice(const int64_t a, const int64_t b)
{
    cmb_assert_release(a < b);

    const uint64_t span = (uint64_t)b - (uint64_t)a;
    cmb_assert_release(span != UINT64_MAX);
    const uint64_t x = cmb_random_discrete_uniform(span + 1u);
    const int64_t r = (int64_t)((uint64_t)a + x);

    cmb_assert_debug((r >= a) && (r <= b));
    return r;
}

extern uint64_t cmb_random_discrete_nonuniform(uint64_t n, const double *pa);

CMB_MAYBE_UNUSED
static inline int64_t cmb_random_loaded_dice(const int64_t a, const int64_t b,
                                             const double *pa)
{
    cmb_assert_release(a < b);
    cmb_assert_release(pa != NULL);

    /* Avoid potential overflow in edge cases, assuming a two's complement machine */
    const uint64_t span = (uint64_t)b - (uint64_t)a;
    cmb_assert_release(span < UINT64_MAX);
    const uint64_t n = span + 1u;
    const uint64_t offset = cmb_random_discrete_nonuniform(n, pa);
    const int64_t r = (int64_t)((uint64_t)a + offset);

    cmb_assert_debug((r >= a) && (r <= b));
    return r;
}

struct cmb_random_alias {
    uint64_t n;
    uint64_t *uprob;
    uint64_t *alias;
    struct cmi_memregistry_item destroy;
};

extern struct cmb_random_alias *cmb_random_alias_create(uint64_t n,
                                                        const double *pa);

extern uint64_t cmb_random_alias_sample(const struct cmb_random_alias *ap);

extern void cmb_random_alias_destroy(struct cmb_random_alias *ap);

#endif /* CIMBA_CMB_RANDOM_H */