prng/xs128p.d

   1 /* Written in 2016 by David Blackman and Sebastiano Vigna (vigna@acm.org)
   2  * To the extent possible under law, the author has dedicated all copyright
   3  * and related and neighboring rights to this software to the public domain
   4  * worldwide. This software is distributed without any warranty.
   5  *
   6  * See <http://creativecommons.org/publicdomain/zero/1.0/>.
   7  */
   8 // D port by Ketmar // Invisible Vector
   9 // http://xoroshiro.di.unimi.it/xoroshiro128plus.c
  10 module iv.prng.xs128p /*is aliced*/;
  11 import iv.alice;
  12
  13 /* This is the successor to xorshift128+. It is the fastest full-period
  14  * generator passing BigCrush without systematic failures, but due to the
  15  * relatively short period it is acceptable only for applications with a
  16  * mild amount of parallelism; otherwise, use a xorshift1024* generator.
  17  *
  18  * Beside passing BigCrush, this generator passes the PractRand test suite
  19  * up to (and included) 16TB, with the exception of binary rank tests,
  20  * which fail due to the lowest bit being an LFSR; all other bits pass all
  21  * tests.
  22  *
  23  * Note that the generator uses a simulated rotate operation, which most C
  24  * compilers will turn into a single instruction. In Java, you can use
  25  * Long.rotateLeft(). In languages that do not make low-level rotation
  26  * instructions accessible xorshift128+ could be faster.
  27  *
  28  * The state must be seeded so that it is not everywhere zero. If you have
  29  * a 64-bit seed, we suggest to seed a splitmix64 generator and use its
  30  * output to fill s. */
  31 struct XS128P {
  32 private:
  33   enum si0 = 0x4d00000000a6829buL;
  34   enum si1 = 0x000029a000000000uL;
  35   ulong[2] s = [si0,si1];//[0x29a,0];
  36
  37 public:
  38 pure nothrow @trusted @nogc:
  39   enum bool isUniformRandom = true;
  40   // tnx to Joseph Rushton Wakeling
  41   enum ulong min = ulong.min;
  42   enum ulong max = ulong.max;
  43
  44   enum bool empty = false;
  45
  46   this (ulong s0, ulong s1=0) { seed(s0, s1); }
  47   this() (auto ref ulong[2] as) { seed(as[]); }
  48
  49   @property ulong front () const { pragma(inline, true); return s.ptr[0]+s.ptr[1]; }
  50
  51   auto save () const { pragma(inline, true); return XS128P(s.ptr[0], s.ptr[1]); }
  52
  53   void popFront () {
  54     immutable ulong s1 = s.ptr[1]^s.ptr[0];
  55     s.ptr[0] = mixin(rol!("s.ptr[0]", 55))^s1^(s1<<14); // a, b
  56     s.ptr[1] = mixin(rol!("s1", 36)); // c
  57   }
  58
  59   void seed (ulong s0, ulong s1=0) { pragma(inline, true); s.ptr[0] = s0; s.ptr[1] = s1; if (!s.ptr[0] && !s.ptr[1]) { s.ptr[0] = si0; s.ptr[1] = si1; } }
  60   void seed() (auto ref ulong[2] as) { pragma(inline, true); s[] = as[]; if (!s.ptr[0] && !s.ptr[1]) { s.ptr[0] = si0; s.ptr[1] = si1; } }
  61
  62   /* This is the jump function for the generator. It is equivalent
  63      to 2^64 calls to next(); it can be used to generate 2^64
  64      non-overlapping subsequences for parallel computations. */
  65   void jump () {
  66     ulong s0 = 0;
  67     ulong s1 = 0;
  68     foreach (ulong jmp; [0xbeac0467eba5facb, 0xd86b048b86aa9922]) {
  69       foreach (immutable b; 0..64) {
  70         if (cast(ulong)(jmp&1)<<b) {
  71           s0 ^= s.ptr[0];
  72           s1 ^= s.ptr[1];
  73         }
  74         popFront();
  75       }
  76     }
  77     s.ptr[0] = s0;
  78     s.ptr[1] = s1;
  79   }
  80
  81 private:
  82   static enum rol(string var, int count) /*if (count > 0 && count < 64)*/ = "("~var~"<<"~count.stringof~")|("~var~">>(64-"~count.stringof~"))";
  83 }
  84
  85
  86 version(test_xs128p) unittest {
  87   static immutable ulong[8] checkValues = [
  88     5548480508102869659uL,
  89     1528641388659518426uL,
  90     5608521991067537321uL,
  91     7624901879973463697uL,
  92     10149332692067980019uL,
  93     17875739269509510643uL,
  94     3809796149818962200uL,
  95     13627210250571023489uL,
  96   ];
  97   {
  98     auto rng = XS128P(0);
  99     foreach (ulong v; checkValues) {
 100       if (v != rng.front) assert(0);
 101       rng.popFront();
 102     }
 103   }
 104   // std.random test
 105   {
 106     import std.random : uniform;
 107     auto rng = XS128P(0);
 108     foreach (immutable _; 0..8) {
 109       import std.stdio;
 110       auto v = uniform!"[)"(0, 4, rng);
 111       writeln(v, "uL");
 112     }
 113   }
 114 }