1 // Written in the D programming language. 2 /* 3 NYSL Version 0.9982 4 5 A. This software is "Everyone'sWare". It means: 6 Anybody who has this software can use it as if he/she is 7 the author. 8 9 A-1. Freeware. No fee is required. 10 A-2. You can freely redistribute this software. 11 A-3. You can freely modify this software. And the source 12 may be used in any software with no limitation. 13 A-4. When you release a modified version to public, you 14 must publish it with your name. 15 16 B. The author is not responsible for any kind of damages or loss 17 while using or misusing this software, which is distributed 18 "AS IS". No warranty of any kind is expressed or implied. 19 You use AT YOUR OWN RISK. 20 21 C. Copyrighted to Kazuki KOMATSU 22 23 D. Above three clauses are applied both to source and binary 24 form of this software. 25 */ 26 27 /** 28 このモジュールは、標準ライブラリのstd.randomを強化します。 29 */ 30 31 module carbon.random; 32 33 34 import std.random; 35 36 // http://www.iro.umontreal.ca/~lecuyer/myftp/papers/lfsr04.pdf 37 // http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt 38 private struct WELLConstants_t(UInt) 39 { 40 alias UIntType = UInt; 41 42 uint wordSize, regSize, rotN; 43 44 uint[3] m; 45 46 string[8] ts; 47 48 bool doTempering = false; 49 UIntType[2] temperingBC; 50 51 enum UIntType[8] a = 52 [0, 53 0xda442d24, 54 0xd3e43ffd, 55 0x8bdcb91e, 56 0x86a9d87e, 57 0xa8c296d1, 58 0x5d6b45cc, 59 0xb729fcec]; 60 61 UIntType M(uint n : 0)(UIntType x) const { return 0; } 62 UIntType M(uint n : 1)(UIntType x) const { return x; } 63 UIntType M(uint n : 2, int t)(UIntType x) const 64 { 65 static if(t >= 0) 66 return x >> t; 67 else 68 return x << (-t); 69 } 70 71 UIntType M(uint n : 3, int t)(UIntType x) const { return x ^ M!(2, t)(x); } 72 UIntType M(uint n : 4, uint a)(UIntType x) const { return (x & 1u) ? ((v >> 1) ^ a) : (v >> 1); } 73 UIntType M(uint n : 5, int t, uint b)(UIntType x) const 74 { 75 return x ^ (M!(2, t)(x) & b); 76 } 77 78 UIntType M(uint n : 6, uint r, uint s, uint t, uint a)(UIntType x) const 79 { 80 immutable ds = ~(1 << (wordSize - 1 - s)), 81 dt = (1 << (wordSize - 1 - t)), 82 rot = (x << r) ^ (x >> (wordSize - r)); 83 84 if(x & dt) 85 return (rot & ds) ^ a; 86 else 87 return rot & ds; 88 } 89 90 91 static UIntType T(alias thisValue, uint n)(UIntType x) 92 if(is(typeof(thisValue) == typeof(this))) 93 { 94 with(thisValue) 95 { 96 return mixin("thisValue." ~ thisValue.ts[n] ~ "(x)"); 97 } 98 } 99 } 100 101 102 private auto _makeStructConstant(S, string expr)() 103 { 104 mixin(`S ret = {` ~ expr ~ `};`); 105 106 return ret; 107 } 108 109 110 enum WELLConstants_t!UIntType[string] 111 WELLConstants(UIntType) = 112 ["512a" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 113 wordSize : 32, regSize : 16, rotN : 0, 114 m : [13, 9, 5], 115 ts : [ 116 "M!(3, -16)", "M!(3, -15)", "M!(3, +11)", "M!(0)", 117 "M!(3, -2)", "M!(3, -18)", "M!(2, -28)", "M!(5, -5, a[1])", 118 ] 119 }), 120 121 "521a" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 122 wordSize : 32, regSize : 17, rotN : 23, 123 m : [13, 11, 10], 124 ts : [ 125 "M!(3, -13)", "M!(3, -15)", "M!(1)", "M!(2, -21)", 126 "M!(3, -13)", "M!(2, 1)", "M!(0)", "M!(3, 11)" 127 ] 128 }), 129 130 "521b" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 131 wordSize : 32, regSize : 17, rotN : 23, 132 m : [11, 10, 7], 133 ts : [ 134 "M!(3, -21)", "M!(3, 6)", "M!(0)", "M!(3, -13)", 135 "M!(3, 13)", "M!(2, -10)", "M!(2, -5)", "M!(3, 13)" 136 ] 137 }), 138 139 "607a" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 140 wordSize : 32, regSize : 19, rotN : 1, 141 m : [16, 15, 14], 142 ts : [ 143 "M!(3, 19)", "M!(3, 11)", "M!(3, -14)", "M!(1)", 144 "M!(3, 18)", "M!(1)", "M!(0)", "M!(3, -15)" 145 ] 146 }), 147 148 "607b" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 149 wordSize : 32, regSize : 25, rotN : 0, 150 m : [16, 8, 13], 151 ts : [ 152 "M!(3, -18)", "M!(3, -14)", "M!(0)", "M!(3, 18)", 153 "M!(3, -24)", "M!(3, 5)", "M!(3, -1)", "M!(0)" 154 ] 155 }), 156 157 "800a" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 158 wordSize : 32, regSize : 25, rotN : 0, 159 m : [14, 18, 17], 160 ts : [ 161 "M!(1)", "M!(3, -15)", "M!(3, 10)", "M!(3, -11)", 162 "M!(3, 16)", "M!(2, 20)", "M!(1)", "M!(3, -28)" 163 ] 164 }), 165 166 "800b" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 167 wordSize : 32, regSize : 25, rotN : 0, 168 m : [9, 4, 22], 169 ts : [ 170 "M!(3, -29)", "M!(2, -14)", "M!(1)", "M!(2, 19)", 171 "M!(1)", "M!(3, 10)", "M!(4, a[2])", "M!(3, -25)" 172 ] 173 }), 174 175 "1024a" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 176 wordSize : 32, regSize : 32, rotN : 0, 177 m : [3, 24, 10], 178 ts : [ 179 "M!(1)", "M!(3, 8)", "M!(3, -19)", "M!(3, -14)", 180 "M!(3, -11)", "M!(3, -7)", "M!(3, -13)", "M!(0)" 181 ] 182 }), 183 184 "1024b" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 185 wordSize : 32, regSize : 32, rotN : 0, 186 m : [22, 25, 26], 187 ts : [ 188 "M!(3, -21)", "M!(3, 17)", "M!(4, a[3])", "M!(3, 15)", 189 "M!(3, -14)", "M!(3, -21)", "M!(1)", "M!(0)" 190 ] 191 }), 192 193 "19937a" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 194 wordSize : 32, regSize : 624, rotN : 31, 195 m : [70, 179, 449], 196 ts : [ 197 "M!(3, -25)", "M!(3, 27)", "M!(2, 9)", "M!(3, 1)", 198 "M!(1)", "M!(3, -9)", "M!(3, -21)", "M!(3, 21)" 199 ] 200 }), 201 202 "19937b" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 203 wordSize : 32, regSize : 624, rotN : 31, 204 m : [203, 613, 123], 205 ts : [ 206 "M!(3, 7)", "M!(1)", "M!(3, 12)", "M!(3, -10)", 207 "M!(3, -19)", "M!(2, -11)", "M!(3, 4)", "M!(3, 10)" 208 ] 209 }), 210 211 "19937c" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 212 wordSize : 32, regSize : 624, rotN : 31, 213 m : [70, 179, 449], 214 ts : [ 215 "M!(3, -25)", "M!(3, 27)", "M!(2, 9)", "M!(3, 1)", 216 "M!(1)", "M!(3, -9)", "M!(3, -21)", "M!(3, 21)" 217 ], 218 doTempering : true, 219 temperingBC : [0xe46e1700, 0x9b868000] 220 }), 221 222 "21701a" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 223 wordSize : 32, regSize : 679, rotN : 27, 224 m : [151, 327, 84], 225 ts : [ 226 "M!(1)", "M!(3, -26)", "M!(3, 19)", "M!(0)", 227 "M!(3, 27)", "M!(3, -11)", "M!(6, 15, 27, 10, a[4])", "M!(3, -16)" 228 ] 229 }), 230 231 "23209a" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 232 wordSize : 32, regSize : 726, rotN : 23, 233 m : [667, 43, 462], 234 ts : [ 235 "M!(3, 28)", "M!(1)", "M!(3, 18)", "M!(3, 3)", 236 "M!(3, 21)", "M!(3, -17)", "M!(3, -28)", "M!(3, -1)" 237 ] 238 }), 239 240 "23209b" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 241 wordSize : 32, regSize : 726, rotN : 23, 242 m : [610, 175, 662], 243 ts : [ 244 "M!(4, a[5])", "M!(1)", "M!(6, 15, 15, 30, a[6])", "M!(3, -24)", 245 "M!(3, -26)", "M!(1)", "M!(0)", "M!(3, 16)" 246 ] 247 }), 248 249 "44497a" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 250 wordSize : 32, regSize : 1391, rotN : 15, 251 m : [23, 481, 229], 252 ts : [ 253 "M!(3, -24)", "M!(3, 30)", "M!(3, -10)", "M!(2, -26)", 254 "M!(1)", "M!(3, 20)", "M!(6, 9, 5, 14, a[7])", "M!(1)" 255 ] 256 }), 257 258 "44497b" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 259 wordSize : 32, regSize : 1391, rotN : 15, 260 m : [23, 481, 229], 261 ts : [ 262 "M!(3, -24)", "M!(3, 30)", "M!(3, 10)", "M!(2, -26)", 263 "M!(1)", "M!(3, 20)", "M!(6, 9, 5, 14, a[7])", "M!(1)" 264 ], 265 doTempering : true, 266 temperingBC : [0x93dd1400, 0xfa118000] 267 }) 268 ]; 269 270 /** WELL(512a) Random Number Generator. 271 See: http://www.iro.umontreal.ca/~panneton/WELLRNG.html 272 */ 273 alias WELLEngine(string name) = WELLEngine!(uint, name); 274 275 struct WELLEngine(UIntType, string name) 276 if( name == "512a" 277 || name == "521a" || name == "521b" 278 || name == "607a" || name == "607b" 279 || name == "800a" || name == "800b" 280 || name == "1024a" || name == "1024b" 281 || name == "19937a" || name == "19937b" || name == "19937c" 282 || name == "21701a" 283 || name == "23209a" || name == "23209b" 284 || name == "44497a" || name == "44497b") 285 { 286 enum Constant = WELLConstants!UIntType[name]; 287 enum size_t _stateSize = Constant.regSize; 288 289 enum uint MASKU = (Constant.rotN == 0) ? 0 : ((~0U) >> (Constant.wordSize - Constant.rotN)); 290 enum uint MASKL = ~MASKU; 291 292 public: 293 /// Mark as Random Number Generator 294 enum isUniformRandom = true; 295 296 /// 取り得る最大値 297 enum UIntType max = uint.max; 298 299 /// 取り得る最小値 300 enum UIntType min = 0; 301 302 /** 303 */ 304 this(uint value) 305 { 306 seed(value); 307 } 308 309 310 /** 311 */ 312 void seed(uint value) 313 { 314 _state[0] = value; 315 316 // from std.range.XorshiftEngine. 317 foreach(uint i; 1 .. cast(uint)_state.length) 318 _state[i] = cast(UIntType)(1812433253UL * (_state[i-1] ^ (_state[i-1] >> (Constant.wordSize - 2))) + i + 1); 319 320 popFront(); 321 } 322 323 324 /// range primitives 325 void popFront() 326 { 327 auto V0 = &(_state[_stateIdx]), 328 VM1 = &(_state[(_stateIdx + Constant.m[0]) % Constant.regSize]), 329 VM2 = &(_state[(_stateIdx + Constant.m[1]) % Constant.regSize]), 330 VM3 = &(_state[(_stateIdx + Constant.m[2]) % Constant.regSize]), 331 VRm1 = &(_state[(_stateIdx + Constant.regSize - 1) % Constant.regSize]), 332 VRm2 = &(_state[(_stateIdx + Constant.regSize - 2) % Constant.regSize]), 333 newV0 = VRm1, 334 newV1 = V0, 335 newVRm1 = VRm2; 336 337 immutable z0 = (*VRm1 & MASKL) | (*VRm2 & MASKU), 338 z1 = Constant.T!(Constant, 0)(*V0) ^ Constant.T!(Constant, 1)(*VM1), 339 z2 = Constant.T!(Constant, 2)(*VM2) ^ Constant.T!(Constant, 3)(*VM3); 340 341 *newV1 = z1 ^ z2; 342 *newV0 = Constant.T!(Constant, 4)(z0) ^ Constant.T!(Constant, 5)(z1) 343 ^ Constant.T!(Constant, 6)(z2) ^ Constant.T!(Constant, 7)(*newV1); 344 345 _stateIdx = (_stateIdx + Constant.regSize - 1) % Constant.regSize; 346 } 347 348 349 /// ditto 350 @property uint front() pure nothrow @safe 351 { 352 static if(Constant.doTempering) 353 { 354 immutable x = _state[_stateIdx], 355 y = x ^ ((x << 7) & Constant.temperingBC[0]); 356 357 return y ^ ((y << 15) & Constant.temperingBC[1]); 358 } 359 else 360 return _state[_stateIdx]; 361 } 362 363 364 /// ditto 365 enum bool empty = false; 366 367 368 /// ditto 369 @property typeof(this) save() pure nothrow @safe 370 { 371 return this; 372 } 373 374 375 private: 376 uint[_stateSize] _state; 377 size_t _stateIdx; 378 } 379 380 /// 381 unittest 382 { 383 import std.algorithm; 384 import std.range; 385 import std.stdio; 386 387 WELLEngine!"512a" rng; 388 389 static assert(isUniformRNG!(typeof(rng))); 390 static assert(isSeedable!(typeof(rng))); 391 392 rng.seed(100); 393 394 // http://sc.yutopp.net/entries/5335700643f75e10dc0014c9 395 assert(equal(rng.save.take(8), 396 [ 2230636158, 397 1842930638, 398 155680193, 399 1855495099, 400 2311897807, 401 3102313483, 402 3970788677, 403 3720522367,])); 404 405 // save test 406 auto saved1 = rng.save; 407 auto saved2 = rng.save; 408 409 rng.popFrontN(100); 410 assert(equal(saved1.save.take(64), saved2.save.take(64))); 411 412 saved1.popFrontN(100); 413 assert(equal(rng.save.take(64), saved1.save.take(64))); 414 415 // http://sc.yutopp.net/entries/53361aca43f75e10dc0014fb 416 assert(rng.front == 1947823519); 417 rng.popFront(); 418 rng.popFrontN(10000); 419 assert(rng.front == 2831551372); 420 } 421 422 unittest 423 { 424 import std.algorithm; 425 import std.range; 426 import std.stdio; 427 428 WELLEngine!"1024a" rng; 429 rng.seed(100); 430 431 static assert(isUniformRNG!(typeof(rng))); 432 static assert(isSeedable!(typeof(rng))); 433 434 // http://sc.yutopp.net/entries/53361e0043f75e10dc001514 435 assert(equal(rng.save.take(8), 436 [ 1729689691, 437 963076657, 438 888412938, 439 181100396, 440 3310127585, 441 3649309487, 442 2484075420, 443 1423389279,])); 444 445 rng.popFrontN(100); 446 assert(rng.front == 725664384); 447 rng.popFront(); 448 449 rng.popFrontN(10000); 450 assert(rng.front == 3953644315); 451 } 452 453 unittest 454 { 455 import std.algorithm; 456 import std.range; 457 import std.stdio; 458 459 WELLEngine!"19937a" rng; 460 rng.seed(100); 461 462 static assert(isUniformRNG!(typeof(rng))); 463 static assert(isSeedable!(typeof(rng))); 464 465 // http://sc.yutopp.net/entries/53365f2b43f75e10dc001532 466 assert(equal(rng.save.take(8), 467 [ 3859347685, 468 3376854944, 469 2220854319, 470 1533421060, 471 3247527917, 472 1794400208, 473 2014239377, 474 1401918048,])); 475 476 rng.popFrontN(100); 477 assert(rng.front == 2444980538); 478 rng.popFront(); 479 480 rng.popFrontN(10000); 481 assert(rng.front == 490674394); 482 } 483 484 unittest 485 { 486 import std.algorithm; 487 import std.range; 488 import std.stdio; 489 490 WELLEngine!"44497a" rng; 491 rng.seed(100); 492 493 static assert(isUniformRNG!(typeof(rng))); 494 static assert(isSeedable!(typeof(rng))); 495 496 // http://sc.yutopp.net/entries/533662e443f75e10dc00154b 497 assert(equal(rng.save.take(8), 498 [ 1904938054, 499 1236099671, 500 761528580, 501 261553665, 502 3145325643, 503 603047593, 504 3491142409, 505 496221207,])); 506 507 rng.popFrontN(100); 508 assert(rng.front == 738539296); 509 rng.popFront(); 510 511 rng.popFrontN(10000); 512 assert(rng.front == 2913233053); 513 }