1 // Written in the D programming language. 2 3 /** 4 Facilities for random number generation. 5 6 $(SCRIPT inhibitQuickIndex = 1;) 7 $(DIVC quickindex, 8 $(BOOKTABLE, 9 $(TR $(TH Category) $(TH Functions)) 10 $(TR $(TD Uniform sampling) $(TD 11 $(LREF uniform) 12 $(LREF uniform01) 13 $(LREF uniformDistribution) 14 )) 15 $(TR $(TD Element sampling) $(TD 16 $(LREF choice) 17 $(LREF dice) 18 )) 19 $(TR $(TD Range sampling) $(TD 20 $(LREF randomCover) 21 $(LREF randomSample) 22 )) 23 $(TR $(TD Default Random Engines) $(TD 24 $(LREF rndGen) 25 $(LREF Random) 26 $(LREF unpredictableSeed) 27 )) 28 $(TR $(TD Linear Congruential Engines) $(TD 29 $(LREF MinstdRand) 30 $(LREF MinstdRand0) 31 $(LREF LinearCongruentialEngine) 32 )) 33 $(TR $(TD Mersenne Twister Engines) $(TD 34 $(LREF Mt19937) 35 $(LREF Mt19937_64) 36 $(LREF MersenneTwisterEngine) 37 )) 38 $(TR $(TD Xorshift Engines) $(TD 39 $(LREF Xorshift) 40 $(LREF XorshiftEngine) 41 $(LREF Xorshift32) 42 $(LREF Xorshift64) 43 $(LREF Xorshift96) 44 $(LREF Xorshift128) 45 $(LREF Xorshift160) 46 $(LREF Xorshift192) 47 )) 48 $(TR $(TD Shuffle) $(TD 49 $(LREF partialShuffle) 50 $(LREF randomShuffle) 51 )) 52 $(TR $(TD Traits) $(TD 53 $(LREF isSeedable) 54 $(LREF isUniformRNG) 55 )) 56 )) 57 58 $(RED Disclaimer:) The random number generators and API provided in this 59 module are not designed to be cryptographically secure, and are therefore 60 unsuitable for cryptographic or security-related purposes such as generating 61 authentication tokens or network sequence numbers. For such needs, please use a 62 reputable cryptographic library instead. 63 64 The new-style generator objects hold their own state so they are 65 immune of threading issues. The generators feature a number of 66 well-known and well-documented methods of generating random 67 numbers. An overall fast and reliable means to generate random numbers 68 is the $(D_PARAM Mt19937) generator, which derives its name from 69 "$(LINK2 https://en.wikipedia.org/wiki/Mersenne_Twister, Mersenne Twister) 70 with a period of 2 to the power of 71 19937". In memory-constrained situations, 72 $(LINK2 https://en.wikipedia.org/wiki/Linear_congruential_generator, 73 linear congruential generators) such as `MinstdRand0` and `MinstdRand` might be 74 useful. The standard library provides an alias $(D_PARAM Random) for 75 whichever generator it considers the most fit for the target 76 environment. 77 78 In addition to random number generators, this module features 79 distributions, which skew a generator's output statistical 80 distribution in various ways. So far the uniform distribution for 81 integers and real numbers have been implemented. 82 83 Source: $(PHOBOSSRC std/random.d) 84 85 Macros: 86 87 Copyright: Copyright Andrei Alexandrescu 2008 - 2009, Joseph Rushton Wakeling 2012. 88 License: $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0). 89 Authors: $(HTTP erdani.org, Andrei Alexandrescu) 90 Masahiro Nakagawa (Xorshift random generator) 91 $(HTTP braingam.es, Joseph Rushton Wakeling) (Algorithm D for random sampling) 92 Ilya Yaroshenko (Mersenne Twister implementation, adapted from $(HTTPS github.com/libmir/mir-random, mir-random)) 93 Credits: The entire random number library architecture is derived from the 94 excellent $(HTTP open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2461.pdf, C++0X) 95 random number facility proposed by Jens Maurer and contributed to by 96 researchers at the Fermi laboratory (excluding Xorshift). 97 */ 98 /* 99 Copyright Andrei Alexandrescu 2008 - 2009. 100 Distributed under the Boost Software License, Version 1.0. 101 (See accompanying file LICENSE_1_0.txt or copy at 102 http://www.boost.org/LICENSE_1_0.txt) 103 */ 104 module std.random; 105 106 107 import std.range.primitives; 108 import std.traits; 109 110 /// 111 @safe unittest 112 { 113 import std.algorithm.comparison : among, equal; 114 import std.range : iota; 115 116 // seed a random generator with a constant 117 auto rnd = Random(42); 118 119 // Generate a uniformly-distributed integer in the range [0, 14] 120 // If no random generator is passed, the global `rndGen` would be used 121 auto i = uniform(0, 15, rnd); 122 assert(i >= 0 && i < 15); 123 124 // Generate a uniformly-distributed real in the range [0, 100) 125 auto r = uniform(0.0L, 100.0L, rnd); 126 assert(r >= 0 && r < 100); 127 128 // Sample from a custom type 129 enum Fruit { apple, mango, pear } 130 auto f = rnd.uniform!Fruit; 131 with(Fruit) 132 assert(f.among(apple, mango, pear)); 133 134 // Generate a 32-bit random number 135 auto u = uniform!uint(rnd); 136 static assert(is(typeof(u) == uint)); 137 138 // Generate a random number in the range in the range [0, 1) 139 auto u2 = uniform01(rnd); 140 assert(u2 >= 0 && u2 < 1); 141 142 // Select an element randomly 143 auto el = 10.iota.choice(rnd); 144 assert(0 <= el && el < 10); 145 146 // Throw a dice with custom proportions 147 // 0: 20%, 1: 10%, 2: 60% 148 auto val = rnd.dice(0.2, 0.1, 0.6); 149 assert(0 <= val && val <= 2); 150 151 auto rnd2 = MinstdRand0(42); 152 153 // Select a random subsample from a range 154 assert(10.iota.randomSample(3, rnd2).equal([7, 8, 9])); 155 156 // Cover all elements in an array in random order 157 version (D_LP64) // https://issues.dlang.org/show_bug.cgi?id=15147 158 assert(10.iota.randomCover(rnd2).equal([7, 4, 2, 0, 1, 6, 8, 3, 9, 5])); 159 else 160 assert(10.iota.randomCover(rnd2).equal([4, 8, 7, 3, 5, 9, 2, 6, 0, 1])); 161 162 // Shuffle an array 163 version (D_LP64) // https://issues.dlang.org/show_bug.cgi?id=15147 164 assert([0, 1, 2, 4, 5].randomShuffle(rnd2).equal([2, 0, 4, 5, 1])); 165 else 166 assert([0, 1, 2, 4, 5].randomShuffle(rnd2).equal([4, 2, 5, 0, 1])); 167 } 168 169 version (OSX) 170 version = Darwin; 171 else version (iOS) 172 version = Darwin; 173 else version (TVOS) 174 version = Darwin; 175 else version (WatchOS) 176 version = Darwin; 177 178 version (D_InlineAsm_X86) version = InlineAsm_X86_Any; 179 version (D_InlineAsm_X86_64) version = InlineAsm_X86_Any; 180 181 version (StdUnittest) 182 { 183 static import std.meta; 184 package alias Xorshift64_64 = XorshiftEngine!(ulong, 64, -12, 25, -27); 185 package alias Xorshift128_64 = XorshiftEngine!(ulong, 128, 23, -18, -5); 186 package alias PseudoRngTypes = std.meta.AliasSeq!(MinstdRand0, MinstdRand, Mt19937, Xorshift32, Xorshift64, 187 Xorshift96, Xorshift128, Xorshift160, Xorshift192, 188 Xorshift64_64, Xorshift128_64); 189 } 190 191 // Segments of the code in this file Copyright (c) 1997 by Rick Booth 192 // From "Inner Loops" by Rick Booth, Addison-Wesley 193 194 // Work derived from: 195 196 /* 197 A C-program for MT19937, with initialization improved 2002/1/26. 198 Coded by Takuji Nishimura and Makoto Matsumoto. 199 200 Before using, initialize the state by using init_genrand(seed) 201 or init_by_array(init_key, key_length). 202 203 Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, 204 All rights reserved. 205 206 Redistribution and use in source and binary forms, with or without 207 modification, are permitted provided that the following conditions 208 are met: 209 210 1. Redistributions of source code must retain the above copyright 211 notice, this list of conditions and the following disclaimer. 212 213 2. Redistributions in binary form must reproduce the above copyright 214 notice, this list of conditions and the following disclaimer in the 215 documentation and/or other materials provided with the distribution. 216 217 3. The names of its contributors may not be used to endorse or promote 218 products derived from this software without specific prior written 219 permission. 220 221 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 222 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 223 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 224 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 225 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 226 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 227 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 228 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 229 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 230 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 231 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 232 233 234 Any feedback is very welcome. 235 http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html 236 email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) 237 */ 238 239 /** 240 * Test if Rng is a random-number generator. The overload 241 * taking a ElementType also makes sure that the Rng generates 242 * values of that type. 243 * 244 * A random-number generator has at least the following features: 245 * $(UL 246 * $(LI it's an InputRange) 247 * $(LI it has a 'bool isUniformRandom' field readable in CTFE) 248 * ) 249 */ 250 template isUniformRNG(Rng, ElementType) 251 { 252 enum bool isUniformRNG = .isUniformRNG!Rng && 253 is(std.range.primitives.ElementType!Rng == ElementType); 254 } 255 256 /** 257 * ditto 258 */ 259 template isUniformRNG(Rng) 260 { 261 enum bool isUniformRNG = isInputRange!Rng && 262 is(typeof( 263 { 264 static assert(Rng.isUniformRandom); //tag 265 })); 266 } 267 268 /// 269 @safe unittest 270 { 271 struct NoRng 272 { 273 @property uint front() {return 0;} 274 @property bool empty() {return false;} 275 void popFront() {} 276 } 277 static assert(!isUniformRNG!(NoRng)); 278 279 struct validRng 280 { 281 @property uint front() {return 0;} 282 @property bool empty() {return false;} 283 void popFront() {} 284 285 enum isUniformRandom = true; 286 } 287 static assert(isUniformRNG!(validRng, uint)); 288 static assert(isUniformRNG!(validRng)); 289 } 290 291 @safe unittest 292 { 293 // two-argument predicate should not require @property on `front` 294 // https://issues.dlang.org/show_bug.cgi?id=19837 295 struct Rng 296 { 297 float front() {return 0;} 298 void popFront() {} 299 enum empty = false; 300 enum isUniformRandom = true; 301 } 302 static assert(isUniformRNG!(Rng, float)); 303 } 304 305 /** 306 * Test if Rng is seedable. The overload 307 * taking a SeedType also makes sure that the Rng can be seeded with SeedType. 308 * 309 * A seedable random-number generator has the following additional features: 310 * $(UL 311 * $(LI it has a 'seed(ElementType)' function) 312 * ) 313 */ 314 template isSeedable(Rng, SeedType) 315 { 316 enum bool isSeedable = isUniformRNG!(Rng) && 317 is(typeof( 318 { 319 Rng r = void; // can define a Rng object 320 SeedType s = void; 321 r.seed(s); // can seed a Rng 322 })); 323 } 324 325 ///ditto 326 template isSeedable(Rng) 327 { 328 enum bool isSeedable = isUniformRNG!Rng && 329 is(typeof( 330 { 331 Rng r = void; // can define a Rng object 332 alias SeedType = typeof(r.front); 333 SeedType s = void; 334 r.seed(s); // can seed a Rng 335 })); 336 } 337 338 /// 339 @safe unittest 340 { 341 struct validRng 342 { 343 @property uint front() {return 0;} 344 @property bool empty() {return false;} 345 void popFront() {} 346 347 enum isUniformRandom = true; 348 } 349 static assert(!isSeedable!(validRng, uint)); 350 static assert(!isSeedable!(validRng)); 351 352 struct seedRng 353 { 354 @property uint front() {return 0;} 355 @property bool empty() {return false;} 356 void popFront() {} 357 void seed(uint val){} 358 enum isUniformRandom = true; 359 } 360 static assert(isSeedable!(seedRng, uint)); 361 static assert(!isSeedable!(seedRng, ulong)); 362 static assert(isSeedable!(seedRng)); 363 } 364 365 @safe @nogc pure nothrow unittest 366 { 367 struct NoRng 368 { 369 @property uint front() {return 0;} 370 @property bool empty() {return false;} 371 void popFront() {} 372 } 373 static assert(!isUniformRNG!(NoRng, uint)); 374 static assert(!isUniformRNG!(NoRng)); 375 static assert(!isSeedable!(NoRng, uint)); 376 static assert(!isSeedable!(NoRng)); 377 378 struct NoRng2 379 { 380 @property uint front() {return 0;} 381 @property bool empty() {return false;} 382 void popFront() {} 383 384 enum isUniformRandom = false; 385 } 386 static assert(!isUniformRNG!(NoRng2, uint)); 387 static assert(!isUniformRNG!(NoRng2)); 388 static assert(!isSeedable!(NoRng2, uint)); 389 static assert(!isSeedable!(NoRng2)); 390 391 struct NoRng3 392 { 393 @property bool empty() {return false;} 394 void popFront() {} 395 396 enum isUniformRandom = true; 397 } 398 static assert(!isUniformRNG!(NoRng3, uint)); 399 static assert(!isUniformRNG!(NoRng3)); 400 static assert(!isSeedable!(NoRng3, uint)); 401 static assert(!isSeedable!(NoRng3)); 402 403 struct validRng 404 { 405 @property uint front() {return 0;} 406 @property bool empty() {return false;} 407 void popFront() {} 408 409 enum isUniformRandom = true; 410 } 411 static assert(isUniformRNG!(validRng, uint)); 412 static assert(isUniformRNG!(validRng)); 413 static assert(!isSeedable!(validRng, uint)); 414 static assert(!isSeedable!(validRng)); 415 416 struct seedRng 417 { 418 @property uint front() {return 0;} 419 @property bool empty() {return false;} 420 void popFront() {} 421 void seed(uint val){} 422 enum isUniformRandom = true; 423 } 424 static assert(isUniformRNG!(seedRng, uint)); 425 static assert(isUniformRNG!(seedRng)); 426 static assert(isSeedable!(seedRng, uint)); 427 static assert(!isSeedable!(seedRng, ulong)); 428 static assert(isSeedable!(seedRng)); 429 } 430 431 /** 432 Linear Congruential generator. When m = 0, no modulus is used. 433 */ 434 struct LinearCongruentialEngine(UIntType, UIntType a, UIntType c, UIntType m) 435 if (isUnsigned!UIntType) 436 { 437 ///Mark this as a Rng 438 enum bool isUniformRandom = true; 439 /// Does this generator have a fixed range? ($(D_PARAM true)). 440 enum bool hasFixedRange = true; 441 /// Lowest generated value (`1` if $(D c == 0), `0` otherwise). 442 enum UIntType min = ( c == 0 ? 1 : 0 ); 443 /// Highest generated value ($(D modulus - 1)). 444 enum UIntType max = m - 1; 445 /** 446 The parameters of this distribution. The random number is $(D_PARAM x 447 = (x * multipler + increment) % modulus). 448 */ 449 enum UIntType multiplier = a; 450 ///ditto 451 enum UIntType increment = c; 452 ///ditto 453 enum UIntType modulus = m; 454 455 static assert(isIntegral!(UIntType)); 456 static assert(m == 0 || a < m); 457 static assert(m == 0 || c < m); 458 static assert(m == 0 || 459 (cast(ulong) a * (m-1) + c) % m == (c < a ? c - a + m : c - a)); 460 461 // Check for maximum range 462 private static ulong gcd(ulong a, ulong b) @safe pure nothrow @nogc 463 { 464 while (b) 465 { 466 auto t = b; 467 b = a % b; 468 a = t; 469 } 470 return a; 471 } 472 473 private static ulong primeFactorsOnly(ulong n) @safe pure nothrow @nogc 474 { 475 ulong result = 1; 476 ulong iter = 2; 477 for (; n >= iter * iter; iter += 2 - (iter == 2)) 478 { 479 if (n % iter) continue; 480 result *= iter; 481 do 482 { 483 n /= iter; 484 } while (n % iter == 0); 485 } 486 return result * n; 487 } 488 489 @safe pure nothrow unittest 490 { 491 static assert(primeFactorsOnly(100) == 10); 492 //writeln(primeFactorsOnly(11)); 493 static assert(primeFactorsOnly(11) == 11); 494 static assert(primeFactorsOnly(7 * 7 * 7 * 11 * 15 * 11) == 7 * 11 * 15); 495 static assert(primeFactorsOnly(129 * 2) == 129 * 2); 496 // enum x = primeFactorsOnly(7 * 7 * 7 * 11 * 15); 497 // static assert(x == 7 * 11 * 15); 498 } 499 500 private static bool properLinearCongruentialParameters(ulong m, 501 ulong a, ulong c) @safe pure nothrow @nogc 502 { 503 if (m == 0) 504 { 505 static if (is(UIntType == uint)) 506 { 507 // Assume m is uint.max + 1 508 m = (1uL << 32); 509 } 510 else 511 { 512 return false; 513 } 514 } 515 // Bounds checking 516 if (a == 0 || a >= m || c >= m) return false; 517 // c and m are relatively prime 518 if (c > 0 && gcd(c, m) != 1) return false; 519 // a - 1 is divisible by all prime factors of m 520 if ((a - 1) % primeFactorsOnly(m)) return false; 521 // if a - 1 is multiple of 4, then m is a multiple of 4 too. 522 if ((a - 1) % 4 == 0 && m % 4) return false; 523 // Passed all tests 524 return true; 525 } 526 527 // check here 528 static assert(c == 0 || properLinearCongruentialParameters(m, a, c), 529 "Incorrect instantiation of LinearCongruentialEngine"); 530 531 /** 532 Constructs a $(D_PARAM LinearCongruentialEngine) generator seeded with 533 `x0`. 534 */ 535 this(UIntType x0) @safe pure nothrow @nogc 536 { 537 seed(x0); 538 } 539 540 /** 541 (Re)seeds the generator. 542 */ 543 void seed(UIntType x0 = 1) @safe pure nothrow @nogc 544 { 545 _x = modulus ? (x0 % modulus) : x0; 546 static if (c == 0) 547 { 548 //Necessary to prevent generator from outputting an endless series of zeroes. 549 if (_x == 0) 550 _x = max; 551 } 552 popFront(); 553 } 554 555 /** 556 Advances the random sequence. 557 */ 558 void popFront() @safe pure nothrow @nogc 559 { 560 static if (m) 561 { 562 static if (is(UIntType == uint) && m == uint.max) 563 { 564 immutable ulong 565 x = (cast(ulong) a * _x + c), 566 v = x >> 32, 567 w = x & uint.max; 568 immutable y = cast(uint)(v + w); 569 _x = (y < v || y == uint.max) ? (y + 1) : y; 570 } 571 else static if (is(UIntType == uint) && m == int.max) 572 { 573 immutable ulong 574 x = (cast(ulong) a * _x + c), 575 v = x >> 31, 576 w = x & int.max; 577 immutable uint y = cast(uint)(v + w); 578 _x = (y >= int.max) ? (y - int.max) : y; 579 } 580 else 581 { 582 _x = cast(UIntType) ((cast(ulong) a * _x + c) % m); 583 } 584 } 585 else 586 { 587 _x = a * _x + c; 588 } 589 } 590 591 /** 592 Returns the current number in the random sequence. 593 */ 594 @property UIntType front() const @safe pure nothrow @nogc 595 { 596 return _x; 597 } 598 599 /// 600 @property typeof(this) save() const @safe pure nothrow @nogc 601 { 602 return this; 603 } 604 605 /** 606 Always `false` (random generators are infinite ranges). 607 */ 608 enum bool empty = false; 609 610 // https://issues.dlang.org/show_bug.cgi?id=21610 611 static if (m) 612 { 613 private UIntType _x = (a + c) % m; 614 } 615 else 616 { 617 private UIntType _x = a + c; 618 } 619 } 620 621 /// Declare your own linear congruential engine 622 @safe unittest 623 { 624 alias CPP11LCG = LinearCongruentialEngine!(uint, 48271, 0, 2_147_483_647); 625 626 // seed with a constant 627 auto rnd = CPP11LCG(42); 628 auto n = rnd.front; // same for each run 629 assert(n == 2027382); 630 } 631 632 /// Declare your own linear congruential engine 633 @safe unittest 634 { 635 // glibc's LCG 636 alias GLibcLCG = LinearCongruentialEngine!(uint, 1103515245, 12345, 2_147_483_648); 637 638 // Seed with an unpredictable value 639 auto rnd = GLibcLCG(unpredictableSeed); 640 auto n = rnd.front; // different across runs 641 } 642 643 /// Declare your own linear congruential engine 644 @safe unittest 645 { 646 // Visual C++'s LCG 647 alias MSVCLCG = LinearCongruentialEngine!(uint, 214013, 2531011, 0); 648 649 // seed with a constant 650 auto rnd = MSVCLCG(1); 651 auto n = rnd.front; // same for each run 652 assert(n == 2745024); 653 } 654 655 // Ensure that unseeded LCGs produce correct values 656 @safe unittest 657 { 658 alias LGE = LinearCongruentialEngine!(uint, 10000, 19682, 19683); 659 660 LGE rnd; 661 assert(rnd.front == 9999); 662 } 663 664 /** 665 Define $(D_PARAM LinearCongruentialEngine) generators with well-chosen 666 parameters. `MinstdRand0` implements Park and Miller's "minimal 667 standard" $(HTTP 668 wikipedia.org/wiki/Park%E2%80%93Miller_random_number_generator, 669 generator) that uses 16807 for the multiplier. `MinstdRand` 670 implements a variant that has slightly better spectral behavior by 671 using the multiplier 48271. Both generators are rather simplistic. 672 */ 673 alias MinstdRand0 = LinearCongruentialEngine!(uint, 16_807, 0, 2_147_483_647); 674 /// ditto 675 alias MinstdRand = LinearCongruentialEngine!(uint, 48_271, 0, 2_147_483_647); 676 677 /// 678 @safe @nogc unittest 679 { 680 // seed with a constant 681 auto rnd0 = MinstdRand0(1); 682 auto n = rnd0.front; 683 // same for each run 684 assert(n == 16807); 685 686 // Seed with an unpredictable value 687 rnd0.seed(unpredictableSeed); 688 n = rnd0.front; // different across runs 689 } 690 691 @safe @nogc unittest 692 { 693 import std.range; 694 static assert(isForwardRange!MinstdRand); 695 static assert(isUniformRNG!MinstdRand); 696 static assert(isUniformRNG!MinstdRand0); 697 static assert(isUniformRNG!(MinstdRand, uint)); 698 static assert(isUniformRNG!(MinstdRand0, uint)); 699 static assert(isSeedable!MinstdRand); 700 static assert(isSeedable!MinstdRand0); 701 static assert(isSeedable!(MinstdRand, uint)); 702 static assert(isSeedable!(MinstdRand0, uint)); 703 704 // The correct numbers are taken from The Database of Integer Sequences 705 // http://www.research.att.com/~njas/sequences/eisBTfry00128.txt 706 enum ulong[20] checking0 = [ 707 16807UL,282475249,1622650073,984943658,1144108930,470211272, 708 101027544,1457850878,1458777923,2007237709,823564440,1115438165, 709 1784484492,74243042,114807987,1137522503,1441282327,16531729, 710 823378840,143542612 ]; 711 //auto rnd0 = MinstdRand0(1); 712 MinstdRand0 rnd0; 713 714 foreach (e; checking0) 715 { 716 assert(rnd0.front == e); 717 rnd0.popFront(); 718 } 719 // Test the 10000th invocation 720 // Correct value taken from: 721 // http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2461.pdf 722 rnd0.seed(); 723 popFrontN(rnd0, 9999); 724 assert(rnd0.front == 1043618065); 725 726 // Test MinstdRand 727 enum ulong[6] checking = [48271UL,182605794,1291394886,1914720637,2078669041, 728 407355683]; 729 //auto rnd = MinstdRand(1); 730 MinstdRand rnd; 731 foreach (e; checking) 732 { 733 assert(rnd.front == e); 734 rnd.popFront(); 735 } 736 737 // Test the 10000th invocation 738 // Correct value taken from: 739 // http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2461.pdf 740 rnd.seed(); 741 popFrontN(rnd, 9999); 742 assert(rnd.front == 399268537); 743 744 // Check .save works 745 static foreach (Type; std.meta.AliasSeq!(MinstdRand0, MinstdRand)) 746 {{ 747 auto rnd1 = Type(123_456_789); 748 rnd1.popFront(); 749 // https://issues.dlang.org/show_bug.cgi?id=15853 750 auto rnd2 = ((const ref Type a) => a.save())(rnd1); 751 assert(rnd1 == rnd2); 752 // Enable next test when RNGs are reference types 753 version (none) { assert(rnd1 !is rnd2); } 754 for (auto i = 0; i < 3; i++, rnd1.popFront, rnd2.popFront) 755 assert(rnd1.front() == rnd2.front()); 756 }} 757 } 758 759 @safe @nogc unittest 760 { 761 auto rnd0 = MinstdRand0(MinstdRand0.modulus); 762 auto n = rnd0.front; 763 rnd0.popFront(); 764 assert(n != rnd0.front); 765 } 766 767 /** 768 The $(LINK2 https://en.wikipedia.org/wiki/Mersenne_Twister, Mersenne Twister) generator. 769 */ 770 struct MersenneTwisterEngine(UIntType, size_t w, size_t n, size_t m, size_t r, 771 UIntType a, size_t u, UIntType d, size_t s, 772 UIntType b, size_t t, 773 UIntType c, size_t l, UIntType f) 774 if (isUnsigned!UIntType) 775 { 776 static assert(0 < w && w <= UIntType.sizeof * 8); 777 static assert(1 <= m && m <= n); 778 static assert(0 <= r && 0 <= u && 0 <= s && 0 <= t && 0 <= l); 779 static assert(r <= w && u <= w && s <= w && t <= w && l <= w); 780 static assert(0 <= a && 0 <= b && 0 <= c); 781 static assert(n <= ptrdiff_t.max); 782 783 ///Mark this as a Rng 784 enum bool isUniformRandom = true; 785 786 /** 787 Parameters for the generator. 788 */ 789 enum size_t wordSize = w; 790 enum size_t stateSize = n; /// ditto 791 enum size_t shiftSize = m; /// ditto 792 enum size_t maskBits = r; /// ditto 793 enum UIntType xorMask = a; /// ditto 794 enum size_t temperingU = u; /// ditto 795 enum UIntType temperingD = d; /// ditto 796 enum size_t temperingS = s; /// ditto 797 enum UIntType temperingB = b; /// ditto 798 enum size_t temperingT = t; /// ditto 799 enum UIntType temperingC = c; /// ditto 800 enum size_t temperingL = l; /// ditto 801 enum UIntType initializationMultiplier = f; /// ditto 802 803 /// Smallest generated value (0). 804 enum UIntType min = 0; 805 /// Largest generated value. 806 enum UIntType max = UIntType.max >> (UIntType.sizeof * 8u - w); 807 // note, `max` also serves as a bitmask for the lowest `w` bits 808 static assert(a <= max && b <= max && c <= max && f <= max); 809 810 /// The default seed value. 811 enum UIntType defaultSeed = 5489u; 812 813 // Bitmasks used in the 'twist' part of the algorithm 814 private enum UIntType lowerMask = (cast(UIntType) 1u << r) - 1, 815 upperMask = (~lowerMask) & max; 816 817 /* 818 Collection of all state variables 819 used by the generator 820 */ 821 private struct State 822 { 823 /* 824 State array of the generator. This 825 is iterated through backwards (from 826 last element to first), providing a 827 few extra compiler optimizations by 828 comparison to the forward iteration 829 used in most implementations. 830 */ 831 UIntType[n] data; 832 833 /* 834 Cached copy of most recently updated 835 element of `data` state array, ready 836 to be tempered to generate next 837 `front` value 838 */ 839 UIntType z; 840 841 /* 842 Most recently generated random variate 843 */ 844 UIntType front; 845 846 /* 847 Index of the entry in the `data` 848 state array that will be twisted 849 in the next `popFront()` call 850 */ 851 size_t index; 852 } 853 854 /* 855 State variables used by the generator; 856 initialized to values equivalent to 857 explicitly seeding the generator with 858 `defaultSeed` 859 */ 860 private State state = defaultState(); 861 /* NOTE: the above is a workaround to ensure 862 backwards compatibility with the original 863 implementation, which permitted implicit 864 construction. With `@disable this();` 865 it would not be necessary. */ 866 867 /** 868 Constructs a MersenneTwisterEngine object. 869 */ 870 this(UIntType value) @safe pure nothrow @nogc 871 { 872 seed(value); 873 } 874 875 /** 876 Generates the default initial state for a Mersenne 877 Twister; equivalent to the internal state obtained 878 by calling `seed(defaultSeed)` 879 */ 880 private static State defaultState() @safe pure nothrow @nogc 881 { 882 if (!__ctfe) assert(false); 883 State mtState; 884 seedImpl(defaultSeed, mtState); 885 return mtState; 886 } 887 888 /** 889 Seeds a MersenneTwisterEngine object. 890 Note: 891 This seed function gives 2^w starting points (the lowest w bits of 892 the value provided will be used). To allow the RNG to be started 893 in any one of its internal states use the seed overload taking an 894 InputRange. 895 */ 896 void seed()(UIntType value = defaultSeed) @safe pure nothrow @nogc 897 { 898 this.seedImpl(value, this.state); 899 } 900 901 /** 902 Implementation of the seeding mechanism, which 903 can be used with an arbitrary `State` instance 904 */ 905 private static void seedImpl(UIntType value, ref State mtState) @nogc 906 { 907 mtState.data[$ - 1] = value; 908 static if (max != UIntType.max) 909 { 910 mtState.data[$ - 1] &= max; 911 } 912 913 foreach_reverse (size_t i, ref e; mtState.data[0 .. $ - 1]) 914 { 915 e = f * (mtState.data[i + 1] ^ (mtState.data[i + 1] >> (w - 2))) + cast(UIntType)(n - (i + 1)); 916 static if (max != UIntType.max) 917 { 918 e &= max; 919 } 920 } 921 922 mtState.index = n - 1; 923 924 /* double popFront() to guarantee both `mtState.z` 925 and `mtState.front` are derived from the newly 926 set values in `mtState.data` */ 927 MersenneTwisterEngine.popFrontImpl(mtState); 928 MersenneTwisterEngine.popFrontImpl(mtState); 929 } 930 931 /** 932 Seeds a MersenneTwisterEngine object using an InputRange. 933 934 Throws: 935 `Exception` if the InputRange didn't provide enough elements to seed the generator. 936 The number of elements required is the 'n' template parameter of the MersenneTwisterEngine struct. 937 */ 938 void seed(T)(T range) if (isInputRange!T && is(immutable ElementType!T == immutable UIntType)) 939 { 940 this.seedImpl(range, this.state); 941 } 942 943 /** 944 Implementation of the range-based seeding mechanism, 945 which can be used with an arbitrary `State` instance 946 */ 947 private static void seedImpl(T)(T range, ref State mtState) 948 if (isInputRange!T && is(immutable ElementType!T == immutable UIntType)) 949 { 950 size_t j; 951 for (j = 0; j < n && !range.empty; ++j, range.popFront()) 952 { 953 ptrdiff_t idx = n - j - 1; 954 mtState.data[idx] = range.front; 955 } 956 957 mtState.index = n - 1; 958 959 if (range.empty && j < n) 960 { 961 import core.internal.string : UnsignedStringBuf, unsignedToTempString; 962 963 UnsignedStringBuf buf = void; 964 string s = "MersenneTwisterEngine.seed: Input range didn't provide enough elements: Need "; 965 s ~= unsignedToTempString(n, buf) ~ " elements."; 966 throw new Exception(s); 967 } 968 969 /* double popFront() to guarantee both `mtState.z` 970 and `mtState.front` are derived from the newly 971 set values in `mtState.data` */ 972 MersenneTwisterEngine.popFrontImpl(mtState); 973 MersenneTwisterEngine.popFrontImpl(mtState); 974 } 975 976 /** 977 Advances the generator. 978 */ 979 void popFront() @safe pure nothrow @nogc 980 { 981 this.popFrontImpl(this.state); 982 } 983 984 /* 985 Internal implementation of `popFront()`, which 986 can be used with an arbitrary `State` instance 987 */ 988 private static void popFrontImpl(ref State mtState) @nogc 989 { 990 /* This function blends two nominally independent 991 processes: (i) calculation of the next random 992 variate `mtState.front` from the cached previous 993 `data` entry `z`, and (ii) updating the value 994 of `data[index]` and `mtState.z` and advancing 995 the `index` value to the next in sequence. 996 997 By interweaving the steps involved in these 998 procedures, rather than performing each of 999 them separately in sequence, the variables 1000 are kept 'hot' in CPU registers, allowing 1001 for significantly faster performance. */ 1002 ptrdiff_t index = mtState.index; 1003 ptrdiff_t next = index - 1; 1004 if (next < 0) 1005 next = n - 1; 1006 auto z = mtState.z; 1007 ptrdiff_t conj = index - m; 1008 if (conj < 0) 1009 conj = index - m + n; 1010 1011 static if (d == UIntType.max) 1012 { 1013 z ^= (z >> u); 1014 } 1015 else 1016 { 1017 z ^= (z >> u) & d; 1018 } 1019 1020 auto q = mtState.data[index] & upperMask; 1021 auto p = mtState.data[next] & lowerMask; 1022 z ^= (z << s) & b; 1023 auto y = q | p; 1024 auto x = y >> 1; 1025 z ^= (z << t) & c; 1026 if (y & 1) 1027 x ^= a; 1028 auto e = mtState.data[conj] ^ x; 1029 z ^= (z >> l); 1030 mtState.z = mtState.data[index] = e; 1031 mtState.index = next; 1032 1033 /* technically we should take the lowest `w` 1034 bits here, but if the tempering bitmasks 1035 `b` and `c` are set correctly, this should 1036 be unnecessary */ 1037 mtState.front = z; 1038 } 1039 1040 /** 1041 Returns the current random value. 1042 */ 1043 @property UIntType front() @safe const pure nothrow @nogc 1044 { 1045 return this.state.front; 1046 } 1047 1048 /// 1049 @property typeof(this) save() @safe const pure nothrow @nogc 1050 { 1051 return this; 1052 } 1053 1054 /** 1055 Always `false`. 1056 */ 1057 enum bool empty = false; 1058 } 1059 1060 /// 1061 @safe unittest 1062 { 1063 // seed with a constant 1064 Mt19937 gen; 1065 auto n = gen.front; // same for each run 1066 assert(n == 3499211612); 1067 1068 // Seed with an unpredictable value 1069 gen.seed(unpredictableSeed); 1070 n = gen.front; // different across runs 1071 } 1072 1073 /** 1074 A `MersenneTwisterEngine` instantiated with the parameters of the 1075 original engine $(HTTP math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html, 1076 MT19937), generating uniformly-distributed 32-bit numbers with a 1077 period of 2 to the power of 19937. Recommended for random number 1078 generation unless memory is severely restricted, in which case a $(LREF 1079 LinearCongruentialEngine) would be the generator of choice. 1080 */ 1081 alias Mt19937 = MersenneTwisterEngine!(uint, 32, 624, 397, 31, 1082 0x9908b0df, 11, 0xffffffff, 7, 1083 0x9d2c5680, 15, 1084 0xefc60000, 18, 1_812_433_253); 1085 1086 /// 1087 @safe @nogc unittest 1088 { 1089 // seed with a constant 1090 Mt19937 gen; 1091 auto n = gen.front; // same for each run 1092 assert(n == 3499211612); 1093 1094 // Seed with an unpredictable value 1095 gen.seed(unpredictableSeed); 1096 n = gen.front; // different across runs 1097 } 1098 1099 @safe nothrow unittest 1100 { 1101 import std.algorithm; 1102 import std.range; 1103 static assert(isUniformRNG!Mt19937); 1104 static assert(isUniformRNG!(Mt19937, uint)); 1105 static assert(isSeedable!Mt19937); 1106 static assert(isSeedable!(Mt19937, uint)); 1107 static assert(isSeedable!(Mt19937, typeof(map!((a) => unpredictableSeed)(repeat(0))))); 1108 Mt19937 gen; 1109 assert(gen.front == 3499211612); 1110 popFrontN(gen, 9999); 1111 assert(gen.front == 4123659995); 1112 try { gen.seed(iota(624u)); } catch (Exception) { assert(false); } 1113 assert(gen.front == 3708921088u); 1114 popFrontN(gen, 9999); 1115 assert(gen.front == 165737292u); 1116 } 1117 1118 /** 1119 A `MersenneTwisterEngine` instantiated with the parameters of the 1120 original engine $(HTTP en.wikipedia.org/wiki/Mersenne_Twister, 1121 MT19937-64), generating uniformly-distributed 64-bit numbers with a 1122 period of 2 to the power of 19937. 1123 */ 1124 alias Mt19937_64 = MersenneTwisterEngine!(ulong, 64, 312, 156, 31, 1125 0xb5026f5aa96619e9, 29, 0x5555555555555555, 17, 1126 0x71d67fffeda60000, 37, 1127 0xfff7eee000000000, 43, 6_364_136_223_846_793_005); 1128 1129 /// 1130 @safe @nogc unittest 1131 { 1132 // Seed with a constant 1133 auto gen = Mt19937_64(12345); 1134 auto n = gen.front; // same for each run 1135 assert(n == 6597103971274460346); 1136 1137 // Seed with an unpredictable value 1138 gen.seed(unpredictableSeed!ulong); 1139 n = gen.front; // different across runs 1140 } 1141 1142 @safe nothrow unittest 1143 { 1144 import std.algorithm; 1145 import std.range; 1146 static assert(isUniformRNG!Mt19937_64); 1147 static assert(isUniformRNG!(Mt19937_64, ulong)); 1148 static assert(isSeedable!Mt19937_64); 1149 static assert(isSeedable!(Mt19937_64, ulong)); 1150 static assert(isSeedable!(Mt19937_64, typeof(map!((a) => unpredictableSeed!ulong)(repeat(0))))); 1151 Mt19937_64 gen; 1152 assert(gen.front == 14514284786278117030uL); 1153 popFrontN(gen, 9999); 1154 assert(gen.front == 9981545732273789042uL); 1155 try { gen.seed(iota(312uL)); } catch (Exception) { assert(false); } 1156 assert(gen.front == 14660652410669508483uL); 1157 popFrontN(gen, 9999); 1158 assert(gen.front == 15956361063660440239uL); 1159 } 1160 1161 @safe unittest 1162 { 1163 import std.algorithm; 1164 import std.exception; 1165 import std.range; 1166 1167 Mt19937 gen; 1168 1169 assertThrown(gen.seed(map!((a) => 123_456_789U)(repeat(0, 623)))); 1170 1171 gen.seed(123_456_789U.repeat(624)); 1172 //infinite Range 1173 gen.seed(123_456_789U.repeat); 1174 } 1175 1176 @safe @nogc pure nothrow unittest 1177 { 1178 uint a, b; 1179 { 1180 Mt19937 gen; 1181 a = gen.front; 1182 } 1183 { 1184 Mt19937 gen; 1185 gen.popFront(); 1186 //popFrontN(gen, 1); // skip 1 element 1187 b = gen.front; 1188 } 1189 assert(a != b); 1190 } 1191 1192 @safe @nogc unittest 1193 { 1194 // Check .save works 1195 static foreach (Type; std.meta.AliasSeq!(Mt19937, Mt19937_64)) 1196 {{ 1197 auto gen1 = Type(123_456_789); 1198 gen1.popFront(); 1199 // https://issues.dlang.org/show_bug.cgi?id=15853 1200 auto gen2 = ((const ref Type a) => a.save())(gen1); 1201 assert(gen1 == gen2); // Danger, Will Robinson -- no opEquals for MT 1202 // Enable next test when RNGs are reference types 1203 version (none) { assert(gen1 !is gen2); } 1204 for (auto i = 0; i < 100; i++, gen1.popFront, gen2.popFront) 1205 assert(gen1.front() == gen2.front()); 1206 }} 1207 } 1208 1209 // https://issues.dlang.org/show_bug.cgi?id=11690 1210 @safe @nogc pure nothrow unittest 1211 { 1212 alias MT(UIntType, uint w) = MersenneTwisterEngine!(UIntType, w, 624, 397, 31, 1213 0x9908b0df, 11, 0xffffffff, 7, 1214 0x9d2c5680, 15, 1215 0xefc60000, 18, 1812433253); 1216 1217 static immutable ulong[] expectedFirstValue = [3499211612uL, 3499211612uL, 1218 171143175841277uL, 1145028863177033374uL]; 1219 1220 static immutable ulong[] expected10kValue = [4123659995uL, 4123659995uL, 1221 51991688252792uL, 3031481165133029945uL]; 1222 1223 static foreach (i, R; std.meta.AliasSeq!(MT!(uint, 32), MT!(ulong, 32), MT!(ulong, 48), MT!(ulong, 64))) 1224 {{ 1225 auto a = R(); 1226 a.seed(a.defaultSeed); // checks that some alternative paths in `seed` are utilized 1227 assert(a.front == expectedFirstValue[i]); 1228 a.popFrontN(9999); 1229 assert(a.front == expected10kValue[i]); 1230 }} 1231 } 1232 1233 /++ 1234 Xorshift generator. 1235 Implemented according to $(HTTP www.jstatsoft.org/v08/i14/paper, Xorshift RNGs) 1236 (Marsaglia, 2003) when the size is small. For larger sizes the generator 1237 uses Sebastino Vigna's optimization of using an index to avoid needing 1238 to rotate the internal array. 1239 1240 Period is `2 ^^ nbits - 1` except for a legacy 192-bit uint version (see 1241 note below). 1242 1243 Params: 1244 UIntType = Word size of this xorshift generator and the return type 1245 of `opCall`. 1246 nbits = The number of bits of state of this generator. This must be 1247 a positive multiple of the size in bits of UIntType. If 1248 nbits is large this struct may occupy slightly more memory 1249 than this so it can use a circular counter instead of 1250 shifting the entire array. 1251 sa = The direction and magnitude of the 1st shift. Positive 1252 means left, negative means right. 1253 sb = The direction and magnitude of the 2nd shift. Positive 1254 means left, negative means right. 1255 sc = The direction and magnitude of the 3rd shift. Positive 1256 means left, negative means right. 1257 1258 Note: 1259 For historical compatibility when `nbits == 192` and `UIntType` is `uint` 1260 a legacy hybrid PRNG is used consisting of a 160-bit xorshift combined 1261 with a 32-bit counter. This combined generator has period equal to the 1262 least common multiple of `2^^160 - 1` and `2^^32`. 1263 1264 Previous versions of `XorshiftEngine` did not provide any mechanism to specify 1265 the directions of the shifts, taking each shift as an unsigned magnitude. 1266 For backwards compatibility, because three shifts in the same direction 1267 cannot result in a full-period XorshiftEngine, when all three of `sa`, `sb`, 1268 `sc, are positive `XorshiftEngine` treats them as unsigned magnitudes and 1269 uses shift directions to match the old behavior of `XorshiftEngine`. 1270 1271 Not every set of shifts results in a full-period xorshift generator. 1272 The template does not currently at compile-time perform a full check 1273 for maximum period but in a future version might reject parameters 1274 resulting in shorter periods. 1275 +/ 1276 struct XorshiftEngine(UIntType, uint nbits, int sa, int sb, int sc) 1277 if (isUnsigned!UIntType && !(sa > 0 && sb > 0 && sc > 0)) 1278 { 1279 static assert(nbits > 0 && nbits % (UIntType.sizeof * 8) == 0, 1280 "nbits must be an even multiple of "~UIntType.stringof 1281 ~".sizeof * 8, not "~nbits.stringof~"."); 1282 1283 static assert(!((sa >= 0) == (sb >= 0) && (sa >= 0) == (sc >= 0)) 1284 && (sa * sb * sc != 0), 1285 "shifts cannot be zero and cannot all be in same direction: cannot be [" 1286 ~sa.stringof~", "~sb.stringof~", "~sc.stringof~"]."); 1287 1288 static assert(sa != sb && sb != sc, 1289 "consecutive shifts with the same magnitude and direction would partially or completely cancel!"); 1290 1291 static assert(UIntType.sizeof == uint.sizeof || UIntType.sizeof == ulong.sizeof, 1292 "XorshiftEngine currently does not support type " ~ UIntType.sizeof 1293 ~ " because it does not have a `seed` implementation for arrays " 1294 ~ " of element types other than uint and ulong."); 1295 1296 public: 1297 ///Mark this as a Rng 1298 enum bool isUniformRandom = true; 1299 /// Always `false` (random generators are infinite ranges). 1300 enum empty = false; 1301 /// Smallest generated value. 1302 enum UIntType min = _state.length == 1 ? 1 : 0; 1303 /// Largest generated value. 1304 enum UIntType max = UIntType.max; 1305 1306 1307 private: 1308 // Legacy 192 bit uint hybrid counter/xorshift. 1309 enum bool isLegacy192Bit = UIntType.sizeof == uint.sizeof && nbits == 192; 1310 1311 // Shift magnitudes. 1312 enum a = (sa < 0 ? -sa : sa); 1313 enum b = (sb < 0 ? -sb : sb); 1314 enum c = (sc < 0 ? -sc : sc); 1315 1316 // Shift expressions to mix in. 1317 enum shiftA(string expr) = `((`~expr~`) `~(sa > 0 ? `<< a)` : ` >>> a)`); 1318 enum shiftB(string expr) = `((`~expr~`) `~(sb > 0 ? `<< b)` : ` >>> b)`); 1319 enum shiftC(string expr) = `((`~expr~`) `~(sc > 0 ? `<< c)` : ` >>> c)`); 1320 1321 enum N = nbits / (UIntType.sizeof * 8); 1322 1323 // For N <= 2 it is strictly worse to use an index. 1324 // Informal third-party benchmarks suggest that for `ulong` it is 1325 // faster to use an index when N == 4. For `uint` we err on the side 1326 // of not increasing the struct's size and only switch to the other 1327 // implementation when N > 4. 1328 enum useIndex = !isLegacy192Bit && (UIntType.sizeof >= ulong.sizeof ? N > 3 : N > 4); 1329 static if (useIndex) 1330 { 1331 enum initialIndex = N - 1; 1332 uint _index = initialIndex; 1333 } 1334 1335 static if (N == 1 && UIntType.sizeof <= uint.sizeof) 1336 { 1337 UIntType[N] _state = [cast(UIntType) 2_463_534_242]; 1338 } 1339 else static if (isLegacy192Bit) 1340 { 1341 UIntType[N] _state = [123_456_789, 362_436_069, 521_288_629, 88_675_123, 5_783_321, 6_615_241]; 1342 UIntType value_; 1343 } 1344 else static if (N <= 5 && UIntType.sizeof <= uint.sizeof) 1345 { 1346 UIntType[N] _state = [ 1347 cast(UIntType) 123_456_789, 1348 cast(UIntType) 362_436_069, 1349 cast(UIntType) 521_288_629, 1350 cast(UIntType) 88_675_123, 1351 cast(UIntType) 5_783_321][0 .. N]; 1352 } 1353 else 1354 { 1355 UIntType[N] _state = () 1356 { 1357 static if (UIntType.sizeof < ulong.sizeof) 1358 { 1359 uint x0 = 123_456_789; 1360 enum uint m = 1_812_433_253U; 1361 } 1362 else static if (UIntType.sizeof <= ulong.sizeof) 1363 { 1364 ulong x0 = 123_456_789; 1365 enum ulong m = 6_364_136_223_846_793_005UL; 1366 } 1367 else 1368 { 1369 static assert(0, "Phobos Error: Xorshift has no instantiation rule for " 1370 ~ UIntType.stringof); 1371 } 1372 enum uint rshift = typeof(x0).sizeof * 8 - 2; 1373 UIntType[N] result = void; 1374 foreach (i, ref e; result) 1375 { 1376 e = cast(UIntType) (x0 = (m * (x0 ^ (x0 >>> rshift)) + i + 1)); 1377 if (e == 0) 1378 e = cast(UIntType) (i + 1); 1379 } 1380 return result; 1381 }(); 1382 } 1383 1384 1385 public: 1386 /++ 1387 Constructs a `XorshiftEngine` generator seeded with $(D_PARAM x0). 1388 1389 Params: 1390 x0 = value used to deterministically initialize internal state 1391 +/ 1392 this()(UIntType x0) @safe pure nothrow @nogc 1393 { 1394 seed(x0); 1395 } 1396 1397 1398 /++ 1399 (Re)seeds the generator. 1400 1401 Params: 1402 x0 = value used to deterministically initialize internal state 1403 +/ 1404 void seed()(UIntType x0) @safe pure nothrow @nogc 1405 { 1406 static if (useIndex) 1407 _index = initialIndex; 1408 1409 static if (UIntType.sizeof == uint.sizeof) 1410 { 1411 // Initialization routine from MersenneTwisterEngine. 1412 foreach (uint i, ref e; _state) 1413 { 1414 e = (x0 = (1_812_433_253U * (x0 ^ (x0 >> 30)) + i + 1)); 1415 // Xorshift requires merely that not every word of the internal 1416 // array is 0. For historical compatibility the 32-bit word version 1417 // has the stronger requirement that not any word of the state 1418 // array is 0 after initial seeding. 1419 if (e == 0) 1420 e = (i + 1); 1421 } 1422 } 1423 else static if (UIntType.sizeof == ulong.sizeof) 1424 { 1425 static if (N > 1) 1426 { 1427 // Initialize array using splitmix64 as recommended by Sebastino Vigna. 1428 // By construction the array will not be all zeroes. 1429 // http://xoroshiro.di.unimi.it/splitmix64.c 1430 foreach (ref e; _state) 1431 { 1432 ulong z = (x0 += 0x9e37_79b9_7f4a_7c15UL); 1433 z = (z ^ (z >>> 30)) * 0xbf58_476d_1ce4_e5b9UL; 1434 z = (z ^ (z >>> 27)) * 0x94d0_49bb_1331_11ebUL; 1435 e = z ^ (z >>> 31); 1436 } 1437 } 1438 else 1439 { 1440 // Apply a transformation when N == 1 instead of just copying x0 1441 // directly because it's not unlikely that a user might initialize 1442 // a PRNG with small counting numbers (e.g. 1, 2, 3) that have the 1443 // statistically rare property of having only 1 or 2 non-zero bits. 1444 // Additionally we need to ensure that the internal state is not 1445 // entirely zero. 1446 if (x0 != 0) 1447 _state[0] = x0 * 6_364_136_223_846_793_005UL; 1448 else 1449 _state[0] = typeof(this).init._state[0]; 1450 } 1451 } 1452 else 1453 { 1454 static assert(0, "Phobos Error: Xorshift has no `seed` implementation for " 1455 ~ UIntType.stringof); 1456 } 1457 1458 popFront(); 1459 } 1460 1461 1462 /** 1463 * Returns the current number in the random sequence. 1464 */ 1465 @property 1466 UIntType front() const @safe pure nothrow @nogc 1467 { 1468 static if (isLegacy192Bit) 1469 return value_; 1470 else static if (!useIndex) 1471 return _state[N-1]; 1472 else 1473 return _state[_index]; 1474 } 1475 1476 1477 /** 1478 * Advances the random sequence. 1479 */ 1480 void popFront() @safe pure nothrow @nogc 1481 { 1482 alias s = _state; 1483 static if (isLegacy192Bit) 1484 { 1485 auto x = _state[0] ^ mixin(shiftA!`s[0]`); 1486 static foreach (i; 0 .. N-2) 1487 s[i] = s[i + 1]; 1488 s[N-2] = s[N-2] ^ mixin(shiftC!`s[N-2]`) ^ x ^ mixin(shiftB!`x`); 1489 value_ = s[N-2] + (s[N-1] += 362_437); 1490 } 1491 else static if (N == 1) 1492 { 1493 s[0] ^= mixin(shiftA!`s[0]`); 1494 s[0] ^= mixin(shiftB!`s[0]`); 1495 s[0] ^= mixin(shiftC!`s[0]`); 1496 } 1497 else static if (!useIndex) 1498 { 1499 auto x = s[0] ^ mixin(shiftA!`s[0]`); 1500 static foreach (i; 0 .. N-1) 1501 s[i] = s[i + 1]; 1502 s[N-1] = s[N-1] ^ mixin(shiftC!`s[N-1]`) ^ x ^ mixin(shiftB!`x`); 1503 } 1504 else 1505 { 1506 assert(_index < N); // Invariant. 1507 const sIndexMinus1 = s[_index]; 1508 static if ((N & (N - 1)) == 0) 1509 _index = (_index + 1) & typeof(_index)(N - 1); 1510 else 1511 { 1512 if (++_index >= N) 1513 _index = 0; 1514 } 1515 auto x = s[_index]; 1516 x ^= mixin(shiftA!`x`); 1517 s[_index] = sIndexMinus1 ^ mixin(shiftC!`sIndexMinus1`) ^ x ^ mixin(shiftB!`x`); 1518 } 1519 } 1520 1521 1522 /** 1523 * Captures a range state. 1524 */ 1525 @property 1526 typeof(this) save() const @safe pure nothrow @nogc 1527 { 1528 return this; 1529 } 1530 1531 private: 1532 // Workaround for a DScanner bug. If we remove this `private:` DScanner 1533 // gives erroneous warnings about missing documentation for public symbols 1534 // later in the module. 1535 } 1536 1537 /// ditto 1538 template XorshiftEngine(UIntType, int bits, int a, int b, int c) 1539 if (isUnsigned!UIntType && a > 0 && b > 0 && c > 0) 1540 { 1541 // Compatibility with old parameterizations without explicit shift directions. 1542 static if (bits == UIntType.sizeof * 8) 1543 alias XorshiftEngine = .XorshiftEngine!(UIntType, bits, a, -b, c);//left, right, left 1544 else static if (bits == 192 && UIntType.sizeof == uint.sizeof) 1545 alias XorshiftEngine = .XorshiftEngine!(UIntType, bits, -a, b, c);//right, left, left 1546 else 1547 alias XorshiftEngine = .XorshiftEngine!(UIntType, bits, a, -b, -c);//left, right, right 1548 } 1549 1550 /// 1551 @safe unittest 1552 { 1553 alias Xorshift96 = XorshiftEngine!(uint, 96, 10, 5, 26); 1554 auto rnd = Xorshift96(42); 1555 auto num = rnd.front; // same for each run 1556 assert(num == 2704588748); 1557 } 1558 1559 1560 /** 1561 * Define `XorshiftEngine` generators with well-chosen parameters. See each bits examples of "Xorshift RNGs". 1562 * `Xorshift` is a Xorshift128's alias because 128bits implementation is mostly used. 1563 */ 1564 alias Xorshift32 = XorshiftEngine!(uint, 32, 13, 17, 15) ; 1565 alias Xorshift64 = XorshiftEngine!(uint, 64, 10, 13, 10); /// ditto 1566 alias Xorshift96 = XorshiftEngine!(uint, 96, 10, 5, 26); /// ditto 1567 alias Xorshift128 = XorshiftEngine!(uint, 128, 11, 8, 19); /// ditto 1568 alias Xorshift160 = XorshiftEngine!(uint, 160, 2, 1, 4); /// ditto 1569 alias Xorshift192 = XorshiftEngine!(uint, 192, 2, 1, 4); /// ditto 1570 alias Xorshift = Xorshift128; /// ditto 1571 1572 /// 1573 @safe @nogc unittest 1574 { 1575 // Seed with a constant 1576 auto rnd = Xorshift(1); 1577 auto num = rnd.front; // same for each run 1578 assert(num == 1405313047); 1579 1580 // Seed with an unpredictable value 1581 rnd.seed(unpredictableSeed); 1582 num = rnd.front; // different across rnd 1583 } 1584 1585 @safe @nogc unittest 1586 { 1587 import std.range; 1588 static assert(isForwardRange!Xorshift); 1589 static assert(isUniformRNG!Xorshift); 1590 static assert(isUniformRNG!(Xorshift, uint)); 1591 static assert(isSeedable!Xorshift); 1592 static assert(isSeedable!(Xorshift, uint)); 1593 1594 static assert(Xorshift32.min == 1); 1595 1596 // Result from reference implementation. 1597 static ulong[][] checking = [ 1598 [2463534242UL, 901999875, 3371835698, 2675058524, 1053936272, 3811264849, 1599 472493137, 3856898176, 2131710969, 2312157505], 1600 [362436069UL, 2113136921, 19051112, 3010520417, 951284840, 1213972223, 1601 3173832558, 2611145638, 2515869689, 2245824891], 1602 [521288629UL, 1950277231, 185954712, 1582725458, 3580567609, 2303633688, 1603 2394948066, 4108622809, 1116800180, 3357585673], 1604 [88675123UL, 3701687786, 458299110, 2500872618, 3633119408, 516391518, 1605 2377269574, 2599949379, 717229868, 137866584], 1606 [5783321UL, 393427209, 1947109840, 565829276, 1006220149, 971147905, 1607 1436324242, 2800460115, 1484058076, 3823330032], 1608 [0UL, 246875399, 3690007200, 1264581005, 3906711041, 1866187943, 2481925219, 1609 2464530826, 1604040631, 3653403911], 1610 [16749904790159980466UL, 14489774923612894650UL, 148813570191443371UL, 1611 6529783008780612886UL, 10182425759614080046UL, 16549659571055687844UL, 1612 542957868271744939UL, 9459127085596028450UL, 16001049981702441780UL, 1613 7351634712593111741], 1614 [14750058843113249055UL, 17731577314455387619UL, 1314705253499959044UL, 1615 3113030620614841056UL, 9468075444678629182UL, 13962152036600088141UL, 1616 9030205609946043947UL, 1856726150434672917UL, 8098922200110395314UL, 1617 2772699174618556175UL], 1618 ]; 1619 1620 alias XorshiftTypes = std.meta.AliasSeq!(Xorshift32, Xorshift64, Xorshift96, 1621 Xorshift128, Xorshift160, Xorshift192, Xorshift64_64, Xorshift128_64); 1622 1623 foreach (I, Type; XorshiftTypes) 1624 { 1625 Type rnd; 1626 1627 foreach (e; checking[I]) 1628 { 1629 assert(rnd.front == e); 1630 rnd.popFront(); 1631 } 1632 } 1633 1634 // Check .save works 1635 foreach (Type; XorshiftTypes) 1636 { 1637 auto rnd1 = Type(123_456_789); 1638 rnd1.popFront(); 1639 // https://issues.dlang.org/show_bug.cgi?id=15853 1640 auto rnd2 = ((const ref Type a) => a.save())(rnd1); 1641 assert(rnd1 == rnd2); 1642 // Enable next test when RNGs are reference types 1643 version (none) { assert(rnd1 !is rnd2); } 1644 for (auto i = 0; i <= Type.sizeof / 4; i++, rnd1.popFront, rnd2.popFront) 1645 assert(rnd1.front() == rnd2.front()); 1646 } 1647 } 1648 1649 1650 /* A complete list of all pseudo-random number generators implemented in 1651 * std.random. This can be used to confirm that a given function or 1652 * object is compatible with all the pseudo-random number generators 1653 * available. It is enabled only in unittest mode. 1654 */ 1655 @safe @nogc unittest 1656 { 1657 foreach (Rng; PseudoRngTypes) 1658 { 1659 static assert(isUniformRNG!Rng); 1660 auto rng = Rng(123_456_789); 1661 } 1662 } 1663 1664 version (CRuntime_Bionic) 1665 version = SecureARC4Random; // ChaCha20 1666 version (Darwin) 1667 version = SecureARC4Random; // AES 1668 version (OpenBSD) 1669 version = SecureARC4Random; // ChaCha20 1670 version (NetBSD) 1671 version = SecureARC4Random; // ChaCha20 1672 1673 version (CRuntime_UClibc) 1674 version = LegacyARC4Random; // ARC4 1675 version (FreeBSD) 1676 version = LegacyARC4Random; // ARC4 1677 version (DragonFlyBSD) 1678 version = LegacyARC4Random; // ARC4 1679 version (BSD) 1680 version = LegacyARC4Random; // Unknown implementation 1681 1682 // For the current purpose of unpredictableSeed the difference between 1683 // a secure arc4random implementation and a legacy implementation is 1684 // unimportant. The source code documents this distinction in case in the 1685 // future Phobos is altered to require cryptographically secure sources 1686 // of randomness, and also so other people reading this source code (as 1687 // Phobos is often looked to as an example of good D programming practices) 1688 // do not mistakenly use insecure versions of arc4random in contexts where 1689 // cryptographically secure sources of randomness are needed. 1690 1691 // Performance note: ChaCha20 is about 70% faster than ARC4, contrary to 1692 // what one might assume from it being more secure. 1693 1694 version (SecureARC4Random) 1695 version = AnyARC4Random; 1696 version (LegacyARC4Random) 1697 version = AnyARC4Random; 1698 1699 version (AnyARC4Random) 1700 { 1701 extern(C) private @nogc nothrow 1702 { 1703 uint arc4random() @safe; 1704 void arc4random_buf(scope void* buf, size_t nbytes) @system; 1705 } 1706 } 1707 else 1708 { 1709 private ulong bootstrapSeed() @nogc nothrow 1710 { 1711 // https://issues.dlang.org/show_bug.cgi?id=19580 1712 // previously used `ulong result = void` to start with an arbitary value 1713 // but using an uninitialized variable's value is undefined behavior 1714 // and enabled unwanted optimizations on non-DMD compilers. 1715 ulong result; 1716 enum ulong m = 0xc6a4_a793_5bd1_e995UL; // MurmurHash2_64A constant. 1717 void updateResult(ulong x) 1718 { 1719 x *= m; 1720 x = (x ^ (x >>> 47)) * m; 1721 result = (result ^ x) * m; 1722 } 1723 import core.thread : getpid, Thread; 1724 import core.time : MonoTime; 1725 1726 updateResult(cast(ulong) cast(void*) Thread.getThis()); 1727 updateResult(cast(ulong) getpid()); 1728 updateResult(cast(ulong) MonoTime.currTime.ticks); 1729 result = (result ^ (result >>> 47)) * m; 1730 return result ^ (result >>> 47); 1731 } 1732 1733 // If we don't have arc4random and we don't have RDRAND fall back to this. 1734 private ulong fallbackSeed() @nogc nothrow 1735 { 1736 // Bit avalanche function from MurmurHash3. 1737 static ulong fmix64(ulong k) @nogc nothrow pure @safe 1738 { 1739 k = (k ^ (k >>> 33)) * 0xff51afd7ed558ccd; 1740 k = (k ^ (k >>> 33)) * 0xc4ceb9fe1a85ec53; 1741 return k ^ (k >>> 33); 1742 } 1743 // Using SplitMix algorithm with constant gamma. 1744 // Chosen gamma is the odd number closest to 2^^64 1745 // divided by the silver ratio (1.0L + sqrt(2.0L)). 1746 enum gamma = 0x6a09e667f3bcc909UL; 1747 import core.atomic : has64BitCAS; 1748 static if (has64BitCAS) 1749 { 1750 import core.atomic : MemoryOrder, atomicLoad, atomicOp, atomicStore, cas; 1751 shared static ulong seed; 1752 shared static bool initialized; 1753 if (0 == atomicLoad!(MemoryOrder.raw)(initialized)) 1754 { 1755 cas(&seed, 0UL, fmix64(bootstrapSeed())); 1756 atomicStore!(MemoryOrder.rel)(initialized, true); 1757 } 1758 return fmix64(atomicOp!"+="(seed, gamma)); 1759 } 1760 else 1761 { 1762 static ulong seed; 1763 static bool initialized; 1764 if (!initialized) 1765 { 1766 seed = fmix64(bootstrapSeed()); 1767 initialized = true; 1768 } 1769 return fmix64(seed += gamma); 1770 } 1771 } 1772 } 1773 1774 /** 1775 A "good" seed for initializing random number engines. Initializing 1776 with $(D_PARAM unpredictableSeed) makes engines generate different 1777 random number sequences every run. 1778 1779 Returns: 1780 A single unsigned integer seed value, different on each successive call 1781 Note: 1782 In general periodically 'reseeding' a PRNG does not improve its quality 1783 and in some cases may harm it. For an extreme example the Mersenne 1784 Twister has `2 ^^ 19937 - 1` distinct states but after `seed(uint)` is 1785 called it can only be in one of `2 ^^ 32` distinct states regardless of 1786 how excellent the source of entropy is. 1787 */ 1788 @property uint unpredictableSeed() @trusted nothrow @nogc 1789 { 1790 version (AnyARC4Random) 1791 { 1792 return arc4random(); 1793 } 1794 else 1795 { 1796 version (InlineAsm_X86_Any) 1797 { 1798 import core.cpuid : hasRdrand; 1799 if (hasRdrand) 1800 { 1801 uint result; 1802 asm @nogc nothrow 1803 { 1804 db 0x0f, 0xc7, 0xf0; // rdrand EAX 1805 jnc LnotUsingRdrand; 1806 mov result, EAX; 1807 // Some AMD CPUs shipped with bugs where RDRAND could fail 1808 // but still set the carry flag to 1. In those cases the 1809 // output will be -1. 1810 cmp EAX, 0xffff_ffff; 1811 jne LusingRdrand; 1812 // If result was -1 verify RDAND isn't constantly returning -1. 1813 db 0x0f, 0xc7, 0xf0; 1814 jnc LusingRdrand; 1815 cmp EAX, 0xffff_ffff; 1816 je LnotUsingRdrand; 1817 } 1818 LusingRdrand: 1819 return result; 1820 } 1821 LnotUsingRdrand: 1822 } 1823 return cast(uint) fallbackSeed(); 1824 } 1825 } 1826 1827 /// ditto 1828 template unpredictableSeed(UIntType) 1829 if (isUnsigned!UIntType) 1830 { 1831 static if (is(UIntType == uint)) 1832 alias unpredictableSeed = .unpredictableSeed; 1833 else static if (!is(Unqual!UIntType == UIntType)) 1834 alias unpredictableSeed = .unpredictableSeed!(Unqual!UIntType); 1835 else 1836 /// ditto 1837 @property UIntType unpredictableSeed() @nogc nothrow @trusted 1838 { 1839 version (AnyARC4Random) 1840 { 1841 static if (UIntType.sizeof <= uint.sizeof) 1842 { 1843 return cast(UIntType) arc4random(); 1844 } 1845 else 1846 { 1847 UIntType result = void; 1848 arc4random_buf(&result, UIntType.sizeof); 1849 return result; 1850 } 1851 } 1852 else 1853 { 1854 version (InlineAsm_X86_Any) 1855 { 1856 import core.cpuid : hasRdrand; 1857 if (hasRdrand) 1858 { 1859 static if (UIntType.sizeof <= uint.sizeof) 1860 { 1861 uint result; 1862 asm @nogc nothrow 1863 { 1864 db 0x0f, 0xc7, 0xf0; // rdrand EAX 1865 jnc LnotUsingRdrand; 1866 mov result, EAX; 1867 // Some AMD CPUs shipped with bugs where RDRAND could fail 1868 // but still set the carry flag to 1. In those cases the 1869 // output will be -1. 1870 cmp EAX, 0xffff_ffff; 1871 jne LusingRdrand; 1872 // If result was -1 verify RDAND isn't constantly returning -1. 1873 db 0x0f, 0xc7, 0xf0; 1874 jnc LusingRdrand; 1875 cmp EAX, 0xffff_ffff; 1876 je LnotUsingRdrand; 1877 } 1878 LusingRdrand: 1879 return cast(UIntType) result; 1880 } 1881 else version (D_InlineAsm_X86_64) 1882 { 1883 ulong result; 1884 asm @nogc nothrow 1885 { 1886 db 0x48, 0x0f, 0xc7, 0xf0; // rdrand RAX 1887 jnc LnotUsingRdrand; 1888 mov result, RAX; 1889 // Some AMD CPUs shipped with bugs where RDRAND could fail 1890 // but still set the carry flag to 1. In those cases the 1891 // output will be -1. 1892 cmp RAX, 0xffff_ffff_ffff_ffff; 1893 jne LusingRdrand; 1894 // If result was -1 verify RDAND isn't constantly returning -1. 1895 db 0x48, 0x0f, 0xc7, 0xf0; 1896 jnc LusingRdrand; 1897 cmp RAX, 0xffff_ffff_ffff_ffff; 1898 je LnotUsingRdrand; 1899 } 1900 LusingRdrand: 1901 return result; 1902 } 1903 else 1904 { 1905 uint resultLow, resultHigh; 1906 asm @nogc nothrow 1907 { 1908 db 0x0f, 0xc7, 0xf0; // rdrand EAX 1909 jnc LnotUsingRdrand; 1910 mov resultLow, EAX; 1911 db 0x0f, 0xc7, 0xf0; // rdrand EAX 1912 jnc LnotUsingRdrand; 1913 mov resultHigh, EAX; 1914 } 1915 if (resultLow != uint.max || resultHigh != uint.max) // Protect against AMD RDRAND bug. 1916 return ((cast(ulong) resultHigh) << 32) ^ resultLow; 1917 } 1918 } 1919 LnotUsingRdrand: 1920 } 1921 return cast(UIntType) fallbackSeed(); 1922 } 1923 } 1924 } 1925 1926 /// 1927 @safe @nogc unittest 1928 { 1929 auto rnd = Random(unpredictableSeed); 1930 auto n = rnd.front; 1931 static assert(is(typeof(n) == uint)); 1932 } 1933 1934 /** 1935 The "default", "favorite", "suggested" random number generator type on 1936 the current platform. It is an alias for one of the previously-defined 1937 generators. You may want to use it if (1) you need to generate some 1938 nice random numbers, and (2) you don't care for the minutiae of the 1939 method being used. 1940 */ 1941 1942 alias Random = Mt19937; 1943 1944 @safe @nogc unittest 1945 { 1946 static assert(isUniformRNG!Random); 1947 static assert(isUniformRNG!(Random, uint)); 1948 static assert(isSeedable!Random); 1949 static assert(isSeedable!(Random, uint)); 1950 } 1951 1952 /** 1953 Global random number generator used by various functions in this 1954 module whenever no generator is specified. It is allocated per-thread 1955 and initialized to an unpredictable value for each thread. 1956 1957 Returns: 1958 A singleton instance of the default random number generator 1959 */ 1960 @property ref Random rndGen() @safe nothrow @nogc 1961 { 1962 static Random result; 1963 static bool initialized; 1964 if (!initialized) 1965 { 1966 static if (isSeedable!(Random, ulong)) 1967 result.seed(unpredictableSeed!ulong); // Avoid unnecessary copy. 1968 else static if (is(Random : MersenneTwisterEngine!Params, Params...)) 1969 initMTEngine(result); 1970 else static if (isSeedable!(Random, uint)) 1971 result.seed(unpredictableSeed!uint); // Avoid unnecessary copy. 1972 else 1973 result = Random(unpredictableSeed); 1974 initialized = true; 1975 } 1976 return result; 1977 } 1978 1979 /// 1980 @safe nothrow @nogc unittest 1981 { 1982 import std.algorithm.iteration : sum; 1983 import std.range : take; 1984 auto rnd = rndGen; 1985 assert(rnd.take(3).sum > 0); 1986 } 1987 1988 /+ 1989 Initialize a 32-bit MersenneTwisterEngine from 64 bits of entropy. 1990 This is private and accepts no seed as a parameter, freeing the internal 1991 implementaton from any need for stability across releases. 1992 +/ 1993 private void initMTEngine(MTEngine)(scope ref MTEngine mt) 1994 if (is(MTEngine : MersenneTwisterEngine!Params, Params...)) 1995 { 1996 pragma(inline, false); // Called no more than once per thread by rndGen. 1997 ulong seed = unpredictableSeed!ulong; 1998 static if (is(typeof(mt.seed(seed)))) 1999 { 2000 mt.seed(seed); 2001 } 2002 else 2003 { 2004 alias UIntType = typeof(mt.front()); 2005 if (seed == 0) seed = -1; // Any number but 0 is fine. 2006 uint s0 = cast(uint) seed; 2007 uint s1 = cast(uint) (seed >> 32); 2008 foreach (ref e; mt.state.data) 2009 { 2010 //http://xoshiro.di.unimi.it/xoroshiro64starstar.c 2011 const tmp = s0 * 0x9E3779BB; 2012 e = ((tmp << 5) | (tmp >> (32 - 5))) * 5; 2013 static if (MTEngine.max != UIntType.max) { e &= MTEngine.max; } 2014 2015 const tmp1 = s0 ^ s1; 2016 s0 = ((s0 << 26) | (s0 >> (32 - 26))) ^ tmp1 ^ (tmp1 << 9); 2017 s1 = (tmp1 << 13) | (tmp1 >> (32 - 13)); 2018 } 2019 2020 mt.state.index = mt.state.data.length - 1; 2021 // double popFront() to guarantee both `mt.state.z` 2022 // and `mt.state.front` are derived from the newly 2023 // set values in `mt.state.data`. 2024 mt.popFront(); 2025 mt.popFront(); 2026 } 2027 } 2028 2029 /** 2030 Generates a number between `a` and `b`. The `boundaries` 2031 parameter controls the shape of the interval (open vs. closed on 2032 either side). Valid values for `boundaries` are `"[]"`, $(D 2033 "$(LPAREN)]"), `"[$(RPAREN)"`, and `"()"`. The default interval 2034 is closed to the left and open to the right. The version that does not 2035 take `urng` uses the default generator `rndGen`. 2036 2037 Params: 2038 a = lower bound of the _uniform distribution 2039 b = upper bound of the _uniform distribution 2040 urng = (optional) random number generator to use; 2041 if not specified, defaults to `rndGen` 2042 2043 Returns: 2044 A single random variate drawn from the _uniform distribution 2045 between `a` and `b`, whose type is the common type of 2046 these parameters 2047 */ 2048 auto uniform(string boundaries = "[)", T1, T2) 2049 (T1 a, T2 b) 2050 if (!is(CommonType!(T1, T2) == void)) 2051 { 2052 return uniform!(boundaries, T1, T2, Random)(a, b, rndGen); 2053 } 2054 2055 /// 2056 @safe unittest 2057 { 2058 auto rnd = Random(unpredictableSeed); 2059 2060 // Generate an integer in [0, 1023] 2061 auto a = uniform(0, 1024, rnd); 2062 assert(0 <= a && a < 1024); 2063 2064 // Generate a float in [0, 1) 2065 auto b = uniform(0.0f, 1.0f, rnd); 2066 assert(0 <= b && b < 1); 2067 2068 // Generate a float in [0, 1] 2069 b = uniform!"[]"(0.0f, 1.0f, rnd); 2070 assert(0 <= b && b <= 1); 2071 2072 // Generate a float in (0, 1) 2073 b = uniform!"()"(0.0f, 1.0f, rnd); 2074 assert(0 < b && b < 1); 2075 } 2076 2077 /// Create an array of random numbers using range functions and UFCS 2078 @safe unittest 2079 { 2080 import std.array : array; 2081 import std.range : generate, takeExactly; 2082 2083 int[] arr = generate!(() => uniform(0, 100)).takeExactly(10).array; 2084 assert(arr.length == 10); 2085 assert(arr[0] >= 0 && arr[0] < 100); 2086 } 2087 2088 @safe unittest 2089 { 2090 MinstdRand0 gen; 2091 foreach (i; 0 .. 20) 2092 { 2093 auto x = uniform(0.0, 15.0, gen); 2094 assert(0 <= x && x < 15); 2095 } 2096 foreach (i; 0 .. 20) 2097 { 2098 auto x = uniform!"[]"('a', 'z', gen); 2099 assert('a' <= x && x <= 'z'); 2100 } 2101 2102 foreach (i; 0 .. 20) 2103 { 2104 auto x = uniform('a', 'z', gen); 2105 assert('a' <= x && x < 'z'); 2106 } 2107 2108 foreach (i; 0 .. 20) 2109 { 2110 immutable ubyte a = 0; 2111 immutable ubyte b = 15; 2112 auto x = uniform(a, b, gen); 2113 assert(a <= x && x < b); 2114 } 2115 } 2116 2117 // Implementation of uniform for floating-point types 2118 /// ditto 2119 auto uniform(string boundaries = "[)", 2120 T1, T2, UniformRandomNumberGenerator) 2121 (T1 a, T2 b, ref UniformRandomNumberGenerator urng) 2122 if (isFloatingPoint!(CommonType!(T1, T2)) && isUniformRNG!UniformRandomNumberGenerator) 2123 { 2124 import std.conv : text; 2125 import std.exception : enforce; 2126 alias NumberType = Unqual!(CommonType!(T1, T2)); 2127 static if (boundaries[0] == '(') 2128 { 2129 import std.math.operations : nextafter; 2130 NumberType _a = nextafter(cast(NumberType) a, NumberType.infinity); 2131 } 2132 else 2133 { 2134 NumberType _a = a; 2135 } 2136 static if (boundaries[1] == ')') 2137 { 2138 import std.math.operations : nextafter; 2139 NumberType _b = nextafter(cast(NumberType) b, -NumberType.infinity); 2140 } 2141 else 2142 { 2143 NumberType _b = b; 2144 } 2145 enforce(_a <= _b, 2146 text("std.random.uniform(): invalid bounding interval ", 2147 boundaries[0], a, ", ", b, boundaries[1])); 2148 NumberType result = 2149 _a + (_b - _a) * cast(NumberType) (urng.front - urng.min) 2150 / (urng.max - urng.min); 2151 urng.popFront(); 2152 return result; 2153 } 2154 2155 // Implementation of uniform for integral types 2156 /+ Description of algorithm and suggestion of correctness: 2157 2158 The modulus operator maps an integer to a small, finite space. For instance, `x 2159 % 3` will map whatever x is into the range [0 .. 3). 0 maps to 0, 1 maps to 1, 2 2160 maps to 2, 3 maps to 0, and so on infinitely. As long as the integer is 2161 uniformly chosen from the infinite space of all non-negative integers then `x % 2162 3` will uniformly fall into that range. 2163 2164 (Non-negative is important in this case because some definitions of modulus, 2165 namely the one used in computers generally, map negative numbers differently to 2166 (-3 .. 0]. `uniform` does not use negative number modulus, thus we can safely 2167 ignore that fact.) 2168 2169 The issue with computers is that integers have a finite space they must fit in, 2170 and our uniformly chosen random number is picked in that finite space. So, that 2171 method is not sufficient. You can look at it as the integer space being divided 2172 into "buckets" and every bucket after the first bucket maps directly into that 2173 first bucket. `[0, 1, 2]`, `[3, 4, 5]`, ... When integers are finite, then the 2174 last bucket has the chance to be "incomplete": `[uint.max - 3, uint.max - 2, 2175 uint.max - 1]`, `[uint.max]` ... (the last bucket only has 1!). The issue here 2176 is that _every_ bucket maps _completely_ to the first bucket except for that 2177 last one. The last one doesn't have corresponding mappings to 1 or 2, in this 2178 case, which makes it unfair. 2179 2180 So, the answer is to simply "reroll" if you're in that last bucket, since it's 2181 the only unfair one. Eventually you'll roll into a fair bucket. Simply, instead 2182 of the meaning of the last bucket being "maps to `[0]`", it changes to "maps to 2183 `[0, 1, 2]`", which is precisely what we want. 2184 2185 To generalize, `upperDist` represents the size of our buckets (and, thus, the 2186 exclusive upper bound for our desired uniform number). `rnum` is a uniformly 2187 random number picked from the space of integers that a computer can hold (we'll 2188 say `UpperType` represents that type). 2189 2190 We'll first try to do the mapping into the first bucket by doing `offset = rnum 2191 % upperDist`. We can figure out the position of the front of the bucket we're in 2192 by `bucketFront = rnum - offset`. 2193 2194 If we start at `UpperType.max` and walk backwards `upperDist - 1` spaces, then 2195 the space we land on is the last acceptable position where a full bucket can 2196 fit: 2197 2198 --- 2199 bucketFront UpperType.max 2200 v v 2201 [..., 0, 1, 2, ..., upperDist - 1] 2202 ^~~ upperDist - 1 ~~^ 2203 --- 2204 2205 If the bucket starts any later, then it must have lost at least one number and 2206 at least that number won't be represented fairly. 2207 2208 --- 2209 bucketFront UpperType.max 2210 v v 2211 [..., upperDist - 1, 0, 1, 2, ..., upperDist - 2] 2212 ^~~~~~~~ upperDist - 1 ~~~~~~~^ 2213 --- 2214 2215 Hence, our condition to reroll is 2216 `bucketFront > (UpperType.max - (upperDist - 1))` 2217 +/ 2218 auto uniform(string boundaries = "[)", T1, T2, RandomGen) 2219 (T1 a, T2 b, ref RandomGen rng) 2220 if ((isIntegral!(CommonType!(T1, T2)) || isSomeChar!(CommonType!(T1, T2))) && 2221 isUniformRNG!RandomGen) 2222 { 2223 import std.conv : text, unsigned; 2224 import std.exception : enforce; 2225 alias ResultType = Unqual!(CommonType!(T1, T2)); 2226 static if (boundaries[0] == '(') 2227 { 2228 enforce(a < ResultType.max, 2229 text("std.random.uniform(): invalid left bound ", a)); 2230 ResultType lower = cast(ResultType) (a + 1); 2231 } 2232 else 2233 { 2234 ResultType lower = a; 2235 } 2236 2237 static if (boundaries[1] == ']') 2238 { 2239 enforce(lower <= b, 2240 text("std.random.uniform(): invalid bounding interval ", 2241 boundaries[0], a, ", ", b, boundaries[1])); 2242 /* Cannot use this next optimization with dchar, as dchar 2243 * only partially uses its full bit range 2244 */ 2245 static if (!is(ResultType == dchar)) 2246 { 2247 if (b == ResultType.max && lower == ResultType.min) 2248 { 2249 // Special case - all bits are occupied 2250 return std.random.uniform!ResultType(rng); 2251 } 2252 } 2253 auto upperDist = unsigned(b - lower) + 1u; 2254 } 2255 else 2256 { 2257 enforce(lower < b, 2258 text("std.random.uniform(): invalid bounding interval ", 2259 boundaries[0], a, ", ", b, boundaries[1])); 2260 auto upperDist = unsigned(b - lower); 2261 } 2262 2263 assert(upperDist != 0); 2264 2265 alias UpperType = typeof(upperDist); 2266 static assert(UpperType.min == 0); 2267 2268 UpperType offset, rnum, bucketFront; 2269 do 2270 { 2271 rnum = uniform!UpperType(rng); 2272 offset = rnum % upperDist; 2273 bucketFront = rnum - offset; 2274 } // while we're in an unfair bucket... 2275 while (bucketFront > (UpperType.max - (upperDist - 1))); 2276 2277 return cast(ResultType)(lower + offset); 2278 } 2279 2280 @safe unittest 2281 { 2282 import std.conv : to; 2283 auto gen = Mt19937(123_456_789); 2284 static assert(isForwardRange!(typeof(gen))); 2285 2286 auto a = uniform(0, 1024, gen); 2287 assert(0 <= a && a <= 1024); 2288 auto b = uniform(0.0f, 1.0f, gen); 2289 assert(0 <= b && b < 1, to!string(b)); 2290 auto c = uniform(0.0, 1.0); 2291 assert(0 <= c && c < 1); 2292 2293 static foreach (T; std.meta.AliasSeq!(char, wchar, dchar, byte, ubyte, short, ushort, 2294 int, uint, long, ulong, float, double, real)) 2295 {{ 2296 T lo = 0, hi = 100; 2297 2298 // Try tests with each of the possible bounds 2299 { 2300 T init = uniform(lo, hi); 2301 size_t i = 50; 2302 while (--i && uniform(lo, hi) == init) {} 2303 assert(i > 0); 2304 } 2305 { 2306 T init = uniform!"[)"(lo, hi); 2307 size_t i = 50; 2308 while (--i && uniform(lo, hi) == init) {} 2309 assert(i > 0); 2310 } 2311 { 2312 T init = uniform!"(]"(lo, hi); 2313 size_t i = 50; 2314 while (--i && uniform(lo, hi) == init) {} 2315 assert(i > 0); 2316 } 2317 { 2318 T init = uniform!"()"(lo, hi); 2319 size_t i = 50; 2320 while (--i && uniform(lo, hi) == init) {} 2321 assert(i > 0); 2322 } 2323 { 2324 T init = uniform!"[]"(lo, hi); 2325 size_t i = 50; 2326 while (--i && uniform(lo, hi) == init) {} 2327 assert(i > 0); 2328 } 2329 2330 /* Test case with closed boundaries covering whole range 2331 * of integral type 2332 */ 2333 static if (isIntegral!T || isSomeChar!T) 2334 { 2335 foreach (immutable _; 0 .. 100) 2336 { 2337 auto u = uniform!"[]"(T.min, T.max); 2338 static assert(is(typeof(u) == T)); 2339 assert(T.min <= u, "Lower bound violation for uniform!\"[]\" with " ~ T.stringof); 2340 assert(u <= T.max, "Upper bound violation for uniform!\"[]\" with " ~ T.stringof); 2341 } 2342 } 2343 }} 2344 2345 auto reproRng = Xorshift(239842); 2346 2347 static foreach (T; std.meta.AliasSeq!(char, wchar, dchar, byte, ubyte, short, 2348 ushort, int, uint, long, ulong)) 2349 {{ 2350 T lo = T.min + 10, hi = T.max - 10; 2351 T init = uniform(lo, hi, reproRng); 2352 size_t i = 50; 2353 while (--i && uniform(lo, hi, reproRng) == init) {} 2354 assert(i > 0); 2355 }} 2356 2357 { 2358 bool sawLB = false, sawUB = false; 2359 foreach (i; 0 .. 50) 2360 { 2361 auto x = uniform!"[]"('a', 'd', reproRng); 2362 if (x == 'a') sawLB = true; 2363 if (x == 'd') sawUB = true; 2364 assert('a' <= x && x <= 'd'); 2365 } 2366 assert(sawLB && sawUB); 2367 } 2368 2369 { 2370 bool sawLB = false, sawUB = false; 2371 foreach (i; 0 .. 50) 2372 { 2373 auto x = uniform('a', 'd', reproRng); 2374 if (x == 'a') sawLB = true; 2375 if (x == 'c') sawUB = true; 2376 assert('a' <= x && x < 'd'); 2377 } 2378 assert(sawLB && sawUB); 2379 } 2380 2381 { 2382 bool sawLB = false, sawUB = false; 2383 foreach (i; 0 .. 50) 2384 { 2385 immutable int lo = -2, hi = 2; 2386 auto x = uniform!"()"(lo, hi, reproRng); 2387 if (x == (lo+1)) sawLB = true; 2388 if (x == (hi-1)) sawUB = true; 2389 assert(lo < x && x < hi); 2390 } 2391 assert(sawLB && sawUB); 2392 } 2393 2394 { 2395 bool sawLB = false, sawUB = false; 2396 foreach (i; 0 .. 50) 2397 { 2398 immutable ubyte lo = 0, hi = 5; 2399 auto x = uniform(lo, hi, reproRng); 2400 if (x == lo) sawLB = true; 2401 if (x == (hi-1)) sawUB = true; 2402 assert(lo <= x && x < hi); 2403 } 2404 assert(sawLB && sawUB); 2405 } 2406 2407 { 2408 foreach (i; 0 .. 30) 2409 { 2410 assert(i == uniform(i, i+1, reproRng)); 2411 } 2412 } 2413 } 2414 2415 /+ 2416 Generates an unsigned integer in the half-open range `[0, k)`. 2417 Non-public because we locally guarantee `k > 0`. 2418 2419 Params: 2420 k = unsigned exclusive upper bound; caller guarantees this is non-zero 2421 rng = random number generator to use 2422 2423 Returns: 2424 Pseudo-random unsigned integer strictly less than `k`. 2425 +/ 2426 private UInt _uniformIndex(UniformRNG, UInt = size_t)(const UInt k, ref UniformRNG rng) 2427 if (isUnsigned!UInt && isUniformRNG!UniformRNG) 2428 { 2429 alias ResultType = UInt; 2430 alias UpperType = Unsigned!(typeof(k - 0)); 2431 alias upperDist = k; 2432 2433 assert(upperDist != 0); 2434 2435 // For backwards compatibility use same algorithm as uniform(0, k, rng). 2436 UpperType offset, rnum, bucketFront; 2437 do 2438 { 2439 rnum = uniform!UpperType(rng); 2440 offset = rnum % upperDist; 2441 bucketFront = rnum - offset; 2442 } // while we're in an unfair bucket... 2443 while (bucketFront > (UpperType.max - (upperDist - 1))); 2444 2445 return cast(ResultType) offset; 2446 } 2447 2448 pure @safe unittest 2449 { 2450 // For backwards compatibility check that _uniformIndex(k, rng) 2451 // has the same result as uniform(0, k, rng). 2452 auto rng1 = Xorshift(123_456_789); 2453 auto rng2 = rng1.save(); 2454 const size_t k = (1U << 31) - 1; 2455 assert(_uniformIndex(k, rng1) == uniform(0, k, rng2)); 2456 } 2457 2458 /** 2459 Generates a uniformly-distributed number in the range $(D [T.min, 2460 T.max]) for any integral or character type `T`. If no random 2461 number generator is passed, uses the default `rndGen`. 2462 2463 If an `enum` is used as type, the random variate is drawn with 2464 equal probability from any of the possible values of the enum `E`. 2465 2466 Params: 2467 urng = (optional) random number generator to use; 2468 if not specified, defaults to `rndGen` 2469 2470 Returns: 2471 Random variate drawn from the _uniform distribution across all 2472 possible values of the integral, character or enum type `T`. 2473 */ 2474 auto uniform(T, UniformRandomNumberGenerator) 2475 (ref UniformRandomNumberGenerator urng) 2476 if (!is(T == enum) && (isIntegral!T || isSomeChar!T) && isUniformRNG!UniformRandomNumberGenerator) 2477 { 2478 /* dchar does not use its full bit range, so we must 2479 * revert to the uniform with specified bounds 2480 */ 2481 static if (is(immutable T == immutable dchar)) 2482 { 2483 return uniform!"[]"(T.min, T.max, urng); 2484 } 2485 else 2486 { 2487 auto r = urng.front; 2488 urng.popFront(); 2489 static if (T.sizeof <= r.sizeof) 2490 { 2491 return cast(T) r; 2492 } 2493 else 2494 { 2495 static assert(T.sizeof == 8 && r.sizeof == 4); 2496 T r1 = urng.front | (cast(T) r << 32); 2497 urng.popFront(); 2498 return r1; 2499 } 2500 } 2501 } 2502 2503 /// Ditto 2504 auto uniform(T)() 2505 if (!is(T == enum) && (isIntegral!T || isSomeChar!T)) 2506 { 2507 return uniform!T(rndGen); 2508 } 2509 2510 /// 2511 @safe unittest 2512 { 2513 auto rnd = MinstdRand0(42); 2514 2515 assert(rnd.uniform!ubyte == 102); 2516 assert(rnd.uniform!ulong == 4838462006927449017); 2517 2518 enum Fruit { apple, mango, pear } 2519 version (D_LP64) // https://issues.dlang.org/show_bug.cgi?id=15147 2520 assert(rnd.uniform!Fruit == Fruit.mango); 2521 } 2522 2523 @safe unittest 2524 { 2525 // https://issues.dlang.org/show_bug.cgi?id=21383 2526 auto rng1 = Xorshift32(123456789); 2527 auto rng2 = rng1.save; 2528 assert(rng1.uniform!dchar == rng2.uniform!dchar); 2529 // https://issues.dlang.org/show_bug.cgi?id=21384 2530 assert(rng1.uniform!(const shared dchar) <= dchar.max); 2531 // https://issues.dlang.org/show_bug.cgi?id=8671 2532 double t8671 = 1.0 - uniform(0.0, 1.0); 2533 } 2534 2535 @safe unittest 2536 { 2537 static foreach (T; std.meta.AliasSeq!(char, wchar, dchar, byte, ubyte, short, ushort, 2538 int, uint, long, ulong)) 2539 {{ 2540 T init = uniform!T(); 2541 size_t i = 50; 2542 while (--i && uniform!T() == init) {} 2543 assert(i > 0); 2544 2545 foreach (immutable _; 0 .. 100) 2546 { 2547 auto u = uniform!T(); 2548 static assert(is(typeof(u) == T)); 2549 assert(T.min <= u, "Lower bound violation for uniform!" ~ T.stringof); 2550 assert(u <= T.max, "Upper bound violation for uniform!" ~ T.stringof); 2551 } 2552 }} 2553 } 2554 2555 /// ditto 2556 auto uniform(E, UniformRandomNumberGenerator) 2557 (ref UniformRandomNumberGenerator urng) 2558 if (is(E == enum) && isUniformRNG!UniformRandomNumberGenerator) 2559 { 2560 static immutable E[EnumMembers!E.length] members = [EnumMembers!E]; 2561 return members[std.random.uniform(0, members.length, urng)]; 2562 } 2563 2564 /// Ditto 2565 auto uniform(E)() 2566 if (is(E == enum)) 2567 { 2568 return uniform!E(rndGen); 2569 } 2570 2571 @safe unittest 2572 { 2573 enum Fruit { Apple = 12, Mango = 29, Pear = 72 } 2574 foreach (_; 0 .. 100) 2575 { 2576 foreach (f; [uniform!Fruit(), rndGen.uniform!Fruit()]) 2577 { 2578 assert(f == Fruit.Apple || f == Fruit.Mango || f == Fruit.Pear); 2579 } 2580 } 2581 } 2582 2583 /** 2584 * Generates a uniformly-distributed floating point number of type 2585 * `T` in the range [0, 1$(RPAREN). If no random number generator is 2586 * specified, the default RNG `rndGen` will be used as the source 2587 * of randomness. 2588 * 2589 * `uniform01` offers a faster generation of random variates than 2590 * the equivalent $(D uniform!"[$(RPAREN)"(0.0, 1.0)) and so may be preferred 2591 * for some applications. 2592 * 2593 * Params: 2594 * rng = (optional) random number generator to use; 2595 * if not specified, defaults to `rndGen` 2596 * 2597 * Returns: 2598 * Floating-point random variate of type `T` drawn from the _uniform 2599 * distribution across the half-open interval [0, 1$(RPAREN). 2600 * 2601 */ 2602 T uniform01(T = double)() 2603 if (isFloatingPoint!T) 2604 { 2605 return uniform01!T(rndGen); 2606 } 2607 2608 /// ditto 2609 T uniform01(T = double, UniformRNG)(ref UniformRNG rng) 2610 if (isFloatingPoint!T && isUniformRNG!UniformRNG) 2611 out (result) 2612 { 2613 assert(0 <= result); 2614 assert(result < 1); 2615 } 2616 do 2617 { 2618 alias R = typeof(rng.front); 2619 static if (isIntegral!R) 2620 { 2621 enum T factor = 1 / (T(1) + rng.max - rng.min); 2622 } 2623 else static if (isFloatingPoint!R) 2624 { 2625 enum T factor = 1 / (rng.max - rng.min); 2626 } 2627 else 2628 { 2629 static assert(false); 2630 } 2631 2632 while (true) 2633 { 2634 immutable T u = (rng.front - rng.min) * factor; 2635 rng.popFront(); 2636 2637 static if (isIntegral!R && T.mant_dig >= (8 * R.sizeof)) 2638 { 2639 /* If RNG variates are integral and T has enough precision to hold 2640 * R without loss, we're guaranteed by the definition of factor 2641 * that precisely u < 1. 2642 */ 2643 return u; 2644 } 2645 else 2646 { 2647 /* Otherwise we have to check whether u is beyond the assumed range 2648 * because of the loss of precision, or for another reason, a 2649 * floating-point RNG can return a variate that is exactly equal to 2650 * its maximum. 2651 */ 2652 if (u < 1) 2653 { 2654 return u; 2655 } 2656 } 2657 } 2658 2659 // Shouldn't ever get here. 2660 assert(false); 2661 } 2662 2663 /// 2664 @safe @nogc unittest 2665 { 2666 import std.math.operations : feqrel; 2667 2668 auto rnd = MinstdRand0(42); 2669 2670 // Generate random numbers in the range in the range [0, 1) 2671 auto u1 = uniform01(rnd); 2672 assert(u1 >= 0 && u1 < 1); 2673 2674 auto u2 = rnd.uniform01!float; 2675 assert(u2 >= 0 && u2 < 1); 2676 2677 // Confirm that the random values with the initial seed 42 are 0.000328707 and 0.524587 2678 assert(u1.feqrel(0.000328707) > 20); 2679 assert(u2.feqrel(0.524587) > 20); 2680 } 2681 2682 @safe @nogc unittest 2683 { 2684 import std.meta; 2685 static foreach (UniformRNG; PseudoRngTypes) 2686 {{ 2687 2688 static foreach (T; std.meta.AliasSeq!(float, double, real)) 2689 {{ 2690 UniformRNG rng = UniformRNG(123_456_789); 2691 2692 auto a = uniform01(); 2693 assert(is(typeof(a) == double)); 2694 assert(0 <= a && a < 1); 2695 2696 auto b = uniform01(rng); 2697 assert(is(typeof(a) == double)); 2698 assert(0 <= b && b < 1); 2699 2700 auto c = uniform01!T(); 2701 assert(is(typeof(c) == T)); 2702 assert(0 <= c && c < 1); 2703 2704 auto d = uniform01!T(rng); 2705 assert(is(typeof(d) == T)); 2706 assert(0 <= d && d < 1); 2707 2708 T init = uniform01!T(rng); 2709 size_t i = 50; 2710 while (--i && uniform01!T(rng) == init) {} 2711 assert(i > 0); 2712 assert(i < 50); 2713 }} 2714 }} 2715 } 2716 2717 /** 2718 Generates a uniform probability distribution of size `n`, i.e., an 2719 array of size `n` of positive numbers of type `F` that sum to 2720 `1`. If `useThis` is provided, it is used as storage. 2721 */ 2722 F[] uniformDistribution(F = double)(size_t n, F[] useThis = null) 2723 if (isFloatingPoint!F) 2724 { 2725 import std.numeric : normalize; 2726 useThis.length = n; 2727 foreach (ref e; useThis) 2728 { 2729 e = uniform(0.0, 1); 2730 } 2731 normalize(useThis); 2732 return useThis; 2733 } 2734 2735 /// 2736 @safe unittest 2737 { 2738 import std.algorithm.iteration : reduce; 2739 import std.math.operations : isClose; 2740 2741 auto a = uniformDistribution(5); 2742 assert(a.length == 5); 2743 assert(isClose(reduce!"a + b"(a), 1)); 2744 2745 a = uniformDistribution(10, a); 2746 assert(a.length == 10); 2747 assert(isClose(reduce!"a + b"(a), 1)); 2748 } 2749 2750 /** 2751 Returns a random, uniformly chosen, element `e` from the supplied 2752 $(D Range range). If no random number generator is passed, the default 2753 `rndGen` is used. 2754 2755 Params: 2756 range = a random access range that has the `length` property defined 2757 urng = (optional) random number generator to use; 2758 if not specified, defaults to `rndGen` 2759 2760 Returns: 2761 A single random element drawn from the `range`. If it can, it will 2762 return a `ref` to the $(D range element), otherwise it will return 2763 a copy. 2764 */ 2765 auto ref choice(Range, RandomGen = Random)(Range range, ref RandomGen urng) 2766 if (isRandomAccessRange!Range && hasLength!Range && isUniformRNG!RandomGen) 2767 { 2768 assert(range.length > 0, 2769 __PRETTY_FUNCTION__ ~ ": invalid Range supplied. Range cannot be empty"); 2770 2771 return range[uniform(size_t(0), $, urng)]; 2772 } 2773 2774 /// ditto 2775 auto ref choice(Range)(Range range) 2776 { 2777 return choice(range, rndGen); 2778 } 2779 2780 /// ditto 2781 auto ref choice(Range, RandomGen = Random)(ref Range range, ref RandomGen urng) 2782 if (isRandomAccessRange!Range && hasLength!Range && isUniformRNG!RandomGen) 2783 { 2784 assert(range.length > 0, 2785 __PRETTY_FUNCTION__ ~ ": invalid Range supplied. Range cannot be empty"); 2786 return range[uniform(size_t(0), $, urng)]; 2787 } 2788 2789 /// ditto 2790 auto ref choice(Range)(ref Range range) 2791 { 2792 return choice(range, rndGen); 2793 } 2794 2795 /// 2796 @safe unittest 2797 { 2798 auto rnd = MinstdRand0(42); 2799 2800 auto elem = [1, 2, 3, 4, 5].choice(rnd); 2801 version (D_LP64) // https://issues.dlang.org/show_bug.cgi?id=15147 2802 assert(elem == 3); 2803 } 2804 2805 @safe unittest 2806 { 2807 import std.algorithm.searching : canFind; 2808 2809 static class MyTestClass 2810 { 2811 int x; 2812 2813 this(int x) 2814 { 2815 this.x = x; 2816 } 2817 } 2818 2819 MyTestClass[] testClass; 2820 foreach (i; 0 .. 5) 2821 { 2822 testClass ~= new MyTestClass(i); 2823 } 2824 2825 auto elem = choice(testClass); 2826 2827 assert(canFind!((ref MyTestClass a, ref MyTestClass b) => a.x == b.x)(testClass, elem), 2828 "Choice did not return a valid element from the given Range"); 2829 } 2830 2831 @system unittest 2832 { 2833 import std.algorithm.iteration : map; 2834 import std.algorithm.searching : canFind; 2835 2836 auto array = [1, 2, 3, 4, 5]; 2837 auto elemAddr = &choice(array); 2838 2839 assert(array.map!((ref e) => &e).canFind(elemAddr), 2840 "Choice did not return a ref to an element from the given Range"); 2841 assert(array.canFind(*(cast(int *)(elemAddr))), 2842 "Choice did not return a valid element from the given Range"); 2843 } 2844 2845 @safe unittest // https://issues.dlang.org/show_bug.cgi?id=18631 2846 { 2847 auto rng = MinstdRand0(42); 2848 const a = [0,1,2]; 2849 const(int[]) b = [0, 1, 2]; 2850 auto x = choice(a); 2851 auto y = choice(b); 2852 auto z = choice(cast(const)[1, 2, 3]); 2853 auto x1 = choice(a, rng); 2854 auto y1 = choice(b, rng); 2855 auto z1 = choice(cast(const)[1, 2, 3], rng); 2856 } 2857 2858 @safe unittest // Ref range (https://issues.dlang.org/show_bug.cgi?id=18631 PR) 2859 { 2860 struct TestRange 2861 { 2862 int x; 2863 ref int front() return {return x;} 2864 ref int back() return {return x;} 2865 void popFront() {} 2866 void popBack() {} 2867 bool empty = false; 2868 TestRange save() {return this;} 2869 size_t length = 10; 2870 alias opDollar = length; 2871 ref int opIndex(size_t i) return {return x;} 2872 } 2873 2874 TestRange r = TestRange(10); 2875 int* s = &choice(r); 2876 } 2877 2878 /** 2879 Shuffles elements of `r` using `gen` as a shuffler. `r` must be 2880 a random-access range with length. If no RNG is specified, `rndGen` 2881 will be used. 2882 2883 Params: 2884 r = random-access range whose elements are to be shuffled 2885 gen = (optional) random number generator to use; if not 2886 specified, defaults to `rndGen` 2887 Returns: 2888 The shuffled random-access range. 2889 */ 2890 2891 Range randomShuffle(Range, RandomGen)(Range r, ref RandomGen gen) 2892 if (isRandomAccessRange!Range && isUniformRNG!RandomGen) 2893 { 2894 import std.algorithm.mutation : swapAt; 2895 const n = r.length; 2896 foreach (i; 0 .. n) 2897 { 2898 r.swapAt(i, i + _uniformIndex(n - i, gen)); 2899 } 2900 return r; 2901 } 2902 2903 /// ditto 2904 Range randomShuffle(Range)(Range r) 2905 if (isRandomAccessRange!Range) 2906 { 2907 return randomShuffle(r, rndGen); 2908 } 2909 2910 /// 2911 @safe unittest 2912 { 2913 auto rnd = MinstdRand0(42); 2914 2915 auto arr = [1, 2, 3, 4, 5].randomShuffle(rnd); 2916 version (D_LP64) // https://issues.dlang.org/show_bug.cgi?id=15147 2917 assert(arr == [3, 5, 2, 4, 1]); 2918 } 2919 2920 @safe unittest 2921 { 2922 int[10] sa = void; 2923 int[10] sb = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]; 2924 import std.algorithm.sorting : sort; 2925 foreach (RandomGen; PseudoRngTypes) 2926 { 2927 sa[] = sb[]; 2928 auto a = sa[]; 2929 auto b = sb[]; 2930 auto gen = RandomGen(123_456_789); 2931 randomShuffle(a, gen); 2932 sort(a); 2933 assert(a == b); 2934 randomShuffle(a); 2935 sort(a); 2936 assert(a == b); 2937 } 2938 // For backwards compatibility verify randomShuffle(r, gen) 2939 // is equivalent to partialShuffle(r, 0, r.length, gen). 2940 auto gen1 = Xorshift(123_456_789); 2941 auto gen2 = gen1.save(); 2942 sa[] = sb[]; 2943 // @nogc std.random.randomShuffle. 2944 // https://issues.dlang.org/show_bug.cgi?id=19156 2945 () @nogc nothrow pure { randomShuffle(sa[], gen1); }(); 2946 partialShuffle(sb[], sb.length, gen2); 2947 assert(sa[] == sb[]); 2948 } 2949 2950 // https://issues.dlang.org/show_bug.cgi?id=18501 2951 @safe unittest 2952 { 2953 import std.algorithm.comparison : among; 2954 auto r = randomShuffle([0,1]); 2955 assert(r.among([0,1],[1,0])); 2956 } 2957 2958 /** 2959 Partially shuffles the elements of `r` such that upon returning $(D r[0 .. n]) 2960 is a random subset of `r` and is randomly ordered. $(D r[n .. r.length]) 2961 will contain the elements not in $(D r[0 .. n]). These will be in an undefined 2962 order, but will not be random in the sense that their order after 2963 `partialShuffle` returns will not be independent of their order before 2964 `partialShuffle` was called. 2965 2966 `r` must be a random-access range with length. `n` must be less than 2967 or equal to `r.length`. If no RNG is specified, `rndGen` will be used. 2968 2969 Params: 2970 r = random-access range whose elements are to be shuffled 2971 n = number of elements of `r` to shuffle (counting from the beginning); 2972 must be less than `r.length` 2973 gen = (optional) random number generator to use; if not 2974 specified, defaults to `rndGen` 2975 Returns: 2976 The shuffled random-access range. 2977 */ 2978 Range partialShuffle(Range, RandomGen)(Range r, in size_t n, ref RandomGen gen) 2979 if (isRandomAccessRange!Range && isUniformRNG!RandomGen) 2980 { 2981 import std.algorithm.mutation : swapAt; 2982 import std.exception : enforce; 2983 enforce(n <= r.length, "n must be <= r.length for partialShuffle."); 2984 foreach (i; 0 .. n) 2985 { 2986 r.swapAt(i, uniform(i, r.length, gen)); 2987 } 2988 return r; 2989 } 2990 2991 /// ditto 2992 Range partialShuffle(Range)(Range r, in size_t n) 2993 if (isRandomAccessRange!Range) 2994 { 2995 return partialShuffle(r, n, rndGen); 2996 } 2997 2998 /// 2999 @safe unittest 3000 { 3001 auto rnd = MinstdRand0(42); 3002 3003 auto arr = [1, 2, 3, 4, 5, 6]; 3004 arr = arr.dup.partialShuffle(1, rnd); 3005 3006 version (D_LP64) // https://issues.dlang.org/show_bug.cgi?id=15147 3007 assert(arr == [2, 1, 3, 4, 5, 6]); // 1<->2 3008 3009 arr = arr.dup.partialShuffle(2, rnd); 3010 version (D_LP64) // https://issues.dlang.org/show_bug.cgi?id=15147 3011 assert(arr == [1, 4, 3, 2, 5, 6]); // 1<->2, 2<->4 3012 3013 arr = arr.dup.partialShuffle(3, rnd); 3014 version (D_LP64) // https://issues.dlang.org/show_bug.cgi?id=15147 3015 assert(arr == [5, 4, 6, 2, 1, 3]); // 1<->5, 2<->4, 3<->6 3016 } 3017 3018 @safe unittest 3019 { 3020 import std.algorithm; 3021 foreach (RandomGen; PseudoRngTypes) 3022 { 3023 auto a = [0, 1, 1, 2, 3]; 3024 auto b = a.dup; 3025 3026 // Pick a fixed seed so that the outcome of the statistical 3027 // test below is deterministic. 3028 auto gen = RandomGen(12345); 3029 3030 // NUM times, pick LEN elements from the array at random. 3031 immutable int LEN = 2; 3032 immutable int NUM = 750; 3033 int[][] chk; 3034 foreach (step; 0 .. NUM) 3035 { 3036 partialShuffle(a, LEN, gen); 3037 chk ~= a[0 .. LEN].dup; 3038 } 3039 3040 // Check that each possible a[0 .. LEN] was produced at least once. 3041 // For a perfectly random RandomGen, the probability that each 3042 // particular combination failed to appear would be at most 3043 // 0.95 ^^ NUM which is approximately 1,962e-17. 3044 // As long as hardware failure (e.g. bit flip) probability 3045 // is higher, we are fine with this unittest. 3046 sort(chk); 3047 assert(equal(uniq(chk), [ [0,1], [0,2], [0,3], 3048 [1,0], [1,1], [1,2], [1,3], 3049 [2,0], [2,1], [2,3], 3050 [3,0], [3,1], [3,2], ])); 3051 3052 // Check that all the elements are still there. 3053 sort(a); 3054 assert(equal(a, b)); 3055 } 3056 } 3057 3058 /** 3059 Get a random index into a list of weights corresponding to each index 3060 3061 Similar to rolling a die with relative probabilities stored in `proportions`. 3062 Returns the index in `proportions` that was chosen. 3063 3064 Note: 3065 Usually, dice are 'fair', meaning that each side has equal probability 3066 to come up, in which case `1 + uniform(0, 6)` can simply be used. 3067 In future Phobos versions, this function might get renamed to something like 3068 `weightedChoice` to avoid confusion. 3069 3070 Params: 3071 rnd = (optional) random number generator to use; if not 3072 specified, defaults to `rndGen` 3073 proportions = forward range or list of individual values 3074 whose elements correspond to the probabilities 3075 with which to choose the corresponding index 3076 value 3077 3078 Returns: 3079 Random variate drawn from the index values 3080 [0, ... `proportions.length` - 1], with the probability 3081 of getting an individual index value `i` being proportional to 3082 `proportions[i]`. 3083 */ 3084 size_t dice(Rng, Num)(ref Rng rnd, Num[] proportions...) 3085 if (isNumeric!Num && isForwardRange!Rng) 3086 { 3087 return diceImpl(rnd, proportions); 3088 } 3089 3090 /// Ditto 3091 size_t dice(R, Range)(ref R rnd, Range proportions) 3092 if (isForwardRange!Range && isNumeric!(ElementType!Range) && !isArray!Range) 3093 { 3094 return diceImpl(rnd, proportions); 3095 } 3096 3097 /// Ditto 3098 size_t dice(Range)(Range proportions) 3099 if (isForwardRange!Range && isNumeric!(ElementType!Range) && !isArray!Range) 3100 { 3101 return diceImpl(rndGen, proportions); 3102 } 3103 3104 /// Ditto 3105 size_t dice(Num)(Num[] proportions...) 3106 if (isNumeric!Num) 3107 { 3108 return diceImpl(rndGen, proportions); 3109 } 3110 3111 /// 3112 @safe unittest 3113 { 3114 auto d6 = 1 + dice(1, 1, 1, 1, 1, 1); // fair dice roll 3115 auto d6b = 1 + dice(2, 1, 1, 1, 1, 1); // double the chance to roll '1' 3116 3117 auto x = dice(0.5, 0.5); // x is 0 or 1 in equal proportions 3118 auto y = dice(50, 50); // y is 0 or 1 in equal proportions 3119 auto z = dice(70, 20, 10); // z is 0 70% of the time, 1 20% of the time, 3120 // and 2 10% of the time 3121 } 3122 3123 /// 3124 @safe unittest 3125 { 3126 auto rnd = MinstdRand0(42); 3127 auto z = rnd.dice(70, 20, 10); 3128 assert(z == 0); 3129 z = rnd.dice(30, 20, 40, 10); 3130 assert(z == 2); 3131 } 3132 3133 private size_t diceImpl(Rng, Range)(ref Rng rng, scope Range proportions) 3134 if (isForwardRange!Range && isNumeric!(ElementType!Range) && isForwardRange!Rng) 3135 in 3136 { 3137 import std.algorithm.searching : all; 3138 assert(proportions.save.all!"a >= 0"); 3139 } 3140 do 3141 { 3142 import std.algorithm.iteration : reduce; 3143 import std.exception : enforce; 3144 double sum = reduce!"a + b"(0.0, proportions.save); 3145 enforce(sum > 0, "Proportions in a dice cannot sum to zero"); 3146 immutable point = uniform(0.0, sum, rng); 3147 assert(point < sum); 3148 auto mass = 0.0; 3149 3150 size_t i = 0; 3151 foreach (e; proportions) 3152 { 3153 mass += e; 3154 if (point < mass) return i; 3155 i++; 3156 } 3157 // this point should not be reached 3158 assert(false); 3159 } 3160 3161 /// 3162 @safe unittest 3163 { 3164 auto rnd = Xorshift(123_456_789); 3165 auto i = dice(rnd, 0.0, 100.0); 3166 assert(i == 1); 3167 i = dice(rnd, 100.0, 0.0); 3168 assert(i == 0); 3169 3170 i = dice(100U, 0U); 3171 assert(i == 0); 3172 } 3173 3174 /+ @nogc bool array designed for RandomCover. 3175 - constructed with an invariable length 3176 - small length means 0 alloc and bit field (if up to 32(x86) or 64(x64) choices to cover) 3177 - bigger length means non-GC heap allocation(s) and dealloc. +/ 3178 private struct RandomCoverChoices 3179 { 3180 private size_t* buffer; 3181 private immutable size_t _length; 3182 private immutable bool hasPackedBits; 3183 private enum BITS_PER_WORD = typeof(buffer[0]).sizeof * 8; 3184 3185 void opAssign(T)(T) @disable; 3186 3187 this(this) pure nothrow @nogc @trusted 3188 { 3189 import core.stdc.string : memcpy; 3190 import std.internal.memory : enforceMalloc; 3191 3192 if (!hasPackedBits && buffer !is null) 3193 { 3194 const nBytesToAlloc = size_t.sizeof * (_length / BITS_PER_WORD + int(_length % BITS_PER_WORD != 0)); 3195 void* nbuffer = enforceMalloc(nBytesToAlloc); 3196 buffer = cast(size_t*) memcpy(nbuffer, buffer, nBytesToAlloc); 3197 } 3198 } 3199 3200 this(size_t numChoices) pure nothrow @nogc @trusted 3201 { 3202 import std.internal.memory : enforceCalloc; 3203 3204 _length = numChoices; 3205 hasPackedBits = _length <= size_t.sizeof * 8; 3206 if (!hasPackedBits) 3207 { 3208 const nWordsToAlloc = _length / BITS_PER_WORD + int(_length % BITS_PER_WORD != 0); 3209 buffer = cast(size_t*) enforceCalloc(nWordsToAlloc, BITS_PER_WORD / 8); 3210 } 3211 } 3212 3213 size_t length() const pure nothrow @nogc @safe @property {return _length;} 3214 3215 ~this() pure nothrow @nogc @trusted 3216 { 3217 import core.memory : pureFree; 3218 3219 if (!hasPackedBits && buffer !is null) 3220 pureFree(buffer); 3221 } 3222 3223 bool opIndex(size_t index) const pure nothrow @nogc @trusted 3224 { 3225 assert(index < _length); 3226 import core.bitop : bt; 3227 if (!hasPackedBits) 3228 return cast(bool) bt(buffer, index); 3229 else 3230 return ((cast(size_t) buffer) >> index) & size_t(1); 3231 } 3232 3233 void opIndexAssign(bool value, size_t index) pure nothrow @nogc @trusted 3234 { 3235 assert(index < _length); 3236 if (!hasPackedBits) 3237 { 3238 import core.bitop : btr, bts; 3239 if (value) 3240 bts(buffer, index); 3241 else 3242 btr(buffer, index); 3243 } 3244 else 3245 { 3246 if (value) 3247 (*cast(size_t*) &buffer) |= size_t(1) << index; 3248 else 3249 (*cast(size_t*) &buffer) &= ~(size_t(1) << index); 3250 } 3251 } 3252 } 3253 3254 @safe @nogc nothrow unittest 3255 { 3256 static immutable lengths = [3, 32, 65, 256]; 3257 foreach (length; lengths) 3258 { 3259 RandomCoverChoices c = RandomCoverChoices(length); 3260 assert(c.hasPackedBits == (length <= size_t.sizeof * 8)); 3261 c[0] = true; 3262 c[2] = true; 3263 assert(c[0]); 3264 assert(!c[1]); 3265 assert(c[2]); 3266 c[0] = false; 3267 c[1] = true; 3268 c[2] = false; 3269 assert(!c[0]); 3270 assert(c[1]); 3271 assert(!c[2]); 3272 } 3273 } 3274 3275 /** 3276 Covers a given range `r` in a random manner, i.e. goes through each 3277 element of `r` once and only once, just in a random order. `r` 3278 must be a random-access range with length. 3279 3280 If no random number generator is passed to `randomCover`, the 3281 thread-global RNG rndGen will be used internally. 3282 3283 Params: 3284 r = random-access range to cover 3285 rng = (optional) random number generator to use; 3286 if not specified, defaults to `rndGen` 3287 3288 Returns: 3289 Range whose elements consist of the elements of `r`, 3290 in random order. Will be a forward range if both `r` and 3291 `rng` are forward ranges, an 3292 $(REF_ALTTEXT input range, isInputRange, std,range,primitives) otherwise. 3293 */ 3294 struct RandomCover(Range, UniformRNG = void) 3295 if (isRandomAccessRange!Range && (isUniformRNG!UniformRNG || is(UniformRNG == void))) 3296 { 3297 private Range _input; 3298 private RandomCoverChoices _chosen; 3299 private size_t _current; 3300 private size_t _alreadyChosen = 0; 3301 private bool _isEmpty = false; 3302 3303 static if (is(UniformRNG == void)) 3304 { 3305 this(Range input) 3306 { 3307 _input = input; 3308 _chosen = RandomCoverChoices(_input.length); 3309 if (_input.empty) 3310 { 3311 _isEmpty = true; 3312 } 3313 else 3314 { 3315 _current = _uniformIndex(_chosen.length, rndGen); 3316 } 3317 } 3318 } 3319 else 3320 { 3321 private UniformRNG _rng; 3322 3323 this(Range input, ref UniformRNG rng) 3324 { 3325 _input = input; 3326 _rng = rng; 3327 _chosen = RandomCoverChoices(_input.length); 3328 if (_input.empty) 3329 { 3330 _isEmpty = true; 3331 } 3332 else 3333 { 3334 _current = _uniformIndex(_chosen.length, rng); 3335 } 3336 } 3337 3338 this(Range input, UniformRNG rng) 3339 { 3340 this(input, rng); 3341 } 3342 } 3343 3344 static if (hasLength!Range) 3345 { 3346 @property size_t length() 3347 { 3348 return _input.length - _alreadyChosen; 3349 } 3350 } 3351 3352 @property auto ref front() 3353 { 3354 assert(!_isEmpty); 3355 return _input[_current]; 3356 } 3357 3358 void popFront() 3359 { 3360 assert(!_isEmpty); 3361 3362 size_t k = _input.length - _alreadyChosen - 1; 3363 if (k == 0) 3364 { 3365 _isEmpty = true; 3366 ++_alreadyChosen; 3367 return; 3368 } 3369 3370 size_t i; 3371 foreach (e; _input) 3372 { 3373 if (_chosen[i] || i == _current) { ++i; continue; } 3374 // Roll a dice with k faces 3375 static if (is(UniformRNG == void)) 3376 { 3377 auto chooseMe = _uniformIndex(k, rndGen) == 0; 3378 } 3379 else 3380 { 3381 auto chooseMe = _uniformIndex(k, _rng) == 0; 3382 } 3383 assert(k > 1 || chooseMe); 3384 if (chooseMe) 3385 { 3386 _chosen[_current] = true; 3387 _current = i; 3388 ++_alreadyChosen; 3389 return; 3390 } 3391 --k; 3392 ++i; 3393 } 3394 } 3395 3396 static if (isForwardRange!UniformRNG) 3397 { 3398 @property typeof(this) save() 3399 { 3400 auto ret = this; 3401 ret._input = _input.save; 3402 ret._rng = _rng.save; 3403 return ret; 3404 } 3405 } 3406 3407 @property bool empty() const { return _isEmpty; } 3408 } 3409 3410 /// Ditto 3411 auto randomCover(Range, UniformRNG)(Range r, auto ref UniformRNG rng) 3412 if (isRandomAccessRange!Range && isUniformRNG!UniformRNG) 3413 { 3414 return RandomCover!(Range, UniformRNG)(r, rng); 3415 } 3416 3417 /// Ditto 3418 auto randomCover(Range)(Range r) 3419 if (isRandomAccessRange!Range) 3420 { 3421 return RandomCover!(Range, void)(r); 3422 } 3423 3424 /// 3425 @safe unittest 3426 { 3427 import std.algorithm.comparison : equal; 3428 import std.range : iota; 3429 auto rnd = MinstdRand0(42); 3430 3431 version (D_LP64) // https://issues.dlang.org/show_bug.cgi?id=15147 3432 assert(10.iota.randomCover(rnd).equal([7, 4, 2, 0, 1, 6, 8, 3, 9, 5])); 3433 } 3434 3435 @safe unittest // cover RandomCoverChoices postblit for heap storage 3436 { 3437 import std.array : array; 3438 import std.range : iota; 3439 auto a = 1337.iota.randomCover().array; 3440 assert(a.length == 1337); 3441 } 3442 3443 @nogc nothrow pure @safe unittest 3444 { 3445 // Optionally @nogc std.random.randomCover 3446 // https://issues.dlang.org/show_bug.cgi?id=14001 3447 auto rng = Xorshift(123_456_789); 3448 static immutable int[] sa = [1, 2, 3, 4, 5]; 3449 auto r = randomCover(sa, rng); 3450 assert(!r.empty); 3451 const x = r.front; 3452 r.popFront(); 3453 assert(!r.empty); 3454 const y = r.front; 3455 assert(x != y); 3456 } 3457 3458 @safe unittest 3459 { 3460 import std.algorithm; 3461 import std.conv; 3462 int[] a = [ 0, 1, 2, 3, 4, 5, 6, 7, 8 ]; 3463 int[] c; 3464 static foreach (UniformRNG; std.meta.AliasSeq!(void, PseudoRngTypes)) 3465 {{ 3466 static if (is(UniformRNG == void)) 3467 { 3468 auto rc = randomCover(a); 3469 static assert(isInputRange!(typeof(rc))); 3470 static assert(!isForwardRange!(typeof(rc))); 3471 } 3472 else 3473 { 3474 auto rng = UniformRNG(123_456_789); 3475 auto rc = randomCover(a, rng); 3476 static assert(isForwardRange!(typeof(rc))); 3477 // check for constructor passed a value-type RNG 3478 auto rc2 = RandomCover!(int[], UniformRNG)(a, UniformRNG(987_654_321)); 3479 static assert(isForwardRange!(typeof(rc2))); 3480 auto rcEmpty = randomCover(c, rng); 3481 assert(rcEmpty.length == 0); 3482 } 3483 3484 int[] b = new int[9]; 3485 uint i; 3486 foreach (e; rc) 3487 { 3488 //writeln(e); 3489 b[i++] = e; 3490 } 3491 sort(b); 3492 assert(a == b, text(b)); 3493 }} 3494 } 3495 3496 @safe unittest 3497 { 3498 // https://issues.dlang.org/show_bug.cgi?id=12589 3499 int[] r = []; 3500 auto rc = randomCover(r); 3501 assert(rc.length == 0); 3502 assert(rc.empty); 3503 3504 // https://issues.dlang.org/show_bug.cgi?id=16724 3505 import std.range : iota; 3506 auto range = iota(10); 3507 auto randy = range.randomCover; 3508 3509 for (int i=1; i <= range.length; i++) 3510 { 3511 randy.popFront; 3512 assert(randy.length == range.length - i); 3513 } 3514 } 3515 3516 // RandomSample 3517 /** 3518 Selects a random subsample out of `r`, containing exactly `n` 3519 elements. The order of elements is the same as in the original 3520 range. The total length of `r` must be known. If `total` is 3521 passed in, the total number of sample is considered to be $(D 3522 total). Otherwise, `RandomSample` uses `r.length`. 3523 3524 Params: 3525 r = range to sample from 3526 n = number of elements to include in the sample; 3527 must be less than or equal to the total number 3528 of elements in `r` and/or the parameter 3529 `total` (if provided) 3530 total = (semi-optional) number of elements of `r` 3531 from which to select the sample (counting from 3532 the beginning); must be less than or equal to 3533 the total number of elements in `r` itself. 3534 May be omitted if `r` has the `.length` 3535 property and the sample is to be drawn from 3536 all elements of `r`. 3537 rng = (optional) random number generator to use; 3538 if not specified, defaults to `rndGen` 3539 3540 Returns: 3541 Range whose elements consist of a randomly selected subset of 3542 the elements of `r`, in the same order as these elements 3543 appear in `r` itself. Will be a forward range if both `r` 3544 and `rng` are forward ranges, an input range otherwise. 3545 3546 `RandomSample` implements Jeffrey Scott Vitter's Algorithm D 3547 (see Vitter $(HTTP dx.doi.org/10.1145/358105.893, 1984), $(HTTP 3548 dx.doi.org/10.1145/23002.23003, 1987)), which selects a sample 3549 of size `n` in O(n) steps and requiring O(n) random variates, 3550 regardless of the size of the data being sampled. The exception 3551 to this is if traversing k elements on the input range is itself 3552 an O(k) operation (e.g. when sampling lines from an input file), 3553 in which case the sampling calculation will inevitably be of 3554 O(total). 3555 3556 RandomSample will throw an exception if `total` is verifiably 3557 less than the total number of elements available in the input, 3558 or if $(D n > total). 3559 3560 If no random number generator is passed to `randomSample`, the 3561 thread-global RNG rndGen will be used internally. 3562 */ 3563 struct RandomSample(Range, UniformRNG = void) 3564 if (isInputRange!Range && (isUniformRNG!UniformRNG || is(UniformRNG == void))) 3565 { 3566 private size_t _available, _toSelect; 3567 private enum ushort _alphaInverse = 13; // Vitter's recommended value. 3568 private double _Vprime; 3569 private Range _input; 3570 private size_t _index; 3571 private enum Skip { None, A, D } 3572 private Skip _skip = Skip.None; 3573 3574 // If we're using the default thread-local random number generator then 3575 // we shouldn't store a copy of it here. UniformRNG == void is a sentinel 3576 // for this. If we're using a user-specified generator then we have no 3577 // choice but to store a copy. 3578 static if (is(UniformRNG == void)) 3579 { 3580 static if (hasLength!Range) 3581 { 3582 this(Range input, size_t howMany) 3583 { 3584 _input = input; 3585 initialize(howMany, input.length); 3586 } 3587 } 3588 3589 this(Range input, size_t howMany, size_t total) 3590 { 3591 _input = input; 3592 initialize(howMany, total); 3593 } 3594 } 3595 else 3596 { 3597 UniformRNG _rng; 3598 3599 static if (hasLength!Range) 3600 { 3601 this(Range input, size_t howMany, ref scope UniformRNG rng) 3602 { 3603 _rng = rng; 3604 _input = input; 3605 initialize(howMany, input.length); 3606 } 3607 3608 this(Range input, size_t howMany, UniformRNG rng) 3609 { 3610 this(input, howMany, rng); 3611 } 3612 } 3613 3614 this(Range input, size_t howMany, size_t total, ref scope UniformRNG rng) 3615 { 3616 _rng = rng; 3617 _input = input; 3618 initialize(howMany, total); 3619 } 3620 3621 this(Range input, size_t howMany, size_t total, UniformRNG rng) 3622 { 3623 this(input, howMany, total, rng); 3624 } 3625 } 3626 3627 private void initialize(size_t howMany, size_t total) 3628 { 3629 import std.conv : text; 3630 import std.exception : enforce; 3631 _available = total; 3632 _toSelect = howMany; 3633 enforce(_toSelect <= _available, 3634 text("RandomSample: cannot sample ", _toSelect, 3635 " items when only ", _available, " are available")); 3636 static if (hasLength!Range) 3637 { 3638 enforce(_available <= _input.length, 3639 text("RandomSample: specified ", _available, 3640 " items as available when input contains only ", 3641 _input.length)); 3642 } 3643 } 3644 3645 private void initializeFront() 3646 { 3647 assert(_skip == Skip.None); 3648 // We can save ourselves a random variate by checking right 3649 // at the beginning if we should use Algorithm A. 3650 if ((_alphaInverse * _toSelect) > _available) 3651 { 3652 _skip = Skip.A; 3653 } 3654 else 3655 { 3656 _skip = Skip.D; 3657 _Vprime = newVprime(_toSelect); 3658 } 3659 prime(); 3660 } 3661 3662 /** 3663 Range primitives. 3664 */ 3665 @property bool empty() const 3666 { 3667 return _toSelect == 0; 3668 } 3669 3670 /// Ditto 3671 @property auto ref front() 3672 { 3673 assert(!empty); 3674 // The first sample point must be determined here to avoid 3675 // having it always correspond to the first element of the 3676 // input. The rest of the sample points are determined each 3677 // time we call popFront(). 3678 if (_skip == Skip.None) 3679 { 3680 initializeFront(); 3681 } 3682 return _input.front; 3683 } 3684 3685 /// Ditto 3686 void popFront() 3687 { 3688 // First we need to check if the sample has 3689 // been initialized in the first place. 3690 if (_skip == Skip.None) 3691 { 3692 initializeFront(); 3693 } 3694 3695 _input.popFront(); 3696 --_available; 3697 --_toSelect; 3698 ++_index; 3699 prime(); 3700 } 3701 3702 /// Ditto 3703 static if (isForwardRange!Range && isForwardRange!UniformRNG) 3704 { 3705 static if (is(typeof(((const UniformRNG* p) => (*p).save)(null)) : UniformRNG) 3706 && is(typeof(((const Range* p) => (*p).save)(null)) : Range)) 3707 { 3708 @property typeof(this) save() const 3709 { 3710 auto ret = RandomSample.init; 3711 foreach (fieldIndex, ref val; this.tupleof) 3712 { 3713 static if (is(typeof(val) == const(Range)) || is(typeof(val) == const(UniformRNG))) 3714 ret.tupleof[fieldIndex] = val.save; 3715 else 3716 ret.tupleof[fieldIndex] = val; 3717 } 3718 return ret; 3719 } 3720 } 3721 else 3722 { 3723 @property typeof(this) save() 3724 { 3725 auto ret = this; 3726 ret._input = _input.save; 3727 ret._rng = _rng.save; 3728 return ret; 3729 } 3730 } 3731 } 3732 3733 /// Ditto 3734 @property size_t length() const 3735 { 3736 return _toSelect; 3737 } 3738 3739 /** 3740 Returns the index of the visited record. 3741 */ 3742 @property size_t index() 3743 { 3744 if (_skip == Skip.None) 3745 { 3746 initializeFront(); 3747 } 3748 return _index; 3749 } 3750 3751 private size_t skip() 3752 { 3753 assert(_skip != Skip.None); 3754 3755 // Step D1: if the number of points still to select is greater 3756 // than a certain proportion of the remaining data points, i.e. 3757 // if n >= alpha * N where alpha = 1/13, we carry out the 3758 // sampling with Algorithm A. 3759 if (_skip == Skip.A) 3760 { 3761 return skipA(); 3762 } 3763 else if ((_alphaInverse * _toSelect) > _available) 3764 { 3765 // We shouldn't get here unless the current selected 3766 // algorithm is D. 3767 assert(_skip == Skip.D); 3768 _skip = Skip.A; 3769 return skipA(); 3770 } 3771 else 3772 { 3773 assert(_skip == Skip.D); 3774 return skipD(); 3775 } 3776 } 3777 3778 /* 3779 Vitter's Algorithm A, used when the ratio of needed sample values 3780 to remaining data values is sufficiently large. 3781 */ 3782 private size_t skipA() 3783 { 3784 size_t s; 3785 double v, quot, top; 3786 3787 if (_toSelect == 1) 3788 { 3789 static if (is(UniformRNG == void)) 3790 { 3791 s = uniform(0, _available); 3792 } 3793 else 3794 { 3795 s = uniform(0, _available, _rng); 3796 } 3797 } 3798 else 3799 { 3800 v = 0; 3801 top = _available - _toSelect; 3802 quot = top / _available; 3803 3804 static if (is(UniformRNG == void)) 3805 { 3806 v = uniform!"()"(0.0, 1.0); 3807 } 3808 else 3809 { 3810 v = uniform!"()"(0.0, 1.0, _rng); 3811 } 3812 3813 while (quot > v) 3814 { 3815 ++s; 3816 quot *= (top - s) / (_available - s); 3817 } 3818 } 3819 3820 return s; 3821 } 3822 3823 /* 3824 Randomly reset the value of _Vprime. 3825 */ 3826 private double newVprime(size_t remaining) 3827 { 3828 static if (is(UniformRNG == void)) 3829 { 3830 double r = uniform!"()"(0.0, 1.0); 3831 } 3832 else 3833 { 3834 double r = uniform!"()"(0.0, 1.0, _rng); 3835 } 3836 3837 return r ^^ (1.0 / remaining); 3838 } 3839 3840 /* 3841 Vitter's Algorithm D. For an extensive description of the algorithm 3842 and its rationale, see: 3843 3844 * Vitter, J.S. (1984), "Faster methods for random sampling", 3845 Commun. ACM 27(7): 703--718 3846 3847 * Vitter, J.S. (1987) "An efficient algorithm for sequential random 3848 sampling", ACM Trans. Math. Softw. 13(1): 58-67. 3849 3850 Variable names are chosen to match those in Vitter's paper. 3851 */ 3852 private size_t skipD() 3853 { 3854 import std.math.traits : isNaN; 3855 import std.math.rounding : trunc; 3856 // Confirm that the check in Step D1 is valid and we 3857 // haven't been sent here by mistake 3858 assert((_alphaInverse * _toSelect) <= _available); 3859 3860 // Now it's safe to use the standard Algorithm D mechanism. 3861 if (_toSelect > 1) 3862 { 3863 size_t s; 3864 size_t qu1 = 1 + _available - _toSelect; 3865 double x, y1; 3866 3867 assert(!_Vprime.isNaN()); 3868 3869 while (true) 3870 { 3871 // Step D2: set values of x and u. 3872 while (1) 3873 { 3874 x = _available * (1-_Vprime); 3875 s = cast(size_t) trunc(x); 3876 if (s < qu1) 3877 break; 3878 _Vprime = newVprime(_toSelect); 3879 } 3880 3881 static if (is(UniformRNG == void)) 3882 { 3883 double u = uniform!"()"(0.0, 1.0); 3884 } 3885 else 3886 { 3887 double u = uniform!"()"(0.0, 1.0, _rng); 3888 } 3889 3890 y1 = (u * (cast(double) _available) / qu1) ^^ (1.0/(_toSelect - 1)); 3891 3892 _Vprime = y1 * ((-x/_available)+1.0) * ( qu1/( (cast(double) qu1) - s ) ); 3893 3894 // Step D3: if _Vprime <= 1.0 our work is done and we return S. 3895 // Otherwise ... 3896 if (_Vprime > 1.0) 3897 { 3898 size_t top = _available - 1, limit; 3899 double y2 = 1.0, bottom; 3900 3901 if (_toSelect > (s+1)) 3902 { 3903 bottom = _available - _toSelect; 3904 limit = _available - s; 3905 } 3906 else 3907 { 3908 bottom = _available - (s+1); 3909 limit = qu1; 3910 } 3911 3912 foreach (size_t t; limit .. _available) 3913 { 3914 y2 *= top/bottom; 3915 top--; 3916 bottom--; 3917 } 3918 3919 // Step D4: decide whether or not to accept the current value of S. 3920 if (_available/(_available-x) < y1 * (y2 ^^ (1.0/(_toSelect-1)))) 3921 { 3922 // If it's not acceptable, we generate a new value of _Vprime 3923 // and go back to the start of the for (;;) loop. 3924 _Vprime = newVprime(_toSelect); 3925 } 3926 else 3927 { 3928 // If it's acceptable we generate a new value of _Vprime 3929 // based on the remaining number of sample points needed, 3930 // and return S. 3931 _Vprime = newVprime(_toSelect-1); 3932 return s; 3933 } 3934 } 3935 else 3936 { 3937 // Return if condition D3 satisfied. 3938 return s; 3939 } 3940 } 3941 } 3942 else 3943 { 3944 // If only one sample point remains to be taken ... 3945 return cast(size_t) trunc(_available * _Vprime); 3946 } 3947 } 3948 3949 private void prime() 3950 { 3951 if (empty) 3952 { 3953 return; 3954 } 3955 assert(_available && _available >= _toSelect); 3956 immutable size_t s = skip(); 3957 assert(s + _toSelect <= _available); 3958 static if (hasLength!Range) 3959 { 3960 assert(s + _toSelect <= _input.length); 3961 } 3962 assert(!_input.empty); 3963 _input.popFrontExactly(s); 3964 _index += s; 3965 _available -= s; 3966 assert(_available > 0); 3967 } 3968 } 3969 3970 /// Ditto 3971 auto randomSample(Range)(Range r, size_t n, size_t total) 3972 if (isInputRange!Range) 3973 { 3974 return RandomSample!(Range, void)(r, n, total); 3975 } 3976 3977 /// Ditto 3978 auto randomSample(Range)(Range r, size_t n) 3979 if (isInputRange!Range && hasLength!Range) 3980 { 3981 return RandomSample!(Range, void)(r, n, r.length); 3982 } 3983 3984 /// Ditto 3985 auto randomSample(Range, UniformRNG)(Range r, size_t n, size_t total, auto ref UniformRNG rng) 3986 if (isInputRange!Range && isUniformRNG!UniformRNG) 3987 { 3988 return RandomSample!(Range, UniformRNG)(r, n, total, rng); 3989 } 3990 3991 /// Ditto 3992 auto randomSample(Range, UniformRNG)(Range r, size_t n, auto ref UniformRNG rng) 3993 if (isInputRange!Range && hasLength!Range && isUniformRNG!UniformRNG) 3994 { 3995 return RandomSample!(Range, UniformRNG)(r, n, r.length, rng); 3996 } 3997 3998 /// 3999 @safe unittest 4000 { 4001 import std.algorithm.comparison : equal; 4002 import std.range : iota; 4003 auto rnd = MinstdRand0(42); 4004 assert(10.iota.randomSample(3, rnd).equal([7, 8, 9])); 4005 } 4006 4007 @system unittest 4008 { 4009 // @system because it takes the address of a local 4010 import std.conv : text; 4011 import std.exception; 4012 import std.range; 4013 // For test purposes, an infinite input range 4014 struct TestInputRange 4015 { 4016 private auto r = recurrence!"a[n-1] + 1"(0); 4017 bool empty() @property const pure nothrow { return r.empty; } 4018 auto front() @property pure nothrow { return r.front; } 4019 void popFront() pure nothrow { r.popFront(); } 4020 } 4021 static assert(isInputRange!TestInputRange); 4022 static assert(!isForwardRange!TestInputRange); 4023 4024 const(int)[] a = [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 ]; 4025 4026 foreach (UniformRNG; PseudoRngTypes) 4027 (){ // avoid workaround optimizations for large functions 4028 // https://issues.dlang.org/show_bug.cgi?id=2396 4029 auto rng = UniformRNG(1234); 4030 /* First test the most general case: randomSample of input range, with and 4031 * without a specified random number generator. 4032 */ 4033 static assert(isInputRange!(typeof(randomSample(TestInputRange(), 5, 10)))); 4034 static assert(isInputRange!(typeof(randomSample(TestInputRange(), 5, 10, rng)))); 4035 static assert(!isForwardRange!(typeof(randomSample(TestInputRange(), 5, 10)))); 4036 static assert(!isForwardRange!(typeof(randomSample(TestInputRange(), 5, 10, rng)))); 4037 // test case with range initialized by direct call to struct 4038 { 4039 auto sample = 4040 RandomSample!(TestInputRange, UniformRNG) 4041 (TestInputRange(), 5, 10, UniformRNG(987_654_321)); 4042 static assert(isInputRange!(typeof(sample))); 4043 static assert(!isForwardRange!(typeof(sample))); 4044 } 4045 4046 /* Now test the case of an input range with length. We ignore the cases 4047 * already covered by the previous tests. 4048 */ 4049 static assert(isInputRange!(typeof(randomSample(TestInputRange().takeExactly(10), 5)))); 4050 static assert(isInputRange!(typeof(randomSample(TestInputRange().takeExactly(10), 5, rng)))); 4051 static assert(!isForwardRange!(typeof(randomSample(TestInputRange().takeExactly(10), 5)))); 4052 static assert(!isForwardRange!(typeof(randomSample(TestInputRange().takeExactly(10), 5, rng)))); 4053 // test case with range initialized by direct call to struct 4054 { 4055 auto sample = 4056 RandomSample!(typeof(TestInputRange().takeExactly(10)), UniformRNG) 4057 (TestInputRange().takeExactly(10), 5, 10, UniformRNG(654_321_987)); 4058 static assert(isInputRange!(typeof(sample))); 4059 static assert(!isForwardRange!(typeof(sample))); 4060 } 4061 4062 // Now test the case of providing a forward range as input. 4063 static assert(!isForwardRange!(typeof(randomSample(a, 5)))); 4064 static if (isForwardRange!UniformRNG) 4065 { 4066 static assert(isForwardRange!(typeof(randomSample(a, 5, rng)))); 4067 // ... and test with range initialized directly 4068 { 4069 auto sample = 4070 RandomSample!(const(int)[], UniformRNG) 4071 (a, 5, UniformRNG(321_987_654)); 4072 static assert(isForwardRange!(typeof(sample))); 4073 } 4074 } 4075 else 4076 { 4077 static assert(isInputRange!(typeof(randomSample(a, 5, rng)))); 4078 static assert(!isForwardRange!(typeof(randomSample(a, 5, rng)))); 4079 // ... and test with range initialized directly 4080 { 4081 auto sample = 4082 RandomSample!(const(int)[], UniformRNG) 4083 (a, 5, UniformRNG(789_123_456)); 4084 static assert(isInputRange!(typeof(sample))); 4085 static assert(!isForwardRange!(typeof(sample))); 4086 } 4087 } 4088 4089 /* Check that randomSample will throw an error if we claim more 4090 * items are available than there actually are, or if we try to 4091 * sample more items than are available. */ 4092 assert(collectExceptionMsg( 4093 randomSample(a, 5, 15) 4094 ) == "RandomSample: specified 15 items as available when input contains only 10"); 4095 assert(collectExceptionMsg( 4096 randomSample(a, 15) 4097 ) == "RandomSample: cannot sample 15 items when only 10 are available"); 4098 assert(collectExceptionMsg( 4099 randomSample(a, 9, 8) 4100 ) == "RandomSample: cannot sample 9 items when only 8 are available"); 4101 assert(collectExceptionMsg( 4102 randomSample(TestInputRange(), 12, 11) 4103 ) == "RandomSample: cannot sample 12 items when only 11 are available"); 4104 4105 /* Check that sampling algorithm never accidentally overruns the end of 4106 * the input range. If input is an InputRange without .length, this 4107 * relies on the user specifying the total number of available items 4108 * correctly. 4109 */ 4110 { 4111 uint i = 0; 4112 foreach (e; randomSample(a, a.length)) 4113 { 4114 assert(e == i); 4115 ++i; 4116 } 4117 assert(i == a.length); 4118 4119 i = 0; 4120 foreach (e; randomSample(TestInputRange(), 17, 17)) 4121 { 4122 assert(e == i); 4123 ++i; 4124 } 4125 assert(i == 17); 4126 } 4127 4128 4129 // Check length properties of random samples. 4130 assert(randomSample(a, 5).length == 5); 4131 assert(randomSample(a, 5, 10).length == 5); 4132 assert(randomSample(a, 5, rng).length == 5); 4133 assert(randomSample(a, 5, 10, rng).length == 5); 4134 assert(randomSample(TestInputRange(), 5, 10).length == 5); 4135 assert(randomSample(TestInputRange(), 5, 10, rng).length == 5); 4136 4137 // ... and emptiness! 4138 assert(randomSample(a, 0).empty); 4139 assert(randomSample(a, 0, 5).empty); 4140 assert(randomSample(a, 0, rng).empty); 4141 assert(randomSample(a, 0, 5, rng).empty); 4142 assert(randomSample(TestInputRange(), 0, 10).empty); 4143 assert(randomSample(TestInputRange(), 0, 10, rng).empty); 4144 4145 /* Test that the (lazy) evaluation of random samples works correctly. 4146 * 4147 * We cover 2 different cases: a sample where the ratio of sample points 4148 * to total points is greater than the threshold for using Algorithm, and 4149 * one where the ratio is small enough (< 1/13) for Algorithm D to be used. 4150 * 4151 * For each, we also cover the case with and without a specified RNG. 4152 */ 4153 { 4154 // Small sample/source ratio, no specified RNG. 4155 uint i = 0; 4156 foreach (e; randomSample(randomCover(a), 5)) 4157 { 4158 ++i; 4159 } 4160 assert(i == 5); 4161 4162 // Small sample/source ratio, specified RNG. 4163 i = 0; 4164 foreach (e; randomSample(randomCover(a), 5, rng)) 4165 { 4166 ++i; 4167 } 4168 assert(i == 5); 4169 4170 // Large sample/source ratio, no specified RNG. 4171 i = 0; 4172 foreach (e; randomSample(TestInputRange(), 123, 123_456)) 4173 { 4174 ++i; 4175 } 4176 assert(i == 123); 4177 4178 // Large sample/source ratio, specified RNG. 4179 i = 0; 4180 foreach (e; randomSample(TestInputRange(), 123, 123_456, rng)) 4181 { 4182 ++i; 4183 } 4184 assert(i == 123); 4185 4186 /* Sample/source ratio large enough to start with Algorithm D, 4187 * small enough to switch to Algorithm A. 4188 */ 4189 i = 0; 4190 foreach (e; randomSample(TestInputRange(), 10, 131)) 4191 { 4192 ++i; 4193 } 4194 assert(i == 10); 4195 } 4196 4197 // Test that the .index property works correctly 4198 { 4199 auto sample1 = randomSample(TestInputRange(), 654, 654_321); 4200 for (; !sample1.empty; sample1.popFront()) 4201 { 4202 assert(sample1.front == sample1.index); 4203 } 4204 4205 auto sample2 = randomSample(TestInputRange(), 654, 654_321, rng); 4206 for (; !sample2.empty; sample2.popFront()) 4207 { 4208 assert(sample2.front == sample2.index); 4209 } 4210 4211 /* Check that it also works if .index is called before .front. 4212 * See: https://issues.dlang.org/show_bug.cgi?id=10322 4213 */ 4214 auto sample3 = randomSample(TestInputRange(), 654, 654_321); 4215 for (; !sample3.empty; sample3.popFront()) 4216 { 4217 assert(sample3.index == sample3.front); 4218 } 4219 4220 auto sample4 = randomSample(TestInputRange(), 654, 654_321, rng); 4221 for (; !sample4.empty; sample4.popFront()) 4222 { 4223 assert(sample4.index == sample4.front); 4224 } 4225 } 4226 4227 /* Test behaviour if .popFront() is called before sample is read. 4228 * This is a rough-and-ready check that the statistical properties 4229 * are in the ballpark -- not a proper validation of statistical 4230 * quality! This incidentally also checks for reference-type 4231 * initialization bugs, as the foreach () loop will operate on a 4232 * copy of the popFronted (and hence initialized) sample. 4233 */ 4234 { 4235 size_t count0, count1, count99; 4236 foreach (_; 0 .. 50_000) 4237 { 4238 auto sample = randomSample(iota(100), 5, &rng); 4239 sample.popFront(); 4240 foreach (s; sample) 4241 { 4242 if (s == 0) 4243 { 4244 ++count0; 4245 } 4246 else if (s == 1) 4247 { 4248 ++count1; 4249 } 4250 else if (s == 99) 4251 { 4252 ++count99; 4253 } 4254 } 4255 } 4256 /* Statistical assumptions here: this is a sequential sampling process 4257 * so (i) 0 can only be the first sample point, so _can't_ be in the 4258 * remainder of the sample after .popFront() is called. (ii) By similar 4259 * token, 1 can only be in the remainder if it's the 2nd point of the 4260 * whole sample, and hence if 0 was the first; probability of 0 being 4261 * first and 1 second is 5/100 * 4/99 (thank you, Algorithm S:-) and 4262 * so the mean count of 1 should be about 202. Finally, 99 can only 4263 * be the _last_ sample point to be picked, so its probability of 4264 * inclusion should be independent of the .popFront() and it should 4265 * occur with frequency 5/100, hence its count should be about 5000. 4266 * Unfortunately we have to set quite a high tolerance because with 4267 * sample size small enough for unittests to run in reasonable time, 4268 * the variance can be quite high. 4269 */ 4270 assert(count0 == 0); 4271 assert(count1 < 150, text("1: ", count1, " > 150.")); 4272 assert(2_200 < count99, text("99: ", count99, " < 2200.")); 4273 assert(count99 < 2_800, text("99: ", count99, " > 2800.")); 4274 } 4275 4276 /* Odd corner-cases: RandomSample has 2 constructors that are not called 4277 * by the randomSample() helper functions, but that can be used if the 4278 * constructor is called directly. These cover the case of the user 4279 * specifying input but not input length. 4280 */ 4281 { 4282 auto input1 = TestInputRange().takeExactly(456_789); 4283 static assert(hasLength!(typeof(input1))); 4284 auto sample1 = RandomSample!(typeof(input1), void)(input1, 789); 4285 static assert(isInputRange!(typeof(sample1))); 4286 static assert(!isForwardRange!(typeof(sample1))); 4287 assert(sample1.length == 789); 4288 assert(sample1._available == 456_789); 4289 uint i = 0; 4290 for (; !sample1.empty; sample1.popFront()) 4291 { 4292 assert(sample1.front == sample1.index); 4293 ++i; 4294 } 4295 assert(i == 789); 4296 4297 auto input2 = TestInputRange().takeExactly(456_789); 4298 static assert(hasLength!(typeof(input2))); 4299 auto sample2 = RandomSample!(typeof(input2), typeof(rng))(input2, 789, rng); 4300 static assert(isInputRange!(typeof(sample2))); 4301 static assert(!isForwardRange!(typeof(sample2))); 4302 assert(sample2.length == 789); 4303 assert(sample2._available == 456_789); 4304 i = 0; 4305 for (; !sample2.empty; sample2.popFront()) 4306 { 4307 assert(sample2.front == sample2.index); 4308 ++i; 4309 } 4310 assert(i == 789); 4311 } 4312 4313 /* Test that the save property works where input is a forward range, 4314 * and RandomSample is using a (forward range) random number generator 4315 * that is not rndGen. 4316 */ 4317 static if (isForwardRange!UniformRNG) 4318 { 4319 auto sample1 = randomSample(a, 5, rng); 4320 // https://issues.dlang.org/show_bug.cgi?id=15853 4321 auto sample2 = ((const ref typeof(sample1) a) => a.save)(sample1); 4322 assert(sample1.array() == sample2.array()); 4323 } 4324 4325 // https://issues.dlang.org/show_bug.cgi?id=8314 4326 { 4327 auto sample(RandomGen)(uint seed) { return randomSample(a, 1, RandomGen(seed)).front; } 4328 4329 // Start from 1 because not all RNGs accept 0 as seed. 4330 immutable fst = sample!UniformRNG(1); 4331 uint n = 1; 4332 while (sample!UniformRNG(++n) == fst && n < n.max) {} 4333 assert(n < n.max); 4334 } 4335 }(); 4336 }