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