src/rand.c

   1 /* Pseudo random number generation functions derived from the drand48()
   2  * function obtained from pysam source code.
   3  *
   4  * This functions are used in order to replace the default math.random()
   5  * Lua implementation with something having exactly the same behavior
   6  * across different systems (by default Lua uses libc's rand() that is not
   7  * required to implement a specific PRNG generating the same sequence
   8  * in different systems if seeded with the same integer).
   9  *
  10  * The original code appears to be under the public domain.
  11  * I modified it removing the non needed functions and all the
  12  * 1960-style C coding stuff...
  13  *
  14  * ----------------------------------------------------------------------------
  15  *
  16  * Copyright (c) 2011, Salvatore Sanfilippo <antirez at gmail dot com>
  17  * All rights reserved.
  18  *
  19  * Redistribution and use in source and binary forms, with or without
  20  * modification, are permitted provided that the following conditions are met:
  21  *
  22  *   * Redistributions of source code must retain the above copyright notice,
  23  *     this list of conditions and the following disclaimer.
  24  *   * Redistributions in binary form must reproduce the above copyright
  25  *     notice, this list of conditions and the following disclaimer in the
  26  *     documentation and/or other materials provided with the distribution.
  27  *   * Neither the name of Redis nor the names of its contributors may be used
  28  *     to endorse or promote products derived from this software without
  29  *     specific prior written permission.
  30  *
  31  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  32  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  33  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  34  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  35  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  36  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  37  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  38  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  39  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  40  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  41  * POSSIBILITY OF SUCH DAMAGE.
  42  */
  43
  44 #include <stdint.h>
  45
  46 #define N       16
  47 #define MASK    ((1 << (N - 1)) + (1 << (N - 1)) - 1)
  48 #define LOW(x)  ((unsigned)(x) & MASK)
  49 #define HIGH(x) LOW((x) >> N)
  50 #define MUL(x, y, z)    { int32_t l = (long)(x) * (long)(y); \
  51                 (z)[0] = LOW(l); (z)[1] = HIGH(l); }
  52 #define CARRY(x, y)     ((int32_t)(x) + (long)(y) > MASK)
  53 #define ADDEQU(x, y, z) (z = CARRY(x, (y)), x = LOW(x + (y)))
  54 #define X0      0x330E
  55 #define X1      0xABCD
  56 #define X2      0x1234
  57 #define A0      0xE66D
  58 #define A1      0xDEEC
  59 #define A2      0x5
  60 #define C       0xB
  61 #define SET3(x, x0, x1, x2)     ((x)[0] = (x0), (x)[1] = (x1), (x)[2] = (x2))
  62 #define SETLOW(x, y, n) SET3(x, LOW((y)[n]), LOW((y)[(n)+1]), LOW((y)[(n)+2]))
  63 #define SEED(x0, x1, x2) (SET3(x, x0, x1, x2), SET3(a, A0, A1, A2), c = C)
  64 #define REST(v) for (i = 0; i < 3; i++) { xsubi[i] = x[i]; x[i] = temp[i]; } \
  65                 return (v);
  66 #define HI_BIT  (1L << (2 * N - 1))
  67
  68 static uint32_t x[3] = { X0, X1, X2 }, a[3] = { A0, A1, A2 }, c = C;
  69 static void next();
  70
  71 int32_t redisLrand48() {
  72     next();
  73     return (((int32_t)x[2] << (N - 1)) + (x[1] >> 1));
  74 }
  75
  76 void redisSrand48(int32_t seedval) {
  77     SEED(X0, LOW(seedval), HIGH(seedval));
  78 }
  79
  80 static void next() {
  81     uint32_t p[2], q[2], r[2], carry0, carry1;
  82
  83     MUL(a[0], x[0], p);
  84     ADDEQU(p[0], c, carry0);
  85     ADDEQU(p[1], carry0, carry1);
  86     MUL(a[0], x[1], q);
  87     ADDEQU(p[1], q[0], carry0);
  88     MUL(a[1], x[0], r);
  89     x[2] = LOW(carry0 + carry1 + CARRY(p[1], r[0]) + q[1] + r[1] +
  90             a[0] * x[2] + a[1] * x[1] + a[2] * x[0]);
  91     x[1] = LOW(p[1] + r[0]);
  92     x[0] = LOW(p[0]);
  93 }