.. _program_listing_file_include_cmb_random.h:

Program Listing for File cmb_random.h
=====================================

|exhale_lsh| :ref:`Return to documentation for file <file_include_cmb_random.h>` (``include/cmb_random.h``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: c

   
   /*
    * 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 */