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