Line data    Source code

       1           1 : /**
       2             :  * @file lib/random/random.c
       3             :  *
       4             :  * @brief Random Number Generators
       5             :  *
       6             :  * Piece-Wise Deterministic Random Number Generators.
       7             :  *
       8             :  * SPDX-FileCopyrightText: 2008-2021 HPDCS Group <rootsim@googlegroups.com>
       9             :  * SPDX-License-Identifier: GPL-3.0-only
      10             :  */
      11             : #include <lib/random/random.h>
      12             : 
      13             : #include <core/intrinsics.h>
      14             : #include <lib/lib_internal.h>
      15             : #include <lib/random/xoroshiro.h>
      16             : 
      17             : #include <math.h>
      18             : #include <memory.h>
      19             : 
      20           0 : void random_lib_lp_init(void)
      21             : {
      22             :         uint64_t seed = global_config.prng_seed;
      23             :         lp_id_t lid = lp_id_get();
      24             :         struct lib_ctx *ctx = lib_ctx_get();
      25             :         random_init(ctx->rng_s, lid, seed);
      26             :         ctx->unif = NAN;
      27             : }
      28             : 
      29           0 : double Random(void)
      30             : {
      31             :         struct lib_ctx *ctx = lib_ctx_get();
      32             :         uint64_t u_val = random_u64(ctx->rng_s);
      33             :         if (unlikely(!u_val))
      34             :                 return 0.0;
      35             : 
      36             :         double ret = 0.0;
      37             :         unsigned lzs = intrinsics_clz(u_val) + 1;
      38             :         u_val <<= lzs;
      39             :         u_val >>= 12;
      40             : 
      41             :         uint64_t exp = 1023 - lzs;
      42             :         u_val |= exp << 52;
      43             : 
      44             :         memcpy(&ret, &u_val, sizeof(double));
      45             :         return ret;
      46             : }
      47             : 
      48           0 : uint64_t RandomU64(void)
      49             : {
      50             :         struct lib_ctx *ctx = lib_ctx_get();
      51             :         return random_u64(ctx->rng_s);
      52             : }
      53             : 
      54             : /**
      55             :  * Return a random number according to an Exponential distribution.
      56             :  * The mean value of the distribution must be passed as the mean value.
      57             :  *
      58             :  * @param mean Mean value of the distribution
      59             :  * @return A random number
      60             :  */
      61           1 : double Expent(double mean)
      62             : {
      63             :         if (unlikely(mean < 0)) {
      64             :                 log_log(LOG_WARN, "Passed a negative mean into Expent()");
      65             :         }
      66             :         return -mean * log(1 - Random());
      67             : }
      68             : 
      69             : /**
      70             :  * Return a random number according to a Standard Normal Distribution
      71             :  *
      72             :  * @return A random number
      73             :  */
      74           1 : double Normal(void)
      75             : {
      76             :         struct lib_ctx *ctx = lib_ctx_get();
      77             :         if (isnan(ctx->unif)) {
      78             :                 double v1, v2, rsq;
      79             :                 do {
      80             :                         v1 = 2.0 * Random() - 1.0;
      81             :                         v2 = 2.0 * Random() - 1.0;
      82             :                         rsq = v1 * v1 + v2 * v2;
      83             :                 } while (rsq >= 1.0 || rsq == 0);
      84             : 
      85             :                 double fac = sqrt(-2.0 * log(rsq) / rsq);
      86             : 
      87             :                 // Perform Box-Muller transformation to get two normal deviates.
      88             :                 // Return one and save the other for next time.
      89             :                 ctx->unif = v1 * fac;
      90             :                 return v2 * fac;
      91             :         } else {
      92             :                 // A deviate is already available
      93             :                 double ret = ctx->unif;
      94             :                 ctx->unif = NAN;
      95             :                 return ret;
      96             :         }
      97             : }
      98             : 
      99           0 : int RandomRange(int min, int max)
     100             : {
     101             :         return (int)floor(Random() * (max - min + 1)) + min;
     102             : }
     103             : 
     104           0 : int RandomRangeNonUniform(int x, int min, int max)
     105             : {
     106             :         return (((RandomRange(0, x) | RandomRange(min, max))) %
     107             :                         (max - min + 1)) + min;
     108             : }
     109             : 
     110             : /**
     111             :  * Return a number in according to a Gamma Distribution of Integer Order ia,
     112             :  * a waiting time to the ia-th event in a Poisson process of unit mean.
     113             :  *
     114             :  * @author D. E. Knuth
     115             :  * @param ia Integer Order of the Gamma Distribution
     116             :  * @return A random number
     117             :  */
     118           1 : double Gamma(unsigned ia)
     119             : {
     120             :         if (unlikely(ia < 1)) {
     121             :                 log_log(LOG_WARN, "Gamma distribution must have a ia "
     122             :                                 "value >= 1. Defaulting to 1...");
     123             :                 ia = 1;
     124             :         }
     125             : 
     126             :         double x;
     127             : 
     128             :         if (ia < 6) {
     129             :                 // Use direct method, adding waiting times
     130             :                 x = 1.0;
     131             :                 while (ia--)
     132             :                         x *= 1 - Random();
     133             :                 x = -log(x);
     134             :         } else {
     135             :                 double am = ia - 1;
     136             :                 double v1, v2, e, y, s;
     137             :                 // Use rejection method
     138             :                 do {
     139             :                         do {
     140             :                                 do {
     141             :                                         v1 = Random();
     142             :                                         v2 = 2.0 * Random() - 1.0;
     143             :                                 } while (v1 * v1 + v2 * v2 > 1.0);
     144             : 
     145             :                                 y = v2 / v1;
     146             :                                 s = sqrt(2.0 * am + 1.0);
     147             :                                 x = s * y + am;
     148             :                         } while (x < 0.0);
     149             : 
     150             :                         e = (1.0 + y * y) * exp(am * log(x / am) - s * y);
     151             :                 } while (Random() > e);
     152             :         }
     153             : 
     154             :         return x;
     155             : }
     156             : 
     157             : /**
     158             :  * Return the waiting time to the next event in a Poisson process of unit mean.
     159             :  *
     160             :  * @return A random number
     161             :  */
     162           1 : double Poisson(void)
     163             : {
     164             :         return -log(1 - Random());
     165             : }
     166             : 
     167             : /**
     168             :  * Return a random sample from a Zipf distribution.
     169             :  * Based on the rejection method by Luc Devroye for sampling:
     170             :  * "Non-Uniform Random Variate Generation, page 550, Springer-Verlag, 1986
     171             :  *
     172             :  * @param skew The skew of the distribution
     173             :  * @param limit The largest sample to retrieve
     174             :  * @return A random number
     175             :  */
     176           1 : unsigned Zipf(double skew, unsigned limit)
     177             : {
     178             :         double b = pow(2., skew - 1.);
     179             :         double x, t;
     180             :         do {
     181             :                 x = floor(pow(Random(), -1. / skew - 1.));
     182             :                 t = pow(1. + 1. / x, skew - 1.);
     183             :         } while (x > limit || Random() * x * (t - 1.) * b > t * (b - 1.));
     184             :         return x;
     185             : }