"Fossies" - the Fresh Open Source Software Archive

Member "txr-217/rand.c" (10 Jun 2019, 10865 Bytes) of package /linux/misc/txr-217.tar.bz2:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "rand.c" see the Fossies "Dox" file reference documentation.

    1 /* Copyright 2010-2019
    2  * Kaz Kylheku <kaz@kylheku.com>
    3  * Vancouver, Canada
    4  * All rights reserved.
    5  *
    6  * Redistribution and use in source and binary forms, with or without
    7  * modification, are permitted provided that the following conditions are met:
    8  *
    9  * 1. Redistributions of source code must retain the above copyright notice, this
   10  *    list of conditions and the following disclaimer.
   11  *
   12  * 2. Redistributions in binary form must reproduce the above copyright notice,
   13  *    this list of conditions and the following disclaimer in the documentation
   14  *    and/or other materials provided with the distribution.
   15  *
   16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
   17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
   18  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
   19  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
   20  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
   21  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
   22  * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
   23  * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
   24  * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
   25  * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
   26  */
   27 
   28 #include <string.h>
   29 #include <wctype.h>
   30 #include <stdarg.h>
   31 #include <wchar.h>
   32 #include <limits.h>
   33 #include <signal.h>
   34 #include "config.h"
   35 #if HAVE_UNISTD_H
   36 #include <unistd.h>
   37 #endif
   38 #include "lib.h"
   39 #include "signal.h"
   40 #include "unwind.h"
   41 #include "arith.h"
   42 #include "rand.h"
   43 #include "eval.h"
   44 
   45 #define random_warmup (deref(lookup_var_l(nil, random_warmup_s)))
   46 
   47 #if CHAR_BIT * SIZEOF_INT == 32
   48 typedef unsigned int rand32_t;
   49 #elif CHAR_BIT * SIZEOF_LONG == 32
   50 typedef unsigned long rand32_t;
   51 #endif
   52 
   53 /*
   54  * The algorithm here is WELL 512.
   55  * (Francois Panneton, Pierre L'Ecuyer.)
   56  */
   57 struct rand_state {
   58   rand32_t state[16];
   59   unsigned cur;
   60 };
   61 
   62 val random_state_s, random_state_var_s, random_warmup_s;
   63 
   64 static struct cobj_ops random_state_ops = cobj_ops_init(eq,
   65                                                         cobj_print_op,
   66                                                         cobj_destroy_free_op,
   67                                                         cobj_mark_op,
   68                                                         cobj_eq_hash_op);
   69 
   70 /* Source: bits from /dev/random on a Linux server */
   71 static rand32_t rand_tab[16] = {
   72   0x2C272ED6U, 0x4DBD5D69U, 0xC5482819U, 0x142AFCDEU,
   73   0xF7ABAEB0U, 0x454B47F1U, 0xFC85D2ADU, 0x1A9DB177U,
   74   0x2619231BU, 0x6B678AE8U, 0xAC450E78U, 0xA0A96B1CU,
   75   0x88A74E05U, 0xC1CBAEC2U, 0x8170BEADU, 0x29FAF776U
   76 };
   77 
   78 static val make_state(void)
   79 {
   80   struct rand_state *r = coerce(struct rand_state *, chk_malloc(sizeof *r));
   81   return cobj(coerce(mem_t *, r), random_state_s, &random_state_ops);
   82 }
   83 
   84 val random_state_p(val obj)
   85 {
   86   return typeof(obj) == random_state_s ? t : nil;
   87 }
   88 
   89 INLINE rand32_t *rstate(struct rand_state *r, int offs)
   90 {
   91   return &r->state[(r->cur + offs) % 16];
   92 }
   93 
   94 static rand32_t rand32(struct rand_state *r)
   95 {
   96   rand32_t s0 = *rstate(r, 0);
   97   rand32_t s9 = *rstate(r, 9);
   98   rand32_t s13 = *rstate(r, 13);
   99   rand32_t s15 = *rstate(r, 15);
  100 
  101   rand32_t r1 = s0 ^ (s0 << 16) ^ s13 ^ (s13 << 15);
  102   rand32_t r2 = s9 ^ (s9 >> 11);
  103 
  104   rand32_t ns0 = *rstate(r, 0) = r1 ^ r2;
  105   rand32_t ns15 = s15 ^ (s15 << 2) ^ r1 ^ (r1 << 18) ^ r2 ^ (r2 << 28) ^
  106                   ((ns0 ^ (ns0 << 5)) & 0xda442d24ul);
  107 
  108   *rstate(r, 15) = ns15;
  109   r->cur = (r->cur + 15) % 16;
  110   return ns15;
  111 }
  112 
  113 val make_random_state(val seed, val warmup)
  114 {
  115   val self = lit("make-random-state");
  116   val rs = make_state();
  117   int i = 0;
  118   struct rand_state *r = coerce(struct rand_state *,
  119                                 cobj_handle(self, rs, random_state_s));
  120 
  121   seed = default_null_arg(seed);
  122   warmup = default_null_arg(warmup);
  123 
  124   if (bignump(seed)) {
  125     mp_size dig, bit;
  126     mp_int *m = mp(seed);
  127 
  128     for (i = 0, dig = 0, bit = 0; i < 16 && dig < m->used; i++) {
  129       r->state[i] = (m->dp[dig] >> bit) & 0xFFFFFFFFul;
  130       bit += 32;
  131       if (bit >= MP_DIGIT_BIT)
  132         dig++, bit = 0;
  133     }
  134   } else if (fixnump(seed)) {
  135     cnum s = c_num(seed) & NUM_MAX;
  136 
  137     r->state[0] = s & 0xFFFFFFFFul;
  138     i = 1;
  139 #if CHAR_BIT * SIZEOF_PTR == 64
  140     s >>= 32;
  141     r->state[1] = s & 0xFFFFFFFFul;
  142     i = 2;
  143 #elif CHAR_BIT * SIZEOF_PTR > 64
  144 #error port me!
  145 #endif
  146   } else if (nilp(seed)) {
  147     val time = time_sec_usec();
  148     r->state[0] = convert(rand32_t, c_num(car(time)));
  149     r->state[1] = convert(rand32_t, c_num(cdr(time)));
  150 #if HAVE_UNISTD_H
  151     r->state[2] = convert(rand32_t, getpid());
  152     i = 3;
  153 #else
  154     i = 2;
  155 #endif
  156   } else if (random_state_p(seed)) {
  157     struct rand_state *rseed = coerce(struct rand_state *,
  158                                       cobj_handle(self, seed, random_state_s));
  159     *r = *rseed;
  160     return rs;
  161   } else if (vectorp(seed)) {
  162     if (length(seed) < num_fast(17))
  163       uw_throwf(error_s, lit("make-random-state: vector ~s too short"),
  164                 seed, nao);
  165 
  166     for (i = 0; i < 16; i++)
  167       r->state[i] = c_unum(seed->v.vec[i]);
  168 
  169     r->cur = c_num(seed->v.vec[i]);
  170     return rs;
  171   } else {
  172     uw_throwf(error_s, lit("make-random-state: seed ~s is not a number"),
  173               seed, nao);
  174   }
  175 
  176   while (i > 0 && r->state[i - 1] == 0)
  177     i--;
  178 
  179   for (; i < 16; i++)
  180     r->state[i] = rand_tab[i];
  181 
  182   r->cur = 0;
  183 
  184   {
  185     uses_or2;
  186     cnum wu = c_num(or2(warmup, random_warmup));
  187 
  188     for (i = 0; i < wu; i++)
  189       (void) rand32(r);
  190   }
  191 
  192   return rs;
  193 }
  194 
  195 val random_state_get_vec(val state)
  196 {
  197   val self = lit("random-state-get-vec");
  198   struct rand_state *r = coerce(struct rand_state *,
  199                                 cobj_handle(self,
  200                                             default_arg(state, random_state),
  201                                             random_state_s));
  202   int i;
  203   val vec = vector(num_fast(17), nil);
  204 
  205   for (i = 0; i < 16; i++)
  206     vec->v.vec[i] = normalize(bignum_from_uintptr(r->state[i]));
  207 
  208   vec->v.vec[i] = num(r->cur);
  209 
  210   return vec;
  211 }
  212 
  213 val random_fixnum(val state)
  214 {
  215   val self = lit("random-fixnum");
  216   struct rand_state *r = coerce(struct rand_state *,
  217                                 cobj_handle(self,
  218                                             default_arg(state, random_state),
  219                                             random_state_s));
  220   return num(rand32(r) & NUM_MAX);
  221 }
  222 
  223 static val random_float(val state)
  224 {
  225   val self = lit("random-float");
  226   struct rand_state *r = coerce(struct rand_state *,
  227                                 cobj_handle(self,
  228                                             default_arg(state, random_state),
  229                                             random_state_s));
  230   union hack {
  231     volatile double d;
  232     struct {
  233 #if HAVE_LITTLE_ENDIAN
  234       volatile rand32_t lo, hi;
  235 #else
  236       volatile rand32_t hi, lo;
  237 #endif
  238     } r;
  239   } h;
  240 
  241   h.r.lo = rand32(r);
  242   h.r.hi = (rand32(r) & 0xFFFFF) | (1023UL << 20);
  243 
  244   /* The least significant bit of the mantissa is always zero after
  245    * this subtraction, reducing us to 51 bits of precision.
  246    * Still; an attractive approach.
  247    */
  248   return flo(h.d - 1.0);
  249 }
  250 
  251 val random(val state, val modulus)
  252 {
  253   val self = lit("random");
  254   struct rand_state *r = coerce(struct rand_state *,
  255                                 cobj_handle(self, state, random_state_s));
  256   mp_int *m;
  257 
  258   if (bignump(modulus) && !mp_isneg(m = mp(modulus))) {
  259     ucnum bits = mp_count_bits(m) - mp_is_pow_two(m);
  260     ucnum rands_needed = (bits + 32 - 1) / 32;
  261     ucnum msb_rand_bits = bits % 32;
  262     rand32_t msb_rand_mask = convert(rand32_t, -1) >> (32 - msb_rand_bits);
  263     val out = make_bignum();
  264     mp_int *om = mp(out);
  265 
  266     for (;;) {
  267       ucnum i;
  268       for (i = 0; i < rands_needed; i++) {
  269         rand32_t rnd = rand32(r);
  270         mp_err mpe = MP_OKAY;
  271 #if MP_DIGIT_SIZE >= 4
  272         if (i > 0)
  273           mpe = mp_mul_2d(om, 32, om);
  274         else
  275           rnd &= msb_rand_mask;
  276         if (mpe == MP_OKAY)
  277           mpe = mp_add_d(om, rnd, om);
  278 #else
  279         if (i > 0)
  280           mpe = mp_mul_2d(om, 16, om);
  281         else
  282           rnd &= msb_rand_mask;
  283         if (mpe == MP_OKAY)
  284           mpe = mp_add_d(om, rnd & 0xFFFF, om);
  285         if (mpe == MP_OKAY)
  286           mpe = mp_mul_2d(om, 16, om);
  287         if (mpe == MP_OKAY)
  288           mp_add_d(om, rnd >> 16, om);
  289 #endif
  290         if (mpe != MP_OKAY)
  291           do_mp_error(self, mpe);
  292       }
  293       if (mp_cmp(om, m) != MP_LT) {
  294         mp_zero(om);
  295         continue;
  296       }
  297       break;
  298     }
  299 
  300     return normalize(out);
  301   } else if (fixnump(modulus)) {
  302     cnum m = c_num(modulus);
  303     if (m == 1) {
  304       return zero;
  305     } else if (m > 1) {
  306       int bits = highest_bit(m - 1);
  307 #if CHAR_BIT * SIZEOF_PTR >= 64
  308       ucnum rands_needed = (bits + 32 - 1) / 32;
  309 #endif
  310       ucnum msb_rand_bits = bits % 32;
  311       rand32_t msb_rand_mask = convert(rand32_t, -1) >> (32 - msb_rand_bits);
  312       for (;;) {
  313         cnum out = 0;
  314 #if CHAR_BIT * SIZEOF_PTR >= 64
  315         ucnum i;
  316 
  317         for (i = 0; i < rands_needed; i++) {
  318           rand32_t rnd = rand32(r);
  319           out <<= 32;
  320           if (i == 0)
  321             rnd &= msb_rand_mask;
  322           out |= rnd;
  323         }
  324 #else
  325         out = rand32(r) & msb_rand_mask;
  326 #endif
  327         if (out >= m)
  328           continue;
  329         return num(out);
  330       }
  331     }
  332   }
  333 
  334   uw_throwf(numeric_error_s, lit("random: invalid modulus ~s"),
  335       modulus, nao);
  336 }
  337 
  338 val rnd(val modulus, val state)
  339 {
  340   state = default_arg(state, random_state);
  341   return random(state, modulus);
  342 }
  343 
  344 void rand_compat_fixup(int compat_ver)
  345 {
  346   if (compat_ver <= 139) {
  347     loc l = lookup_var_l(nil, random_state_var_s);
  348     memset(rand_tab, 0xAA, sizeof rand_tab);
  349     if (compat_ver <= 114)
  350       random_state_s = random_state_var_s;
  351     set(l, make_random_state(num_fast(42), num_fast(8)));
  352   }
  353 }
  354 
  355 void rand_init(void)
  356 {
  357   random_state_var_s = intern(lit("*random-state*"), user_package);
  358   random_state_s = intern(lit("random-state"), user_package);
  359   random_warmup_s = intern(lit("*random-warmup*"), user_package);
  360 
  361   reg_var(random_state_var_s, make_random_state(num_fast(42), num_fast(8)));
  362   reg_var(random_warmup_s, num_fast(8));
  363 
  364   reg_fun(intern(lit("make-random-state"), user_package),
  365           func_n2o(make_random_state, 0));
  366   reg_fun(intern(lit("random-state-get-vec"), user_package),
  367           func_n1o(random_state_get_vec, 0));
  368   reg_fun(intern(lit("random-state-p"), user_package), func_n1(random_state_p));
  369   reg_fun(intern(lit("random-fixnum"), user_package), func_n1o(random_fixnum, 0));
  370   reg_fun(intern(lit("random-float"), user_package), func_n1o(random_float, 0));
  371   reg_fun(intern(lit("random"), user_package), func_n2(random));
  372   reg_fun(intern(lit("rand"), user_package), func_n2o(rnd, 1));
  373 }