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