etc/gentable_XOR/src/msws.h

   1 /**************************************************************************\
   2 *                                                                          *
   3 *  Middle Square Weyl Sequence Random Number Generator                     *
   4 *                                                                          *
   5 *  msws() - returns a 32 bit unsigned int [0,0xffffffff]                   *
   6 *                                                                          *
   7 *  The state vector consists of three 64 bit words:  x, w, and s           *
   8 *                                                                          *
   9 *  x - random output [0,0xffffffff]                                        *
  10 *  w - Weyl sequence (period 2^64)                                         *
  11 *  s - an odd constant                                                     *
  12 *                                                                          *
  13 *  The Weyl sequence is added to the square of x.  The middle is extracted *
  14 *  by shifting right 32 bits:                                              *
  15 *                                                                          *
  16 *  x *= x; x += (w += s); return x = (x>>32) | (x<<32);                    *
  17 *                                                                          *
  18 *  The constant s should be set to a random 64-bit pattern with the upper  *
  19 *  32 bits non-zero and the least significant bit set to 1.  There are     *
  20 *  on the order of 2^63 numbers of this type.  Since the stream length     *
  21 *  is 2^64, this provides about 2^127 random numbers in total.             *
  22 *                                                                          *
  23 *  Note:  This version of the RNG includes an idea proposed by             *
  24 *  Richard P. Brent (creator of the xorgens RNG).  Brent suggested         *
  25 *  adding the Weyl sequence after squaring instead of before squaring.     *
  26 *  This provides a basis for uniformity in the output.                     *
  27 *                                                                          *
  28 *  Copyright (c) 2014-2018 Bernard Widynski                                *
  29 *                                                                          *
  30 *  This code can be used under the terms of the GNU General Public License *
  31 *  as published by the Free Software Foundation, either version 3 of the   *
  32 *  License, or any later version. See the GPL license at URL               *
  33 *  http://www.gnu.org/licenses                                             *
  34 *                                                                          *
  35 \**************************************************************************/
  36
  37 /**************************************************************************\
  38 *                                                                          *
  39 *  Modifications applied by LoRd_MuldeR <mulder2@gmx.de>:                  *
  40 *                                                                          *
  41 *  1. Added function to get random 32-Bit value with upper limit           *
  42 *  2. Added function to get random 64-Bit value                            *
  43 *  3. Added function to get random byte-sequence of arbitrary length       *
  44 *  4. Improved initialization (seeding) function                           *
  45 *  5. Made all functions thread-safe by avoiding the use of static vars    *
  46 *  6. Implemented C++ wrapper class, for convenience                       *
  47 *                                                                          *
  48 \**************************************************************************/
  49
  50 #ifndef _INC_MSWS_H
  51 #define _INC_MSWS_H
  52
  53 #include <stddef.h>
  54 #include <stdint.h>
  55
  56 typedef uint64_t msws_t[3];
  57
  58 inline static uint32_t msws_uint32(msws_t ctx)
  59 {
  60         ctx[0] *= ctx[0]; ctx[0] += (ctx[1] += ctx[2]);
  61         return (uint32_t)(ctx[0] = (ctx[0] >> 32) | (ctx[0] << 32));
  62 }
  63
  64 static inline uint32_t msws_uint32_max(msws_t ctx, const uint32_t max)
  65 {
  66         static const uint32_t DIV[64] =
  67         {
  68                 0xFFFFFFFF, 0xFFFFFFFF, 0x80000000, 0x55555556,
  69                 0x40000000, 0x33333334, 0x2AAAAAAB, 0x24924925,
  70                 0x20000000, 0x1C71C71D, 0x1999999A, 0x1745D175,
  71                 0x15555556, 0x13B13B14, 0x12492493, 0x11111112,
  72                 0x10000000, 0x0F0F0F10, 0x0E38E38F, 0x0D79435F,
  73                 0x0CCCCCCD, 0x0C30C30D, 0x0BA2E8BB, 0x0B21642D,
  74                 0x0AAAAAAB, 0x0A3D70A4, 0x09D89D8A, 0x097B425F,
  75                 0x0924924A, 0x08D3DCB1, 0x08888889, 0x08421085,
  76                 0x08000000, 0x07C1F07D, 0x07878788, 0x07507508,
  77                 0x071C71C8, 0x06EB3E46, 0x06BCA1B0, 0x06906907,
  78                 0x06666667, 0x063E7064, 0x06186187, 0x05F417D1,
  79                 0x05D1745E, 0x05B05B06, 0x0590B217, 0x0572620B,
  80                 0x05555556, 0x0539782A, 0x051EB852, 0x05050506,
  81                 0x04EC4EC5, 0x04D4873F, 0x04BDA130, 0x04A7904B,
  82                 0x04924925, 0x047DC120, 0x0469EE59, 0x0456C798,
  83                 0x04444445, 0x04325C54, 0x04210843, 0x04104105
  84         };
  85         return (max < 64)
  86                 ? (msws_uint32(ctx) / DIV[max])
  87                 : (msws_uint32(ctx) / (UINT32_MAX / max + 1U));
  88 }
  89
  90 inline static uint64_t msws_uint64(msws_t ctx)
  91 {
  92         return (((uint64_t)msws_uint32(ctx)) << 32U) | msws_uint32(ctx);
  93 }
  94
  95 inline static void msws_bytes(msws_t ctx, uint8_t *buffer, size_t len)
  96 {
  97         if (((len & 3U) == 0U) && ((((uintptr_t)buffer) & 3U) == 0U))
  98         {
  99                 uint32_t *ptr = (uint32_t*)buffer;
 100                 for (len >>= 2U; len; --len)
 101                 {
 102                         *ptr++ = msws_uint32(ctx);
 103                 }
 104         }
 105         else
 106         {
 107                 uint32_t tmp;
 108                 for(size_t sft = 0U; len; --len, sft = (sft + 1U) & 3U)
 109                 {
 110                         *buffer++ = (uint8_t)(sft ? (tmp >>= 8U) : (tmp = msws_uint32(ctx)));
 111                 }
 112         }
 113 }
 114
 115 inline static void msws_init(msws_t ctx, const uint32_t seed)
 116 {
 117         ctx[0] = UINT64_C(0); ctx[1] = UINT64_C(0);
 118         ctx[2] = (((uint64_t)seed) << 1U) + 0xB5AD4ECEDA1CE2A9;
 119         for (int i = 0; i < 13; ++i)
 120         {
 121                 volatile uint32_t q = msws_uint32(ctx);
 122         }
 123 }
 124
 125 #endif //_INC_MSWS_H