src/rand/rand.c

   1 /*
   2  * *****************************************************************************
   3  *
   4  * Copyright (c) 2018-2019 Gavin D. Howard and contributors.
   5  *
   6  * All rights reserved.
   7  *
   8  * Redistribution and use in source and binary forms, with or without
   9  * modification, are permitted provided that the following conditions are met:
  10  *
  11  * * Redistributions of source code must retain the above copyright notice, this
  12  *   list of conditions and the following disclaimer.
  13  *
  14  * * Redistributions in binary form must reproduce the above copyright notice,
  15  *   this list of conditions and the following disclaimer in the documentation
  16  *   and/or other materials provided with the distribution.
  17  *
  18  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  19  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  20  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  21  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
  22  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  23  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  24  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  25  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  26  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  27  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  28  * POSSIBILITY OF SUCH DAMAGE.
  29  *
  30  * *****************************************************************************
  31  *
  32  * Parts of this code are adapted from the following:
  33  *
  34  * PCG, A Family of Better Random Number Generators.
  35  *
  36  * You can find the original source code at:
  37  *   https://github.com/imneme/pcg-c
  38  *
  39  * -----------------------------------------------------------------------------
  40  *
  41  * Parts of this code are also under the following license:
  42  *
  43  * Copyright (c) 2014-2017 Melissa O'Neill and PCG Project contributors
  44  *
  45  * Permission is hereby granted, free of charge, to any person obtaining a copy
  46  * of this software and associated documentation files (the "Software"), to deal
  47  * in the Software without restriction, including without limitation the rights
  48  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  49  * copies of the Software, and to permit persons to whom the Software is
  50  * furnished to do so, subject to the following conditions:
  51  *
  52  * The above copyright notice and this permission notice shall be included in
  53  * all copies or substantial portions of the Software.
  54  *
  55  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  56  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  57  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  58  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  59  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  60  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  61  * SOFTWARE.
  62  *
  63  * *****************************************************************************
  64  *
  65  * Code for the RNG.
  66  *
  67  */
  68
  69 #if BC_ENABLE_EXTRA_MATH
  70
  71 #include <assert.h>
  72 #include <stdlib.h>
  73 #include <string.h>
  74 #include <time.h>
  75 #include <fcntl.h>
  76 #include <unistd.h>
  77
  78 #include <status.h>
  79 #include <num.h>
  80 #include <rand.h>
  81 #include <vm.h>
  82
  83 #if !BC_RAND_BUILTIN
  84
  85 static BcRandState bc_rand_addition(uint_fast64_t a, uint_fast64_t b) {
  86
  87         BcRandState res;
  88
  89         res.lo = a + b;
  90         res.hi = (res.lo < a);
  91
  92         return res;
  93 }
  94
  95 static BcRandState bc_rand_addition2(BcRandState a, BcRandState b) {
  96
  97         BcRandState temp, res;
  98
  99         res = bc_rand_addition(a.lo, b.lo);
 100         temp = bc_rand_addition(a.hi, b.hi);
 101         res.hi += temp.lo;
 102
 103         return res;
 104 }
 105
 106 static BcRandState bc_rand_multiply(uint_fast64_t a, uint_fast64_t b) {
 107
 108         uint_fast64_t al, ah, bl, bh, c0, c1, c2, c3;
 109         BcRandState carry, res;
 110
 111         al = BC_RAND_TRUNC32(a);
 112         ah = BC_RAND_CHOP32(a);
 113         bl = BC_RAND_TRUNC32(b);
 114         bh = BC_RAND_CHOP32(b);
 115
 116         c0 = al * bl;
 117         c1 = al * bh;
 118         c2 = ah * bl;
 119         c3 = ah * bh;
 120
 121         carry = bc_rand_addition(c1, c2);
 122
 123         res = bc_rand_addition(c0, (BC_RAND_TRUNC32(carry.lo)) << 32);
 124         res.hi += BC_RAND_CHOP32(carry.lo) + c3 + (carry.hi << 32);
 125
 126         return res;
 127 }
 128
 129 static BcRandState bc_rand_multiply2(BcRandState a, BcRandState b) {
 130
 131         BcRandState c0, c1, c2, carry;
 132
 133         c0 = bc_rand_multiply(a.lo, b.lo);
 134         c1 = bc_rand_multiply(a.lo, b.hi);
 135         c2 = bc_rand_multiply(a.hi, b.lo);
 136
 137         carry = bc_rand_addition2(c1, c2);
 138         carry.hi = carry.lo;
 139         carry.lo = 0;
 140
 141         return bc_rand_addition2(c0, carry);
 142 }
 143
 144 #endif // BC_RAND_BUILTIN
 145
 146 static void bc_rand_setModified(BcRNGData *r) {
 147
 148 #if BC_RAND_BUILTIN
 149         r->inc |= (BcRandState) 1UL;
 150 #else // BC_RAND_BUILTIN
 151         r->inc.lo |= (uint_fast64_t) 1UL;
 152 #endif // BC_RAND_BUILTIN
 153 }
 154
 155 static void bc_rand_clearModified(BcRNGData *r) {
 156
 157 #if BC_RAND_BUILTIN
 158         r->inc &= ~((BcRandState) 1UL);
 159 #else // BC_RAND_BUILTIN
 160         r->inc.lo &= ~(1UL);
 161 #endif // BC_RAND_BUILTIN
 162 }
 163
 164 static void bc_rand_copy(BcRNGData *d, BcRNGData *s) {
 165         bool unmod = BC_RAND_NOTMODIFIED(d);
 166         memcpy(d, s, sizeof(BcRNGData));
 167         if (!unmod) bc_rand_setModified(d);
 168         else if (!BC_RAND_NOTMODIFIED(s)) bc_rand_clearModified(d);
 169 }
 170
 171 static ulong bc_rand_frand(void *ptr) {
 172
 173         ulong buf[1];
 174         int fd;
 175         ssize_t nread;
 176
 177         assert(ptr != NULL);
 178
 179         fd = *((int*) ptr);
 180
 181         nread = read(fd, buf, sizeof(ulong));
 182
 183         if (BC_ERR(nread != sizeof(ulong))) bc_vm_err(BC_ERROR_FATAL_IO_ERR);
 184
 185         return *((ulong*) buf);
 186 }
 187
 188 static ulong bc_rand_rand(void *ptr) {
 189
 190         size_t i;
 191         ulong res = 0;
 192
 193         BC_UNUSED(ptr);
 194
 195         for (i = 0; i < sizeof(ulong); ++i)
 196                 res |= ((ulong) (rand() & BC_RAND_SRAND_BITS)) << (i * CHAR_BIT);
 197
 198         return res;
 199 }
 200
 201 static BcRandState bc_rand_inc(BcRNGData *r) {
 202
 203         BcRandState inc;
 204
 205 #if BC_RAND_BUILTIN
 206         inc = r->inc | 1;
 207 #else // BC_RAND_BUILTIN
 208         inc.lo = r->inc.lo | 1;
 209         inc.hi = r->inc.hi;
 210 #endif // BC_RAND_BUILTIN
 211
 212         return inc;
 213 }
 214
 215 static void bc_rand_setInc(BcRNGData *r) {
 216
 217 #if BC_RAND_BUILTIN
 218         r->inc <<= 1UL;
 219 #else // BC_RAND_BUILTIN
 220         r->inc.hi <<= 1UL;
 221         r->inc.hi |= (r->inc.lo & (1UL << (BC_LONG_BIT - 1))) >> (BC_LONG_BIT - 1);
 222         r->inc.lo <<= 1UL;
 223 #endif // BC_RAND_BUILTIN
 224 }
 225
 226 static void bc_rand_seedState(BcRandState *state, ulong val1, ulong val2) {
 227
 228 #if BC_RAND_BUILTIN
 229         *state = ((BcRandState) val1) | ((BcRandState) val2) << (BC_LONG_BIT);
 230 #else // BC_RAND_BUILTIN
 231         state->lo = val1;
 232         state->hi = val2;
 233 #endif // BC_RAND_BUILTIN
 234 }
 235
 236 static void bc_rand_seedRNG(BcRNGData *r, ulong state1, ulong state2,
 237                             ulong inc1, ulong inc2)
 238 {
 239         bc_rand_seedState(&r->state, state1, state2);
 240         bc_rand_seedState(&r->inc, inc1, inc2);
 241         bc_rand_setInc(r);
 242 }
 243
 244 static void bc_rand_fill(BcRNGData *r, BcRandUlong fulong, void *ptr) {
 245
 246         ulong state1, state2, inc1, inc2;
 247
 248         state1 = fulong(ptr);
 249         state2 = fulong(ptr);
 250
 251         inc1 = fulong(ptr);
 252         inc2 = fulong(ptr);
 253
 254         bc_rand_seedRNG(r, state1, state2, inc1, inc2);
 255 }
 256
 257 static void bc_rand_step(BcRNGData *r) {
 258         BcRandState temp = bc_rand_mul2(r->state, bc_rand_multiplier);
 259         r->state = bc_rand_add2(temp, bc_rand_inc(r));
 260 }
 261
 262 static BcRand bc_rand_output(BcRNGData *r) {
 263         return BC_RAND_ROT(BC_RAND_FOLD(r->state), BC_RAND_ROTAMT(r->state));
 264 }
 265
 266 static void bc_rand_seedZeroes(BcRNG *r, BcRNGData *rng, size_t idx) {
 267
 268         BcRNGData *rng2;
 269
 270         if (r->v.len <= idx) return;
 271
 272         rng2 = bc_vec_item_rev(&r->v, idx);
 273
 274         if (BC_RAND_ZERO(rng2)) {
 275                 size_t i;
 276                 for (i = 1; i < r->v.len; ++i)
 277                         bc_rand_copy(bc_vec_item_rev(&r->v, i), rng);
 278         }
 279 }
 280
 281 static void bc_rand_srand(BcRNGData *rng) {
 282
 283         int fd;
 284
 285         BC_SIG_LOCK;
 286
 287         fd = open("/dev/urandom", O_RDONLY);
 288
 289         if (BC_NO_ERR(fd >= 0)) {
 290                 bc_rand_fill(rng, bc_rand_frand, &fd);
 291                 close(fd);
 292         }
 293
 294         while (BC_ERR(BC_RAND_ZERO(rng))) bc_rand_fill(rng, bc_rand_rand, NULL);
 295
 296         BC_SIG_UNLOCK;
 297 }
 298
 299 static void bc_rand_propagate(BcRNG *r, BcRNGData *rng) {
 300
 301         if (r->v.len <= 1) return;
 302
 303         if (BC_RAND_NOTMODIFIED(rng)) {
 304
 305                 size_t i;
 306                 bool go = true;
 307
 308                 for (i = 1; go && i < r->v.len; ++i) {
 309                         BcRNGData *rng2 = bc_vec_item_rev(&r->v, i);
 310                         go = BC_RAND_NOTMODIFIED(rng2);
 311                         bc_rand_copy(rng2, rng);
 312                 }
 313
 314                 bc_rand_seedZeroes(r, rng, i);
 315         }
 316         else bc_rand_seedZeroes(r, rng, 1);
 317 }
 318
 319 BcRand bc_rand_int(BcRNG *r) {
 320
 321         BcRNGData *rng = bc_vec_top(&r->v);
 322
 323         if (BC_ERR(BC_RAND_ZERO(rng))) bc_rand_srand(rng);
 324
 325         bc_rand_step(rng);
 326         bc_rand_propagate(r, rng);
 327
 328         return bc_rand_output(rng);
 329 }
 330
 331 BcRand bc_rand_bounded(BcRNG *r, BcRand bound) {
 332
 333         BcRand rand, threshold = (0 - bound) % bound;
 334
 335         do {
 336                 rand = bc_rand_int(r);
 337         } while (rand < threshold);
 338
 339         return rand % bound;
 340 }
 341
 342 void bc_rand_seed(BcRNG *r, ulong state1, ulong state2, ulong inc1, ulong inc2)
 343 {
 344         BcRNGData *rng = bc_vec_top(&r->v);
 345
 346         bc_rand_seedState(&rng->inc, inc1, inc2);
 347         bc_rand_setInc(rng);
 348         bc_rand_setModified(rng);
 349
 350         if (!state1 && !state2) {
 351                 memcpy(&rng->state, &rng->inc, sizeof(BcRandState));
 352                 bc_rand_step(rng);
 353         }
 354         else bc_rand_seedState(&rng->state, state1, state2);
 355
 356         bc_rand_propagate(r, rng);
 357 }
 358
 359 static BcRandState bc_rand_getInc(BcRNGData *r) {
 360
 361         BcRandState res;
 362
 363 #if BC_RAND_BUILTIN
 364         res = r->inc >> 1;
 365 #else // BC_RAND_BUILTIN
 366         res = r->inc;
 367         res.lo >>= 1;
 368         res.lo |= (res.hi & 1) << (BC_LONG_BIT - 1);
 369         res.hi >>= 1;
 370 #endif // BC_RAND_BUILTIN
 371
 372         return res;
 373 }
 374
 375 void bc_rand_getRands(BcRNG *r, BcRand *s1, BcRand *s2, BcRand *i1, BcRand *i2)
 376 {
 377         BcRandState inc;
 378         BcRNGData *rng = bc_vec_top(&r->v);
 379
 380         if (BC_ERR(BC_RAND_ZERO(rng))) bc_rand_srand(rng);
 381
 382         inc = bc_rand_getInc(rng);
 383
 384         *s1 = BC_RAND_TRUNC(rng->state);
 385         *s2 = BC_RAND_CHOP(rng->state);
 386
 387         *i1 = BC_RAND_TRUNC(inc);
 388         *i2 = BC_RAND_CHOP(inc);
 389 }
 390
 391 void bc_rand_push(BcRNG *r) {
 392         BcRNGData rng;
 393         memset(&rng, 0, sizeof(BcRNGData));
 394         if (r->v.len > 0) bc_rand_copy(&rng, bc_vec_top(&r->v));
 395         bc_vec_push(&r->v, &rng);
 396 }
 397
 398 void bc_rand_pop(BcRNG *r, bool reset) {
 399         bc_vec_npop(&r->v, reset ? r->v.len - 1 : 1);
 400 }
 401
 402 void bc_rand_init(BcRNG *r) {
 403         BC_SIG_ASSERT_LOCKED;
 404         bc_vec_init(&r->v, sizeof(BcRNGData), NULL);
 405         bc_rand_push(r);
 406 }
 407
 408 #ifndef NDEBUG
 409 void bc_rand_free(BcRNG *r) {
 410         BC_SIG_ASSERT_LOCKED;
 411         bc_vec_free(&r->v);
 412 }
 413 #endif // NDEBUG
 414
 415 #endif // BC_ENABLE_EXTRA_MATH