1 /++ 2 + Permuted Congruential Generator (PCG) 3 + 4 + Implemented as per the C++ version of PCG, $(HTTP _pcg-random.org). 5 + 6 + Paper available $(HTTP _pcg-random.org/paper.html) 7 + 8 + Author: Melissa O'Neill (C++). D translation Nicholas Wilson. 9 + 10 + PCG Random Number Generation for C++ 11 + 12 + Copyright 2014 Melissa O'Neill <oneill@pcg-random.org> 13 + 14 + Licensed under the Apache License, Version 2.0 (the "License"); 15 + you may not use this file except in compliance with the License. 16 + You may obtain a copy of the License at 17 + 18 + $(HTTP www.apache.org/licenses/LICENSE-2.0) 19 + 20 + Unless required by applicable law or agreed to in writing, software 21 + distributed under the License is distributed on an "AS IS" BASIS, 22 + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 23 + See the License for the specific language governing permissions and 24 + limitations under the License. 25 + 26 + For additional information about the PCG random number generation scheme, 27 + including its license and other licensing options, visit 28 + 29 + $(HTTP _pcg-random.org) 30 +/ 31 module mir.random.engine.pcg; 32 33 import mir.random.engine; 34 import std.traits : ReturnType, TemplateArgsOf; 35 36 @safe: 37 nothrow: 38 @nogc: 39 40 /// 41 @nogc nothrow pure @safe version(mir_random_test) unittest 42 { 43 import mir.random /+: isSaturatedRandomEngine, rand+/; 44 import mir.random.engine.pcg : pcg32; 45 import mir.math.common: fabs; 46 47 static assert(isRandomEngine!pcg32); 48 static assert(isSaturatedRandomEngine!pcg32); 49 auto gen = pcg32(1234u);//Seed with constant. 50 assert(gen.rand!double.fabs == 0x1.e5fe29e2942a8p-1);//Generate number from 0 inclusive to 1 exclusive. 51 assert(gen.rand!ulong == 13957536084079755083UL); 52 } 53 54 /// Select the above mixin templates. 55 enum stream_t 56 { 57 /// Increment is cast(size_t) &RNG. 58 unique, 59 /// Increment is 0. 60 none, 61 /// Increment is the default increment. 62 oneseq, 63 /// Increment is runtime setable and defaults to the default (same as oneseq) 64 specific, 65 } 66 67 /// 32-bit output PCGs with 64 bits of state. 68 alias pcg32 = PermutedCongruentialEngine!(xsh_rr!(uint,ulong),stream_t.specific,true); 69 /// ditto 70 alias pcg32_unique = PermutedCongruentialEngine!(xsh_rr!(uint,ulong),stream_t.unique,true); 71 /// ditto 72 alias pcg32_oneseq = PermutedCongruentialEngine!(xsh_rr!(uint,ulong),stream_t.oneseq,true); 73 /// ditto 74 alias pcg32_fast = PermutedCongruentialEngine!(xsh_rr!(uint,ulong),stream_t.none,true); 75 76 static if (__traits(compiles, ucent.max)) 77 { 78 /// 64-bit output PCGs with 128 bits of state. Requires `ucent` type. 79 alias pcg64 = PermutedCongruentialEngine!(xsh_rr!(ulong,ucent),stream_t.specific,true); 80 /// 81 alias pcg64_unique = PermutedCongruentialEngine!(xsh_rr!(ulong,ucent),stream_t.unique,true); 82 /// 83 alias pcg64_oneseq = PermutedCongruentialEngine!(xsh_rr!(ulong,ucent),stream_t.oneseq,true); 84 /// 85 alias pcg64_fast = PermutedCongruentialEngine!(xsh_rr!(ulong,ucent),stream_t.none,true); 86 } 87 88 /// PCGs with n bits output and n bits of state. 89 alias pcg8_once_insecure = PermutedCongruentialEngine!(rxs_m_xs_forward!(ubyte ,ubyte ),stream_t.specific,true); 90 /// ditto 91 alias pcg16_once_insecure = PermutedCongruentialEngine!(rxs_m_xs_forward!(ushort,ushort),stream_t.specific,true); 92 /// ditto 93 alias pcg32_once_insecure = PermutedCongruentialEngine!(rxs_m_xs_forward!(uint ,uint ),stream_t.specific,true); 94 /// ditto 95 alias pcg64_once_insecure = PermutedCongruentialEngine!(rxs_m_xs_forward!(ulong,ulong ),stream_t.specific,true); 96 //alias pcg128_once_insecure = PermutedCongruentialEngine!(rxs_m_xs_forward!(ucent,ucent,stream_t.specific,true); 97 98 /// As above but the increment is not dynamically setable. 99 alias pcg8_oneseq_once_insecure = PermutedCongruentialEngine!(rxs_m_xs_forward!(ubyte ,ubyte ),stream_t.oneseq,true); 100 /// ditto 101 alias pcg16_oneseq_once_insecure = PermutedCongruentialEngine!(rxs_m_xs_forward!(ushort,ushort),stream_t.oneseq,true); 102 /// ditto 103 alias pcg32_oneseq_once_insecure = PermutedCongruentialEngine!(rxs_m_xs_forward!(uint ,uint ),stream_t.oneseq,true); 104 /// ditto 105 alias pcg64_oneseq_once_insecure = PermutedCongruentialEngine!(rxs_m_xs_forward!(ulong,ulong ),stream_t.oneseq,true); 106 /// ditto 107 /// Requires `ucent` type. 108 static if (__traits(compiles, ucent.max)) 109 alias pcg128_oneseq_once_insecure = PermutedCongruentialEngine!(rxs_m_xs_forward!(ucent,ucent ),stream_t.specific,true); 110 111 /++ 112 + The PermutedCongruentialEngine: 113 + Params: 114 + output = should be one of the above functions. 115 + Controls the output permutation of the state. 116 + streamType = one of unique, none, oneseq, specific. 117 + Controls the Increment of the LCG portion of the PCG. 118 + output_previous = if true then the pre-advance version (increasing instruction-level parallelism) 119 + if false then use the post-advance version (reducing register pressure) 120 + mult_ = optionally set the multiplier for the LCG. 121 +/ 122 struct PermutedCongruentialEngine(alias output, // Output function 123 stream_t streamType, // The stream type 124 bool output_previous, 125 mult_...) if (mult_.length <= 1) 126 { 127 /// 128 enum isRandomEngine = true; 129 130 /// 131 alias Uint = TemplateArgsOf!output[1]; 132 133 static if (mult_.length == 0) 134 enum mult = default_multiplier!Uint; 135 else 136 { 137 static assert(is(typeof(mult_[0]) == Uint), 138 "The specified multiplier must be the state type of the output function"); 139 enum mult = mult_[0]; 140 } 141 142 @disable this(this); 143 @disable this(); 144 static if (streamType == stream_t.none) 145 mixin no_stream!Uint; 146 else static if (streamType == stream_t.unique) 147 mixin unique_stream!Uint; 148 else static if (streamType == stream_t.specific) 149 mixin specific_stream!Uint; 150 else static if (streamType == stream_t.oneseq) 151 mixin oneseq_stream!Uint; 152 else 153 static assert(0); 154 155 /// 156 Uint state; 157 158 /// 159 enum period_pow2 = Uint.sizeof*8 - 2*is_mcg; 160 161 /// 162 enum max = (ReturnType!output).max; 163 164 private: 165 166 static if (__traits(compiles, { enum e = mult + increment; })) 167 { 168 static Uint bump()(Uint state_) 169 { 170 return cast(Uint)(state_ * mult + increment); 171 } 172 } 173 else 174 { 175 Uint bump()(Uint state_) 176 { 177 return cast(Uint)(state_ * mult + increment); 178 } 179 } 180 181 Uint base_generate()() 182 { 183 return state = bump(state); 184 } 185 186 Uint base_generate0()() 187 { 188 Uint old_state = state; 189 state = bump(state); 190 return old_state; 191 } 192 193 public: 194 static if (can_specify_stream) 195 /// 196 this()(Uint seed, Uint stream_ = default_increment_unset_stream!Uint) 197 { 198 state = bump(cast(Uint)(seed + increment)); 199 set_stream(stream_); 200 } 201 else 202 /// 203 this()(Uint seed) 204 { 205 static if (is_mcg) 206 state = seed | 3u; 207 else 208 state = bump(cast(Uint)(seed + increment)); 209 } 210 211 /// 212 ReturnType!output opCall()() 213 { 214 static if(output_previous) 215 return output(base_generate0()); 216 else 217 return output(base_generate()); 218 } 219 220 /++ 221 Skip forward in the random sequence in $(BIGOH log n) time. 222 Even though delta is an unsigned integer, we can pass a 223 signed integer to go backwards, it just goes "the long way round". 224 +/ 225 void skip()(Uint delta) 226 { 227 // The method used here is based on Brown, "Random Number Generation 228 // with Arbitrary Stride,", Transactions of the American Nuclear 229 // Society (Nov. 1994). The algorithm is very similar to fast 230 // exponentiation. 231 // 232 // Even though delta is an unsigned integer, we can pass a 233 // signed integer to go backwards, it just goes "the long way round". 234 235 Uint acc_mult = 1, acc_plus = 0; 236 Uint cur_plus = increment, cur_mult = mult; 237 while (delta > 0) 238 { 239 if (delta & 1) 240 { 241 acc_mult *= cur_mult; 242 acc_plus = cast(Uint)(acc_plus * cur_mult + cur_plus); 243 } 244 cur_plus *= cur_mult + 1; 245 cur_mult *= cur_mult; 246 delta >>= 1; 247 } 248 state = cast(Uint)(acc_mult*state + acc_plus); 249 } 250 251 static if (output_previous) 252 { 253 /++ 254 Compatibility with $(LINK2 https://dlang.org/phobos/std_random.html#.isUniformRNG, 255 Phobos library methods). Presents this RNG as an InputRange. 256 Only available if `output_previous == true`. 257 258 The reason that this is enabled when `output_previous == true` is because 259 `front` can be implemented without additional cost. 260 261 `save` is implemented if `streamType` is not [stream_t.unique]. 262 263 This struct disables its default copy constructor and so will only work with 264 Phobos functions that "do the right thing" and take RNGs by reference and 265 do not accidentally make implicit copies. 266 +/ 267 enum bool isUniformRandom = true; 268 /// ditto 269 enum ReturnType!output min = (ReturnType!output).min; 270 /// ditto 271 enum bool empty = false; 272 /// ditto 273 @property ReturnType!output front()() const { return output(state); } 274 /// ditto 275 void popFront()() { state = bump(state); } 276 /// ditto 277 void seed()(Uint seed) { this.__ctor(seed); } 278 /// ditto 279 static if (streamType != stream_t.unique) 280 @property typeof(this) save()() const 281 { 282 typeof(return) result = void; 283 foreach (i, e; this.tupleof) 284 result.tupleof[i] = e; 285 return result; 286 } 287 } 288 } 289 290 @nogc nothrow pure @safe version(mir_random_test) unittest 291 { 292 //Test that the default generators (all having output_previous==true) 293 //can be used as Phobos-style randoms. 294 import std.meta: AliasSeq; 295 import std.random: isSeedable, isPhobosUniformRNG = isUniformRNG; 296 import std.range.primitives: isForwardRange; 297 foreach(RNG; AliasSeq!(pcg32, pcg32_oneseq, pcg32_fast, 298 pcg32_once_insecure, pcg32_oneseq_once_insecure, 299 pcg64_once_insecure, pcg64_oneseq_once_insecure)) 300 { 301 static assert(isPhobosUniformRNG!(RNG, typeof(RNG.max))); 302 static assert(isSeedable!(RNG, RNG.Uint)); 303 static assert(isForwardRange!RNG); 304 auto gen1 = RNG(1); 305 auto gen2 = RNG(2); 306 auto gen3 = gen1.save; 307 gen2.seed(1); 308 assert(gen1 == gen2); 309 immutable a = gen1.front; 310 gen1.popFront(); 311 assert(a == gen2()); 312 assert(gen1.front == gen2()); 313 assert(a == gen3()); 314 assert(gen1.front == gen3()); 315 } 316 317 foreach(RNG; AliasSeq!(pcg32_unique)) 318 { 319 static assert(isPhobosUniformRNG!(RNG, typeof(RNG.max))); 320 static assert(isSeedable!(RNG, RNG.Uint)); 321 } 322 } 323 324 // Default multiplier to use for the LCG portion of the PCG 325 private template default_multiplier(Uint) 326 { 327 static if (is(Uint == ubyte)) 328 enum ubyte default_multiplier = 141u; 329 else static if (is(Uint == ushort)) 330 enum ushort default_multiplier = 12829u; 331 else static if (is(Uint == uint)) 332 enum uint default_multiplier = 747796405u; 333 else static if (is(Uint == ulong)) 334 enum ulong default_multiplier = 6364136223846793005u; 335 else static if (is(ucent) && is(Uint == ucent)) 336 mixin("enum ucent default_multiplier = 0x2360ED051FC65DA44385DF649FCCF645;"); 337 else 338 static assert(0); 339 } 340 341 // Default increment to use for the LCG portion of the PCG 342 private template default_increment(Uint) 343 { 344 static if (is(Uint == ubyte)) 345 enum ubyte default_increment = 77u; 346 else static if (is(Uint == ushort)) 347 enum ushort default_increment = 47989u; 348 else static if (is(Uint == uint)) 349 enum uint default_increment = 2891336453u; 350 else static if (is(Uint == ulong)) 351 enum ulong default_increment = 1442695040888963407u; 352 else static if (is(ucent) && is(Uint == ucent)) 353 mixin("enum ucent default_increment = 0x5851F42D4C957F2D14057B7EF767814F;"); 354 else 355 static assert(0); 356 } 357 358 private template mcg_multiplier(Uint) 359 { 360 static if (is(Uint == ubyte)) 361 enum ubyte mcg_multiplier = 217u; 362 else static if (is(Uint == ushort)) 363 enum ushort mcg_multiplier = 62169u; 364 else static if (is(Uint == uint)) 365 enum uint mcg_multiplier = 277803737u; 366 else static if (is(Uint == ulong)) 367 enum ulong mcg_multiplier = 12605985483714917081u; 368 else static if (is(ucent) && is(Uint == ucent)) 369 mixin("enum ucent mcg_multiplier = 0x6BC8F622C397699CAEF17502108EF2D9;"); 370 else 371 static assert(0); 372 } 373 374 private template mcg_unmultiplier(Uint) 375 { 376 static if (is(Uint == ubyte)) 377 enum ubyte mcg_unmultiplier = 105u; 378 else static if (is(Uint == ushort)) 379 enum ushort mcg_unmultiplier = 28009u; 380 else static if (is(Uint == uint)) 381 enum uint mcg_unmultiplier = 2897767785u; 382 else static if (is(Uint == ulong)) 383 enum ulong mcg_unmultiplier = 15009553638781119849u; 384 else static if (is(ucent) && is(Uint == ucent)) 385 mixin("enum ucent mcg_unmultiplier = 0xC827645E182BC965D04CA582ACB86D69;"); 386 else 387 static assert(0); 388 } 389 390 private template default_increment_unset_stream(Uint) 391 { 392 enum default_increment_unset_stream = (default_increment!Uint & ~1) >> 1; 393 } 394 395 /// Increment for LCG portion of the PCG is the address of the RNG 396 mixin template unique_stream(Uint) 397 { 398 /// 399 enum is_mcg = false; 400 /// 401 @property Uint increment()() const @trusted 402 { 403 return cast(Uint) (&this) | 1; 404 } 405 /// 406 Uint stream()() 407 { 408 return increment >> 1; 409 } 410 /// 411 enum can_specify_stream = false; 412 /// 413 enum size_t streams_pow2 = Uint.sizeof < size_t.sizeof ? Uint.sizeof : size_t.sizeof - 1u; 414 } 415 416 @nogc nothrow pure @system unittest 417 { 418 pcg32_unique gen = pcg32_unique(1); 419 void* address = &gen; 420 assert(gen.increment == (1 | cast(ulong) address)); 421 } 422 423 424 /// Increment is 0. The LCG portion of the PCG is an MCG. 425 mixin template no_stream(Uint) 426 { 427 /// 428 enum is_mcg = true; 429 /// 430 enum Uint increment = 0; 431 432 /// 433 enum can_specify_stream = false; 434 /// 435 enum size_t streams_pow2 = 0; 436 } 437 438 /// Increment of the LCG portion of the PCG is default_increment. 439 mixin template oneseq_stream(Uint) 440 { 441 /// 442 enum is_mcg = false; 443 /// 444 enum Uint increment = default_increment!Uint; 445 /// 446 enum can_specify_stream = false; 447 /// 448 enum size_t streams_pow2 = 0; 449 } 450 451 /// The increment is dynamically settable and defaults to default_increment!T. 452 mixin template specific_stream(Uint) 453 { 454 /// 455 enum is_mcg = false; 456 /// 457 Uint inc_ = default_increment!Uint; 458 /// 459 alias increment = inc_; 460 /// 461 enum can_specify_stream = true; 462 /// 463 void set_stream()(Uint u) 464 { 465 inc_ = cast(Uint)((u << 1) | 1); 466 } 467 /// 468 enum size_t streams_pow2 = size_t.sizeof*8 -1u; 469 } 470 471 /++ 472 + XorShifts are invertible, but they are someting of a pain to invert. 473 + This function backs them out. It's used by the whacky "inside out" 474 + generator defined later. 475 +/ 476 Uint unxorshift(Uint)(Uint x, size_t bits, size_t shift) 477 { 478 if (2*shift >= bits) { 479 return x ^ (x >> shift); 480 } 481 Uint lowmask1 = (itype(1U) << (bits - shift*2)) - 1; 482 Uint highmask1 = ~lowmask1; 483 Uint top1 = x; 484 Uint bottom1 = x & lowmask1; 485 top1 ^= top1 >> shift; 486 top1 &= highmask1; 487 x = top1 | bottom1; 488 Uint lowmask2 = (itype(1U) << (bits - shift)) - 1; 489 Uint bottom2 = x & lowmask2; 490 bottom2 = unxorshift(bottom2, bits - shift, shift); 491 bottom2 &= lowmask1; 492 return top1 | bottom2; 493 } 494 495 496 /++ 497 + OUTPUT FUNCTIONS. 498 + 499 + These are the core of the PCG generation scheme. They specify how to 500 + turn the base LCG's internal state into the output value of the final 501 + generator. 502 + 503 + All of the classes have code that is written to allow it to be applied 504 + at *arbitrary* bit sizes. 505 +/ 506 507 /++ 508 + XSH RS -- high xorshift, followed by a random shift 509 + 510 + Fast. A good performer. 511 +/ 512 O xsh_rs(O, Uint)(Uint state) 513 { 514 enum bits = Uint.sizeof * 8; 515 enum xtypebits = O.sizeof * 8; 516 enum sparebits = bits - xtypebits; 517 enum opbits = sparebits-5 >= 64 ? 5 518 : sparebits-4 >= 32 ? 4 519 : sparebits-3 >= 16 ? 3 520 : sparebits-2 >= 4 ? 2 521 : sparebits-1 >= 1 ? 1 522 : 0; 523 enum mask = (1 << opbits) - 1; 524 enum maxrandshift = mask; 525 enum topspare = opbits; 526 enum bottomspare = sparebits - topspare; 527 enum xshift = topspare + (xtypebits+maxrandshift)/2; 528 529 auto rshift = opbits ? size_t(state >> (bits - opbits)) & mask : 0; 530 state ^= state >> xshift; 531 O result = O(s >> (bottomspare - maxrandshift + rshift)); 532 return result; 533 } 534 /++ 535 + XSH RR -- high xorshift, followed by a random rotate 536 + 537 + Fast. A good performer. Slightly better statistically than XSH RS. 538 +/ 539 O xsh_rr(O, Uint)(Uint state) 540 { 541 enum bits = Uint.sizeof * 8; 542 enum xtypebits = O.sizeof * 8; 543 enum sparebits = bits - xtypebits; 544 enum wantedopbits = xtypebits >= 128 ? 7 545 : xtypebits >= 64 ? 6 546 : xtypebits >= 32 ? 5 547 : xtypebits >= 16 ? 4 548 : 3; 549 enum opbits = sparebits >= wantedopbits ? wantedopbits : sparebits; 550 enum amplifier = wantedopbits - opbits; 551 enum mask = (1 << opbits) - 1; 552 enum topspare = opbits; 553 enum bottomspare = sparebits - topspare; 554 enum xshift = (topspare + xtypebits)/2; 555 556 auto rot = opbits ? size_t(state >> (bits - opbits)) & mask : 0; 557 auto amprot = (rot << amplifier) & mask; 558 state ^= state >> xshift; 559 O result = cast(O)(state >> bottomspare); 560 import core.bitop: ror; 561 result = ror(result, cast(uint)amprot); 562 return result; 563 } 564 /++ 565 + RXS -- random xorshift 566 +/ 567 O rxs(O, Uint)(Uint state) 568 { 569 enum bits = Uint.sizeof * 8; 570 enum xtypebits = O.sizeof * 8; 571 enum shift = bits - xtypebits; 572 enum extrashift = (xtypebits - shift)/2; 573 enum rshift = shift > 64+8 ? (s >> (bits - 6)) & 63 574 : shift > 32+4 ? (s >> (bits - 5)) & 31 575 : shift > 16+2 ? (s >> (bits - 4)) & 15 576 : shift > 8+1 ? (s >> (bits - 3)) & 7 577 : shift > 4+1 ? (s >> (bits - 2)) & 3 578 : shift > 2+1 ? (s >> (bits - 1)) & 1 579 : 0; 580 state ^= state >> (shift + extrashift - rshift); 581 O result = state >> rshift; 582 return result; 583 } 584 /++ 585 + RXS M XS -- random xorshift, mcg multiply, fixed xorshift 586 + 587 + The most statistically powerful generator, but all those steps 588 + make it slower than some of the others. We give it the rottenest jobs. 589 + 590 + Because it's usually used in contexts where the state type and the 591 + result type are the same, it is a permutation and is thus invertible. 592 + We thus provide a function to invert it. This function is used to 593 + for the "inside out" generator used by the extended generator. 594 +/ 595 O rxs_m_xs_forward(O, Uint)(Uint state) 596 if(is(O == Uint)) 597 { 598 enum bits = Uint.sizeof * 8; 599 enum xtypebits = O.sizeof * 8; 600 enum opbits = xtypebits >= 128 ? 6 601 : xtypebits >= 64 ? 5 602 : xtypebits >= 32 ? 4 603 : xtypebits >= 16 ? 3 604 : 2; 605 enum shift = bits - xtypebits; 606 enum mask = (1 << opbits) - 1; 607 size_t rshift = opbits ? (state >> (bits - opbits)) & mask : 0; 608 state ^= state >> (opbits + rshift); 609 state *= mcg_multiplier!Uint; 610 O result = state >> shift; 611 result ^= result >> ((2U*xtypebits+2U)/3U); 612 return result; 613 } 614 /// ditto 615 O rxs_m_xs_reverse(O, Uint)(Uint state) 616 if(is(O == Uint)) 617 { 618 enum bits = Uint.sizeof * 8; 619 enum opbits = bits >= 128 ? 6 620 : bits >= 64 ? 5 621 : bits >= 32 ? 4 622 : bits >= 16 ? 3 623 : 2; 624 enum mask = (1 << opbits) - 1; 625 626 state = unxorshift(state, bits, (2U*bits+2U)/3U); 627 628 state *= mcg_unmultiplier!Uint; 629 630 auto rshift = opbits ? (state >> (bits - opbits)) & mask : 0; 631 state = unxorshift(s, bits, opbits + rshift); 632 633 return s; 634 } 635 /++ 636 + XSL RR -- fixed xorshift (to low bits), random rotate 637 + 638 + Useful for 128-bit types that are split across two CPU registers. 639 +/ 640 O xsl_rr(O, Uint)(Uint state) 641 { 642 enum bits = Uint.sizeof * 8; 643 enum xtypebits = O.sizeof * 8; 644 enum sparebits = bits - xtypebits; 645 enum wantedopbits = xtypebits >= 128 ? 7 646 : xtypebits >= 64 ? 6 647 : xtypebits >= 32 ? 5 648 : xtypebits >= 16 ? 4 649 : 3; 650 enum opbits = sparebits >= wantedopbits ? wantedopbits : sparebits; 651 enum amplifier = wantedopbits - opbits; 652 enum mask = (1 << opbits) - 1; 653 enum topspare = sparebits; 654 enum bottomspare = sparebits - topspare; 655 enum xshift = (topspare + xtypebits) / 2; 656 657 auto rot = opbits ? size_t(state >> (bits - opbits)) & mask : 0; 658 auto amprot = (rot << amplifier) & mask; 659 state ^= state >> xshift; 660 O result = state >> bottomspare; 661 result = rotr(result, amprot); 662 return result; 663 } 664 665 private template half_size(Uint) 666 { 667 static if (is(Uint == ucent)) 668 alias half_size = ulong; 669 else static if (is(Uint == ulong)) 670 alias half_size = uint; 671 else static if (is(Uint == uint)) 672 alias half_size = ushort; 673 else static if (is(Uint == ushort)) 674 alias half_size = ubyte; 675 else 676 static assert(0); 677 } 678 679 /++ 680 + XSL RR RR -- fixed xorshift (to low bits), random rotate (both parts) 681 + 682 + Useful for 128-bit types that are split across two CPU registers. 683 + If you really want an invertible 128-bit RNG, I guess this is the one. 684 +/ 685 686 O xsl_rr_rr(O, Uint)(Uint state) 687 if(is(O == Uint)) 688 { 689 alias H = half_size!Uint; 690 enum htypebits = H.sizeof * 8; 691 enum bits = Uint.sizeof * 8; 692 enum sparebits = bits - htypebits; 693 enum wantedopbits = htypebits >= 128 ? 7 694 : htypebits >= 64 ? 6 695 : htypebits >= 32 ? 5 696 : htypebits >= 16 ? 4 697 : 3; 698 enum opbits = sparebits >= wantedopbits ? wantedopbits : sparebits; 699 enum amplifier = wantedopbits - opbits; 700 enum mask = (1 << opbits) - 1; 701 enum topspare = sparebits; 702 enum xshift = (topspare + htypebits) / 2; 703 704 auto rot = opbits ? size_t(s >> (bits - opbits)) & mask : 0; 705 auto amprot = (rot << amplifier) & mask; 706 state ^= state >> xshift; 707 H lowbits = cast(H)s; 708 lowbits = rotr(lowbits, amprot); 709 H highbits = cast(H)(s >> topspare); 710 auto rot2 = lowbits & mask; 711 auto amprot2 = (rot2 << amplifier) & mask; 712 highbits = rotr(highbits, amprot2); 713 return (O(highbits) << topspare) ^ O(lowbits); 714 715 } 716 717 /++ 718 + XSH -- fixed xorshift (to high bits) 719 + 720 + Not available at 64-bits or less. 721 +/ 722 723 O xsh(O, Uint)(Uint state) if(Uint.sizeof > 8) 724 { 725 enum bits = Uint.sizeof * 8; 726 enum xtypebits = O.sizeof * 8; 727 enum sparebits = bits - xtypebits; 728 enum topspare = 0; 729 enum bottomspare = sparebits - topspare; 730 enum xshift = (topspare + xtypebits) / 2; 731 732 state ^= state >> xshift; 733 O result = state >> bottomspare; 734 return result; 735 } 736 737 /++ 738 + XSL -- fixed xorshift (to low bits) 739 + 740 + Not available at 64-bits or less. 741 +/ 742 743 O xsl(O, Uint)(Uint state) if(Uint.sizeof > 8) 744 { 745 enum bits = Uint.sizeof * 8; 746 enum xtypebits = O.sizeof * 8; 747 enum sparebits = bits - xtypebits; 748 enum topspare = sparebits; 749 enum bottomspare = sparebits - topspare; 750 enum xshift = (topspare + xtypebits) / 2; 751 752 state ^= state >> xshift; 753 O result = state >> bottomspare; 754 return result; 755 } 756 757 private alias AliasSeq(T...) = T; 758 759 @safe version(mir_random_test) unittest 760 { 761 762 foreach(RNG; AliasSeq!(pcg32,pcg32_unique,pcg32_oneseq,pcg32_fast, 763 pcg8_once_insecure,pcg16_once_insecure,pcg32_once_insecure,pcg64_once_insecure, 764 pcg8_oneseq_once_insecure,pcg16_oneseq_once_insecure,pcg32_oneseq_once_insecure, 765 pcg64_oneseq_once_insecure)) 766 { 767 static assert(isSaturatedRandomEngine!RNG); 768 auto gen = RNG(cast(RNG.Uint)unpredictableSeed); 769 gen(); 770 } 771 } 772 @safe version(mir_random_test) unittest 773 { 774 auto gen = pcg32(0x198c8585); 775 gen.skip(1000); 776 assert(gen()== 0xd187a760); 777 auto gen2 = pcg32(0x198c8585); 778 779 foreach(_; 0 .. 1000) 780 gen2(); 781 assert(gen2() == 0xd187a760); 782 assert(gen() == gen2()); 783 } 784 785 @nogc nothrow pure @safe unittest 786 { 787 foreach (ShouldHaveStaticBump; AliasSeq!(pcg32_oneseq, pcg32_fast, pcg32_oneseq_once_insecure)) 788 static assert (__traits(compiles, { enum e = ShouldHaveStaticBump.bump(1u); })); 789 790 foreach (ShouldLackStaticBump; AliasSeq!(pcg32, pcg32_unique, pcg32_once_insecure)) 791 static assert (!__traits(compiles, { enum e = ShouldLackStaticBump.bump(1u); })); 792 }