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 37 // http://www.iro.umontreal.ca/~lecuyer/myftp/papers/lfsr04.pdf 38 // http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt 39 private struct WELLConstants_t(UInt) 40 { 41 alias UIntType = UInt; 42 43 uint wordSize, regSize, rotN; 44 45 uint[3] m; 46 47 string[8] ts; 48 49 bool doTempering = false; 50 UIntType[2] temperingBC; 51 52 enum UIntType[8] a = 53 [0, 54 0xda442d24, 55 0xd3e43ffd, 56 0x8bdcb91e, 57 0x86a9d87e, 58 0xa8c296d1, 59 0x5d6b45cc, 60 0xb729fcec]; 61 62 UIntType M(uint n : 0)(UIntType x) const { return 0; } 63 UIntType M(uint n : 1)(UIntType x) const { return x; } 64 UIntType M(uint n : 2, int t)(UIntType x) const 65 { 66 static if(t >= 0) 67 return x >> t; 68 else 69 return x << (-t); 70 } 71 72 UIntType M(uint n : 3, int t)(UIntType x) const { return x ^ M!(2, t)(x); } 73 UIntType M(uint n : 4, uint a)(UIntType x) const { return (x & 1u) ? ((v >> 1) ^ a) : (v >> 1); } 74 UIntType M(uint n : 5, int t, uint b)(UIntType x) const 75 { 76 return x ^ (M!(2, t)(x) & b); 77 } 78 79 UIntType M(uint n : 6, uint r, uint s, uint t, uint a)(UIntType x) const 80 { 81 immutable ds = ~(1 << (wordSize - 1 - s)), 82 dt = (1 << (wordSize - 1 - t)), 83 rot = (x << r) ^ (x >> (wordSize - r)); 84 85 if(x & dt) 86 return (rot & ds) ^ a; 87 else 88 return rot & ds; 89 } 90 91 92 static UIntType T(alias thisValue, uint n)(UIntType x) 93 if(is(typeof(thisValue) == typeof(this))) 94 { 95 with(thisValue) 96 { 97 return mixin("thisValue." ~ thisValue.ts[n] ~ "(x)"); 98 } 99 } 100 } 101 102 103 private auto _makeStructConstant(S, string expr)() 104 { 105 mixin(`S ret = {` ~ expr ~ `};`); 106 107 return ret; 108 } 109 110 111 enum WELLConstants_t!UIntType[string] 112 WELLConstants(UIntType) = 113 ["512a" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 114 wordSize : 32, regSize : 16, rotN : 0, 115 m : [13, 9, 5], 116 ts : [ 117 "M!(3, -16)", "M!(3, -15)", "M!(3, +11)", "M!(0)", 118 "M!(3, -2)", "M!(3, -18)", "M!(2, -28)", "M!(5, -5, a[1])", 119 ] 120 }), 121 122 "521a" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 123 wordSize : 32, regSize : 17, rotN : 23, 124 m : [13, 11, 10], 125 ts : [ 126 "M!(3, -13)", "M!(3, -15)", "M!(1)", "M!(2, -21)", 127 "M!(3, -13)", "M!(2, 1)", "M!(0)", "M!(3, 11)" 128 ] 129 }), 130 131 "521b" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 132 wordSize : 32, regSize : 17, rotN : 23, 133 m : [11, 10, 7], 134 ts : [ 135 "M!(3, -21)", "M!(3, 6)", "M!(0)", "M!(3, -13)", 136 "M!(3, 13)", "M!(2, -10)", "M!(2, -5)", "M!(3, 13)" 137 ] 138 }), 139 140 "607a" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 141 wordSize : 32, regSize : 19, rotN : 1, 142 m : [16, 15, 14], 143 ts : [ 144 "M!(3, 19)", "M!(3, 11)", "M!(3, -14)", "M!(1)", 145 "M!(3, 18)", "M!(1)", "M!(0)", "M!(3, -15)" 146 ] 147 }), 148 149 "607b" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 150 wordSize : 32, regSize : 25, rotN : 0, 151 m : [16, 8, 13], 152 ts : [ 153 "M!(3, -18)", "M!(3, -14)", "M!(0)", "M!(3, 18)", 154 "M!(3, -24)", "M!(3, 5)", "M!(3, -1)", "M!(0)" 155 ] 156 }), 157 158 "800a" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 159 wordSize : 32, regSize : 25, rotN : 0, 160 m : [14, 18, 17], 161 ts : [ 162 "M!(1)", "M!(3, -15)", "M!(3, 10)", "M!(3, -11)", 163 "M!(3, 16)", "M!(2, 20)", "M!(1)", "M!(3, -28)" 164 ] 165 }), 166 167 "800b" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 168 wordSize : 32, regSize : 25, rotN : 0, 169 m : [9, 4, 22], 170 ts : [ 171 "M!(3, -29)", "M!(2, -14)", "M!(1)", "M!(2, 19)", 172 "M!(1)", "M!(3, 10)", "M!(4, a[2])", "M!(3, -25)" 173 ] 174 }), 175 176 "1024a" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 177 wordSize : 32, regSize : 32, rotN : 0, 178 m : [3, 24, 10], 179 ts : [ 180 "M!(1)", "M!(3, 8)", "M!(3, -19)", "M!(3, -14)", 181 "M!(3, -11)", "M!(3, -7)", "M!(3, -13)", "M!(0)" 182 ] 183 }), 184 185 "1024b" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 186 wordSize : 32, regSize : 32, rotN : 0, 187 m : [22, 25, 26], 188 ts : [ 189 "M!(3, -21)", "M!(3, 17)", "M!(4, a[3])", "M!(3, 15)", 190 "M!(3, -14)", "M!(3, -21)", "M!(1)", "M!(0)" 191 ] 192 }), 193 194 "19937a" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 195 wordSize : 32, regSize : 624, rotN : 31, 196 m : [70, 179, 449], 197 ts : [ 198 "M!(3, -25)", "M!(3, 27)", "M!(2, 9)", "M!(3, 1)", 199 "M!(1)", "M!(3, -9)", "M!(3, -21)", "M!(3, 21)" 200 ] 201 }), 202 203 "19937b" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 204 wordSize : 32, regSize : 624, rotN : 31, 205 m : [203, 613, 123], 206 ts : [ 207 "M!(3, 7)", "M!(1)", "M!(3, 12)", "M!(3, -10)", 208 "M!(3, -19)", "M!(2, -11)", "M!(3, 4)", "M!(3, 10)" 209 ] 210 }), 211 212 "19937c" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 213 wordSize : 32, regSize : 624, rotN : 31, 214 m : [70, 179, 449], 215 ts : [ 216 "M!(3, -25)", "M!(3, 27)", "M!(2, 9)", "M!(3, 1)", 217 "M!(1)", "M!(3, -9)", "M!(3, -21)", "M!(3, 21)" 218 ], 219 doTempering : true, 220 temperingBC : [0xe46e1700U, 0x9b868000U] 221 }), 222 223 "21701a" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 224 wordSize : 32, regSize : 679, rotN : 27, 225 m : [151, 327, 84], 226 ts : [ 227 "M!(1)", "M!(3, -26)", "M!(3, 19)", "M!(0)", 228 "M!(3, 27)", "M!(3, -11)", "M!(6, 15, 27, 10, a[4])", "M!(3, -16)" 229 ] 230 }), 231 232 "23209a" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 233 wordSize : 32, regSize : 726, rotN : 23, 234 m : [667, 43, 462], 235 ts : [ 236 "M!(3, 28)", "M!(1)", "M!(3, 18)", "M!(3, 3)", 237 "M!(3, 21)", "M!(3, -17)", "M!(3, -28)", "M!(3, -1)" 238 ] 239 }), 240 241 "23209b" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 242 wordSize : 32, regSize : 726, rotN : 23, 243 m : [610, 175, 662], 244 ts : [ 245 "M!(4, a[5])", "M!(1)", "M!(6, 15, 15, 30, a[6])", "M!(3, -24)", 246 "M!(3, -26)", "M!(1)", "M!(0)", "M!(3, 16)" 247 ] 248 }), 249 250 "44497a" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 251 wordSize : 32, regSize : 1391, rotN : 15, 252 m : [23, 481, 229], 253 ts : [ 254 "M!(3, -24)", "M!(3, 30)", "M!(3, -10)", "M!(2, -26)", 255 "M!(1)", "M!(3, 20)", "M!(6, 9, 5, 14, a[7])", "M!(1)" 256 ] 257 }), 258 259 "44497b" : _makeStructConstant!(WELLConstants_t!UIntType, q{ 260 wordSize : 32, regSize : 1391, rotN : 15, 261 m : [23, 481, 229], 262 ts : [ 263 "M!(3, -24)", "M!(3, 30)", "M!(3, -10)", "M!(2, -26)", 264 "M!(1)", "M!(3, 20)", "M!(6, 9, 5, 14, a[7])", "M!(1)" 265 ], 266 doTempering : true, 267 temperingBC : [0x93dd1400U, 0xfa118000U] 268 }) 269 ]; 270 271 /** WELL(512a) Random Number Generator. 272 See: http://www.iro.umontreal.ca/~panneton/WELLRNG.html 273 */ 274 alias WELLEngine(string name) = WELLEngine!(uint, name); 275 276 /// ditto 277 struct WELLEngine(UIntType, string name) 278 if( name == "512a" 279 || name == "521a" || name == "521b" 280 || name == "607a" || name == "607b" 281 || name == "800a" || name == "800b" 282 || name == "1024a" || name == "1024b" 283 || name == "19937a" || name == "19937b" || name == "19937c" 284 || name == "21701a" 285 || name == "23209a" || name == "23209b" 286 || name == "44497a" || name == "44497b") 287 { 288 enum Constant = WELLConstants!UIntType[name]; 289 enum size_t _stateSize = Constant.regSize; 290 291 static if(Constant.rotN == 0) 292 enum uint MASKU = 0; 293 else 294 enum uint MASKU = (~0U) >> (Constant.wordSize - Constant.rotN); 295 296 enum uint MASKL = ~MASKU; 297 298 public: 299 /// Mark as Random Number Generator 300 enum isUniformRandom = true; 301 302 /// maximum value 303 enum UIntType max = uint.max; 304 305 /// minimum value 306 enum UIntType min = 0; 307 308 /** 309 */ 310 this(uint value) 311 { 312 seed(value); 313 } 314 315 316 /** 317 */ 318 void seed(uint value) 319 { 320 _state[0] = value; 321 322 // from std.range.XorshiftEngine. 323 foreach(uint i; 1 .. cast(uint)_state.length) 324 _state[i] = cast(UIntType)(1812433253UL * (_state[i-1] ^ (_state[i-1] >> (Constant.wordSize - 2))) + i + 1); 325 326 popFront(); 327 } 328 329 330 /// range primitives 331 void popFront() 332 { 333 auto V0 = &(_state[_stateIdx]), 334 VM1 = &(_state[(_stateIdx + Constant.m[0]) % Constant.regSize]), 335 VM2 = &(_state[(_stateIdx + Constant.m[1]) % Constant.regSize]), 336 VM3 = &(_state[(_stateIdx + Constant.m[2]) % Constant.regSize]), 337 VRm1 = &(_state[(_stateIdx + Constant.regSize - 1) % Constant.regSize]), 338 VRm2 = &(_state[(_stateIdx + Constant.regSize - 2) % Constant.regSize]), 339 newV0 = VRm1, 340 newV1 = V0, 341 newVRm1 = VRm2; 342 343 immutable z0 = (*VRm1 & MASKL) | (*VRm2 & MASKU), 344 z1 = Constant.T!(Constant, 0)(*V0) ^ Constant.T!(Constant, 1)(*VM1), 345 z2 = Constant.T!(Constant, 2)(*VM2) ^ Constant.T!(Constant, 3)(*VM3), 346 z3 = z1 ^ z2; 347 348 *newV1 = z3; 349 *newV0 = Constant.T!(Constant, 4)(z0) ^ Constant.T!(Constant, 5)(z1) 350 ^ Constant.T!(Constant, 6)(z2) ^ Constant.T!(Constant, 7)(z3); 351 352 _stateIdx = (_stateIdx + Constant.regSize - 1) % Constant.regSize; 353 } 354 355 356 /// ditto 357 @property uint front() pure nothrow @safe 358 { 359 static if(Constant.doTempering) 360 { 361 immutable UIntType x = _state[_stateIdx], 362 y = x ^ ((x << 7) & Constant.temperingBC[0]); 363 364 return y ^ ((y << 15) & Constant.temperingBC[1]); 365 } 366 else 367 return _state[_stateIdx]; 368 } 369 370 371 /// ditto 372 enum bool empty = false; 373 374 375 /// ditto 376 @property typeof(this) save() pure nothrow @safe 377 { 378 return this; 379 } 380 381 382 private: 383 uint[_stateSize] _state; 384 size_t _stateIdx; 385 } 386 387 /// 388 unittest 389 { 390 import std.algorithm; 391 import std.range; 392 import std.stdio; 393 394 WELLEngine!"512a" rng; 395 396 static assert(isUniformRNG!(typeof(rng))); 397 static assert(isSeedable!(typeof(rng))); 398 399 rng.seed(100); 400 401 assert(equal(rng.save.take(8), 402 [ 2230636158, 403 1842930638, 404 155680193, 405 1855495099, 406 2311897807, 407 3102313483, 408 3970788677, 409 3720522367,])); 410 411 // save test 412 auto saved1 = rng.save; 413 auto saved2 = rng.save; 414 415 rng.popFrontN(100); 416 assert(equal(saved1.save.take(64), saved2.save.take(64))); 417 418 saved1.popFrontN(100); 419 assert(equal(rng.save.take(64), saved1.save.take(64))); 420 421 assert(rng.front == 1947823519); 422 rng.popFront(); 423 rng.popFrontN(10000); 424 assert(rng.front == 2831551372); 425 } 426 427 unittest 428 { 429 import std.algorithm; 430 import std.range; 431 import std.stdio; 432 433 WELLEngine!"1024a" rng; 434 rng.seed(100); 435 436 static assert(isUniformRNG!(typeof(rng))); 437 static assert(isSeedable!(typeof(rng))); 438 439 assert(equal(rng.save.take(8), 440 [ 1729689691, 441 963076657, 442 888412938, 443 181100396, 444 3310127585, 445 3649309487, 446 2484075420, 447 1423389279,])); 448 449 rng.popFrontN(100); 450 assert(rng.front == 725664384); 451 rng.popFront(); 452 453 rng.popFrontN(10000); 454 assert(rng.front == 3953644315); 455 } 456 457 unittest 458 { 459 import std.algorithm; 460 import std.range; 461 import std.stdio; 462 463 WELLEngine!"19937a" rng; 464 rng.seed(100); 465 466 static assert(isUniformRNG!(typeof(rng))); 467 static assert(isSeedable!(typeof(rng))); 468 469 assert(equal(rng.save.take(8), 470 [ 3859347685, 471 3376854944, 472 2220854319, 473 1533421060, 474 3247527917, 475 1794400208, 476 2014239377, 477 1401918048,])); 478 479 rng.popFrontN(100); 480 assert(rng.front == 2444980538); 481 rng.popFront(); 482 483 rng.popFrontN(10000); 484 assert(rng.front == 490674394); 485 } 486 487 unittest 488 { 489 import std.algorithm; 490 import std.range; 491 import std.stdio; 492 493 WELLEngine!"44497a" rng; 494 rng.seed(100); 495 496 static assert(isUniformRNG!(typeof(rng))); 497 static assert(isSeedable!(typeof(rng))); 498 499 assert(equal(rng.save.take(8), 500 [ 1904938054, 501 1236099671, 502 761528580, 503 261553665, 504 3145325643, 505 603047593, 506 3491142409, 507 496221207,])); 508 509 rng.popFrontN(100); 510 assert(rng.front == 738539296); 511 rng.popFront(); 512 513 rng.popFrontN(10000); 514 assert(rng.front == 2913233053); 515 } 516 517 unittest 518 { 519 import std.algorithm; 520 import std.range; 521 import std.stdio; 522 523 WELLEngine!"44497b" rng; 524 rng.seed(100); 525 526 static assert(isUniformRNG!(typeof(rng))); 527 static assert(isSeedable!(typeof(rng))); 528 529 assert(equal(rng.save.take(8), 530 [1913523270, 531 3946701399, 532 3211002116, 533 1993047553, 534 3283376203, 535 3676328617, 536 425402121, 537 1179532311,])); 538 539 rng.popFrontN(100); 540 assert(rng.front == 1015818016); 541 rng.popFront(); 542 543 rng.popFrontN(10000); 544 assert(rng.front == 229698717); 545 } 546 547 548 auto dist(string name, string boundaries = "[)", T1, T2, Rng)(T1 a, T2 b, ref Rng rng) 549 if((name == "uni" || name == "uniform") && isUniformRNG!Rng) 550 { 551 return uniform!boundaries(a, b, rng); 552 } 553 554 555 auto ref dist(string name, R, B, Rng)(R p, auto ref B t, auto ref B f, ref Rng rng) 556 if((name == "ber" || name == "bernoulli") && is(R : real) && isUniformRNG!Rng) 557 { 558 immutable R r = uniform01!R(rng); 559 if(r < p) 560 return t; 561 else 562 return f; 563 } 564 565 566 R dist(string name, R, Rng)(R lambda, ref Rng rng) 567 if((name == "exp" || name == "exponential") && is(R : real) && isUniformRNG!Rng) 568 { 569 return -log(uniform01!R(rng)) / lambda; 570 } 571 572 573 R dist(string name, R, Rng)(ref Rng rng) 574 if(name == "stdNormal" && is(R : real) && isUniformRNG!Rng) 575 { 576 immutable x = uniform01!R(rng), 577 y = uniform01!R(rng); 578 579 return sqrt(-2 * log(x)) * cos(2 * PI * y); 580 } 581 582 583 R dist(string name, R, Rng)(R mu, R sigma, ref Rng rng) 584 if(name == "normal" && is(R : real) && isUniformRNG!Rng) 585 { 586 return dist!"stdNormal"(rng) * sigma + mu; 587 } 588 589 590 //N dist(string name, N, R, Rng)(N n, R p, ref Rng rng) 591 //if((name == "bin" || name == "binominal") && is(R : real) && is(N : ulong)) 592 //{ 593 594 //} 595 596 597 template dist(string name, Params...) 598 { 599 auto ref dist(T...)(auto ref T args) 600 { 601 return dist!(name, Params)(forward!args, std.random.rndGen); 602 } 603 } 604 605 606 ///** 607 //一様分布 608 //*/ 609 //auto distGenerator(string name, string boundaries = "[)", T1, T2)(T1 a, T2 b) 610 //{ 611 // static struct Generator() 612 // { 613 // alias F = typeof(uniform(T1.init, T2.init, std.random.rndGen)); 614 615 // F gen(alias bs = boundaries, T1, T2, Rng)(T1 a, T2 b, ref Rng rng) 616 // if(isUniformRNG!Rng) 617 // { 618 // return uniform!bs(a, b, rng); 619 // } 620 621 622 // F gen(alias bs = boundaries, T1, T2)(T1 a, T2 b) 623 // if(isUniformRNG!Rng) 624 // { 625 // return gen!bs(a, b, std.random.rndGen); 626 // } 627 628 629 // F gen(alias bs = boundaries, Rng)(ref Rng rng) 630 // if(isUniformRNG!Rng) 631 // { 632 // return gen!bs(_a, _b, rng); 633 // } 634 635 636 // F gen(alias bs = boundaries)() 637 // if(isUniformRNG!Rng) 638 // { 639 // return gen!bs(std.random.rndGen); 640 // } 641 642 643 // F opCall(Rng)(ref Rng rng) 644 // if(isUniformRNG!Rng) 645 // { 646 // return uniform!boundaries(_a, _b, rng); 647 // } 648 649 650 // F opCall() 651 // { 652 // return opCall(std.random.rndGen); 653 // } 654 655 656 // private: 657 // T1 _a; 658 // T2 _b; 659 // } 660 661 662 // return Generator!()(a, b); 663 //} 664 665 666 ///** 667 //ベルヌーイ分布 668 //*/ 669 //auto distGenerator(string name, R, B = bool)(R p, B t = true, B f = false) 670 //if((name == "ber" || name == "bernoulli") && is(R : real)) 671 //{ 672 // static struct Generator() 673 // { 674 // auto ref gen(R, B, Rng)(R p, auto ref B t, auto ref B f, ref Rng rng) 675 // if(is(R : real) && isUniformRNG!Rng) 676 // { 677 // auto r = uniform01!R(rng); 678 // return r < p ? forward!t : forward!f; 679 // } 680 681 682 // auto ref gen(R, B, Rng)(R p, auto ref B t, auto ref B f) 683 // if(is(R : real)) 684 // { 685 // return gen(p, forward!t, forward!f, std.random.rndGen); 686 // } 687 688 689 // auto ref gen(B, Rng)(auto ref B t, auto ref B f, ref Rng rng) 690 // { 691 // return gen(_p, forward!t, forward!f, std.random.rndGen); 692 // } 693 694 695 // B gen(R, Rng)(R p, ref Rng rng) 696 // { 697 // return gen(p, _t, _f, rng); 698 // } 699 700 701 // //B gen(R, Rng) 702 // } 703 //} 704 705 706 //template distRange(string name, TemplateParams...) 707 //{ 708 // auto distRange(RefRng, Params...)(RefRng rng, auto ref Params params) 709 // { 710 // alias D = typeof(distGenerator!(name, TemplateParams)(params)); 711 // alias F = typeof(D.init(rng)); 712 713 714 // static struct DistRange() 715 // { 716 // F front() @property 717 // { 718 // if(_cached) 719 // return _cach; 720 // else{ 721 // _cach = _dGen(_rng); 722 // _cached = true; 723 // return _cach; 724 // } 725 // } 726 727 728 // void popFront() 729 // { 730 // if(_cached) _cached = false; 731 // else{ 732 // _rng.popFront; 733 // } 734 // } 735 736 737 // static if(isInfinite!RefRng) 738 // enum bool empty = false; 739 // else 740 // { 741 // bool empty() @property 742 // { 743 // return _rng.empty; 744 // } 745 // } 746 747 748 // private: 749 // RefRng _rng; 750 // D _dGen; 751 // F _cach; 752 // bool _cached; 753 // } 754 755 // return DistRange!()(rng, distGenerator!(name, Params)(forward!params), F.init, false); 756 // } 757 //} 758 759 //unittest 760 //{ 761 // import carbon.nonametype; 762 763 // WELLEngine!"512a" rng; 764 // rng.seed(100); 765 766 // auto r = rng.scopedRef!"trusted".distRange!("uni", "[)")(0, 10); 767 // static assert(isInputRange!(typeof(r))); 768 // foreach(e; r.take(1000)) 769 // assert(e >= 0 && e < 10); 770 //}