prng/bj.d

   1 /* Invisible Vector Library
   2  * coded by Ketmar // Invisible Vector <ketmar@ketmar.no-ip.org>
   3  * Understanding is not required. Only obedience.
   4  *
   5  * This program is free software: you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation, either version 3 of the License, or
   8  * (at your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful,
  11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13  * GNU General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
  17  */
  18 module iv.prng.bj /*is aliced*/;
  19 import iv.alice;
  20
  21
  22 // ////////////////////////////////////////////////////////////////////////// //
  23 // by Bob Jenkins
  24 // public domain
  25 // http://burtleburtle.net/bob/rand/smallprng.html
  26 struct BJRng {
  27 private:
  28   // seeded with 0xdeadf00du
  29   uint a = 0xe5595c3bu;
  30   uint b = 0xe60c3611u;
  31   uint c = 0x1ceca1b1u;
  32   uint d = 0x32744417u;
  33
  34 nothrow @nogc:
  35 public:
  36   void randomize () @trusted {
  37     version(Windows) {
  38       import win32.windef, win32.winbase;
  39       uint s0 = xyzzyPRNGHashU32(cast(uint)GetCurrentProcessId());
  40       uint s1 = xyzzyPRNGHashU32(cast(uint)GetTickCount());
  41       seed(s0^s1);
  42     } else {
  43       // assume POSIX
  44       import core.sys.posix.fcntl;
  45       import core.sys.posix.unistd;
  46       uint s0 = 0xdeadf00du;
  47       int fd = open("/dev/urandom", O_RDONLY);
  48       if (fd >= 0) {
  49         read(fd, &s0, s0.sizeof);
  50         close(fd);
  51       }
  52       seed(s0);
  53     }
  54   }
  55
  56 pure @safe:
  57   enum bool isUniformRandom = true;
  58   enum uint min = uint.min;
  59   enum uint max = uint.max;
  60
  61   enum bool empty = false;
  62   @property uint front () const { pragma(inline, true); return d; }
  63   alias popFront = next;
  64   @property auto save () const { BJRng res = void; res.a = this.a; res.b = this.b; res.c = this.c; res.d = this.d; return res; }
  65
  66   this (uint aseed) { pragma(inline, true); seed(aseed); }
  67
  68   void seed (uint seed) {
  69     a = 0xf1ea5eed;
  70     b = c = d = seed;
  71     foreach (immutable _; 0..20) next;
  72   }
  73
  74   @property uint next () {
  75     enum ROT(string var, string cnt) = `cast(uint)(((`~var~`)<<(`~cnt~`))|((`~var~`)>>(32-(`~cnt~`))))`;
  76     uint e;
  77     /* original:
  78       e = a-BJPRNG_ROT(b, 27);
  79       a = b^BJPRNG_ROT(c, 17);
  80       b = c+d;
  81       c = d+e;
  82       d = e+a;
  83     */
  84     /* better, but slower at least in idiotic m$vc */
  85     e = a-mixin(ROT!("b", "23"));
  86     a = b^mixin(ROT!("c", "16"));
  87     b = c+mixin(ROT!("d", "11"));
  88     c = d+e;
  89     d = e+a;
  90     return d;
  91   }
  92 }
  93
  94
  95 version(test_bj) unittest {
  96   static immutable uint[8] checkValues = [
  97     3079471771u,
  98     2798213162u,
  99     3360187034u,
 100     3739077647u,
 101     1276142743u,
 102     771570220u,
 103     1864333648u,
 104     1915806440u,
 105   ];
 106   {
 107     auto rng = BJRng(0);
 108     foreach (ulong v; checkValues) {
 109       if (v != rng.front) assert(0);
 110       //import std.stdio; writeln(rng.front, "u");
 111       rng.popFront();
 112     }
 113   }
 114   // std.random test
 115   {
 116     import std.random : uniform;
 117     auto rng = BJRng(0);
 118     foreach (immutable _; 0..8) {
 119       import std.stdio;
 120       auto v = uniform!"[)"(0, 4, rng);
 121       writeln(v, "uL");
 122     }
 123   }
 124 }
 125
 126
 127 // ////////////////////////////////////////////////////////////////////////// //
 128 // One of the by George Marsaglia's prng generators.
 129 // It is based on George Marsaglia's MWC (multiply with carry) generator.
 130 // Although it is very simple, it passes Marsaglia's DIEHARD series of random
 131 // number generator tests.
 132 struct GMRngSeed64 {
 133 private:
 134   // These values are not magical, just the default values Marsaglia used.
 135   // Any pair of unsigned integers should be fine.
 136   enum uint DefaultW = 521288629u;
 137   enum uint DefaultZ = 362436069u;
 138
 139   uint w = DefaultW;
 140   uint z = DefaultZ;
 141   uint lastnum = (DefaultZ<<16)+DefaultW;
 142
 143 nothrow @nogc:
 144 public:
 145   void randomize () @trusted {
 146     version(Windows) {
 147       import win32.windef, win32.winbase;
 148       w = xyzzyPRNGHashU32(cast(uint)GetCurrentProcessId());
 149       z = xyzzyPRNGHashU32(cast(uint)GetTickCount());
 150     } else {
 151       // assume POSIX
 152       import core.sys.posix.fcntl;
 153       import core.sys.posix.unistd;
 154       w = DefaultW;
 155       z = DefaultZ;
 156       int fd = open("/dev/urandom", O_RDONLY);
 157       if (fd >= 0) {
 158         read(fd, &w, w.sizeof);
 159         read(fd, &z, z.sizeof);
 160         close(fd);
 161       }
 162     }
 163   }
 164
 165 pure @safe:
 166   enum bool isUniformRandom = true;
 167   enum uint min = uint.min;
 168   enum uint max = uint.max;
 169
 170   enum bool empty = false;
 171   @property uint front () const { pragma(inline, true); return lastnum; }
 172   alias popFront = next;
 173   @property auto save () const { GMRngSeed64 res = void; res.w = this.w; res.z = this.z; res.lastnum = this.lastnum; return res; }
 174
 175   this (ulong aseed) { pragma(inline, true); seed(aseed); }
 176
 177   void seed (ulong seed) {
 178     z = cast(uint)(seed>>32);
 179     w = cast(uint)(seed&0xffffffffu);
 180     if (w == 0) w = DefaultW;
 181     if (z == 0) z = DefaultZ;
 182     lastnum = (z<<16)+w;
 183   }
 184
 185   @property uint next () {
 186     if (w == 0) w = DefaultW;
 187     if (z == 0) z = DefaultZ;
 188     z = 36969*(z&0xffff)+(z>>16);
 189     w = 18000*(w&0xffff)+(w>>16);
 190     return (lastnum = (z<<16)+w);
 191   }
 192 }
 193
 194 version(test_gms64) unittest {
 195   static immutable uint[8] checkValues = [
 196     1962359733u,
 197     820856226u,
 198     2331188998u,
 199     4033440000u,
 200     3169966213u,
 201     2572821606u,
 202     100826968u,
 203     1697244543u,
 204   ];
 205   {
 206     auto rng = GMRngSeed64(0);
 207     foreach (ulong v; checkValues) {
 208       if (v != rng.front) assert(0);
 209       //import std.stdio; writeln(rng.front, "u");
 210       rng.popFront();
 211     }
 212   }
 213   // std.random test
 214   {
 215     import std.random : uniform;
 216     auto rng = GMRngSeed64(0);
 217     foreach (immutable _; 0..8) {
 218       import std.stdio;
 219       auto v = uniform!"[)"(0, 4, rng);
 220       writeln(v, "uL");
 221     }
 222   }
 223 }
 224
 225
 226 // ////////////////////////////////////////////////////////////////////////// //
 227 // one of the by George Marsaglia's prng generators, period about 2^160. fast.
 228 struct GMRng {
 229 private:
 230   // seeded with 0xdeadf00du
 231   uint x = 0xdeadf00du;
 232   uint y = 0xe3324aa1u;
 233   uint z = 0x7ed1c277u;
 234   uint w = 0x89574524u;
 235   uint v = 0x359c34a7u;
 236   uint lastnum = 0x4d00381eu; // almost arbitrary
 237
 238 nothrow @nogc:
 239 public:
 240   void randomize () @trusted {
 241     version(Windows) {
 242       import win32.windef, win32.winbase;
 243       uint[5] s0;
 244       s0[0] = xyzzyPRNGHashU32(cast(uint)GetCurrentProcessId());
 245       s0[1] = xyzzyPRNGHashU32(cast(uint)GetTickCount());
 246       seed(s0[]);
 247     } else {
 248       // assume POSIX
 249       import core.sys.posix.fcntl;
 250       import core.sys.posix.unistd;
 251       uint[5] s0;
 252       int fd = open("/dev/urandom", O_RDONLY);
 253       if (fd >= 0) {
 254         read(fd, s0.ptr, s0.sizeof);
 255         close(fd);
 256       }
 257       seed(s0[]);
 258     }
 259   }
 260
 261 pure @safe:
 262   enum bool isUniformRandom = true;
 263   enum uint min = uint.min;
 264   enum uint max = uint.max;
 265
 266   enum bool empty = false;
 267   @property uint front () const { pragma(inline, true); return lastnum; }
 268   alias popFront = next;
 269   @property auto save () const { GMRng res = void; res.x = this.x; res.y = this.y; res.z = this.z; res.w = this.w; res.v = this.v; res.lastnum = this.lastnum; return res; }
 270
 271   this (uint aseed) { pragma(inline, true); seed(aseed); }
 272
 273   void seed (uint seed) {
 274     x = (seed ? seed : 0xdeadf00du);
 275     y = xyzzyPRNGHashU32(x+1);
 276     z = xyzzyPRNGHashU32(y+1);
 277     w = xyzzyPRNGHashU32(z+1);
 278     v = xyzzyPRNGHashU32(w+1);
 279     next;
 280   }
 281
 282   void seed (uint s0, uint s1, uint s2, uint s3, uint s4) {
 283     x = s0;
 284     y = s1;
 285     z = s2;
 286     w = s3;
 287     v = s4;
 288     if (x == 0) x = 0xdeadf00du;
 289     if (y == 0) y = xyzzyPRNGHashU32(x+1);
 290     if (z == 0) z = xyzzyPRNGHashU32(y+1);
 291     if (w == 0) w = xyzzyPRNGHashU32(z+1);
 292     if (v == 0) v = xyzzyPRNGHashU32(w+1);
 293     next;
 294   }
 295
 296   void seed (in uint[] seed) {
 297     x = (seed.length > 0 ? seed[0] : 0xdeadf00du);
 298     y = (seed.length > 1 ? seed[1] : xyzzyPRNGHashU32(x+1));
 299     z = (seed.length > 2 ? seed[2] : xyzzyPRNGHashU32(y+1));
 300     w = (seed.length > 3 ? seed[3] : xyzzyPRNGHashU32(z+1));
 301     v = (seed.length > 4 ? seed[4] : xyzzyPRNGHashU32(w+1));
 302     if (x == 0) x = 0xdeadf00du;
 303     if (y == 0) y = xyzzyPRNGHashU32(x+1);
 304     if (z == 0) z = xyzzyPRNGHashU32(y+1);
 305     if (w == 0) w = xyzzyPRNGHashU32(z+1);
 306     if (v == 0) v = xyzzyPRNGHashU32(w+1);
 307     next;
 308   }
 309
 310   @property uint next () {
 311     uint t;
 312     t = (x^(x>>7));
 313     x = y;
 314     y = z;
 315     z = w;
 316     w = v;
 317     v = (v^(v<<6))^(t^(t<<13));
 318     return (lastnum = (y+y+1)*v);
 319   }
 320 }
 321
 322 version(test_gms) unittest {
 323   static immutable uint[8] checkValues = [
 324     2329293526u,
 325     2821973934u,
 326     2640992451u,
 327     1004589151u,
 328     3902251129u,
 329     2922888142u,
 330     3947715136u,
 331     3516368807u,
 332   ];
 333   {
 334     auto rng = GMRng(0);
 335     foreach (ulong v; checkValues) {
 336       if (v != rng.front) assert(0);
 337       //import std.stdio; writeln(rng.front, "u");
 338       rng.popFront();
 339     }
 340   }
 341   // std.random test
 342   {
 343     import std.random : uniform;
 344     auto rng = GMRng(0);
 345     foreach (immutable _; 0..8) {
 346       import std.stdio;
 347       auto v = uniform!"[)"(0, 4, rng);
 348       writeln(v, "uL");
 349     }
 350   }
 351 }
 352
 353
 354 // ////////////////////////////////////////////////////////////////////////// //
 355 version(test_prng) unittest {
 356   import std.range;
 357   template checkRng(T) {
 358     static assert(isInfinite!T, T.stringof~" is not infinite range");
 359     static assert(isInputRange!T, T.stringof~" is not inpute range");
 360     static assert(isForwardRange!T, T.stringof~" is not forward range");
 361     enum checkRng = true;
 362   }
 363   static assert(checkRng!BJRng);
 364   static assert(checkRng!GMRngSeed64);
 365   static assert(checkRng!GMRng);
 366 }
 367
 368
 369 // ////////////////////////////////////////////////////////////////////////// //
 370 private uint xyzzyPRNGHashU32() (uint a) {
 371   a -= (a<<6);
 372   a ^= (a>>17);
 373   a -= (a<<9);
 374   a ^= (a<<4);
 375   a -= (a<<3);
 376   a ^= (a<<10);
 377   a ^= (a>>15);
 378   return a;
 379 }