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 "cmb_assert.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);

[[maybe_unused]]
static inline double cmb_random(void)
{
    return ldexp((double)(cmb_random_sfc64() >> 11), -53);
}

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

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

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

[[maybe_unused]]
static inline double cmb_random_logistic(const double m, const double s)
{
    cmb_assert_release(s > 0.0);
    const double x = cmb_random();

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

[[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(uint64_t u_cand_x);
[[maybe_unused]]
static inline double cmb_random_std_exponential(void)
{
    const uint64_t u_cand_x = cmb_random_sfc64();
    const uint8_t idx = u_cand_x & 0xff;
    double r = (idx <= cmi_random_exp_zig_max) ?
        cmi_random_exp_zig_pdf_x[idx] * (double) u_cand_x :
        cmi_random_exp_not_hot(u_cand_x);

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

[[maybe_unused]]
static inline double cmb_random_exponential(const double mean)
{
    cmb_assert_release(mean > 0.0);

    double r = mean * cmb_random_std_exponential();

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

[[maybe_unused]]
static inline double cmb_random_erlang(const unsigned k, const double m)
{
    cmb_assert_release(k > 0u);
    cmb_assert_release(m > 0.0);

    double x = 0.0;
    for (unsigned i = 0u; i < k; i++) {
        x += cmb_random_exponential(m);
    }

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

[[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(unsigned n,
                                          const double *ma,
                                          const double *pa);

extern double cmb_random_std_gamma(double shape);

[[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 = (shape >= 1.0) ?
        scale * cmb_random_std_gamma(shape) :
        scale * (cmb_random_std_gamma(shape + 1.0)
                 * pow(cmb_random(), 1.0 / shape));

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

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

    const double x = cmb_random_std_gamma(a);
    const double y = cmb_random_std_gamma(b);
    const double r = x / (x + y);

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

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

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

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

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

    const double x = mode / pow(cmb_random(), 1.0 / shape);

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

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

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

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

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

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

[[maybe_unused]]
static inline unsigned cmb_random_bernoulli(const double p)
{
    cmb_assert_release((p >= 0.0) && (p <= 1.0));

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

extern unsigned cmb_random_geometric(double p);

extern unsigned cmb_random_binomial(unsigned n, double p);

extern unsigned cmb_random_negative_binomial(unsigned m, double p);

[[maybe_unused]]
static inline unsigned cmb_random_pascal(const unsigned m, const double p)
{
    return cmb_random_negative_binomial(m, p);
}

extern unsigned cmb_random_poisson(double r);

[[maybe_unused]]
static inline long cmb_random_dice(const long a, const long b)
{
    cmb_assert (a < b);

    const double x = (double)(b - a + 1) * cmb_random();
    return (long)(floor((double)a + x));
}

extern unsigned cmb_random_loaded_dice(unsigned n, const double *pa);

struct cmb_random_alias {
    unsigned n;
    uint64_t *uprob;
    unsigned *alias;
};

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

[[maybe_unused]]
static inline unsigned cmb_random_alias_sample(const struct cmb_random_alias *ap)
{
    cmb_assert_release(ap != NULL);

    const unsigned idx = (unsigned) (floor(ap->n * cmb_random()));
    const bool c = (cmb_random_sfc64() >= ap->uprob[idx]);
    const unsigned r = (c) ? ap->alias[idx] : idx;

    cmb_assert_debug(r < ap->n);
    return r;
}

extern void cmb_random_alias_destroy(struct cmb_random_alias *ap);

#endif /* CIMBA_CMB_RANDOM_H */