1// Written in the D programming language. 2 3/** 4Facilities for random number generation. 5 6$(RED Disclaimer:) The _random number generators and API provided in this 7module are not designed to be cryptographically secure, and are therefore 8unsuitable for cryptographic or security-related purposes such as generating 9authentication tokens or network sequence numbers. For such needs, please use a 10reputable cryptographic library instead. 11 12The new-style generator objects hold their own state so they are 13immune of threading issues. The generators feature a number of 14well-known and well-documented methods of generating random 15numbers. An overall fast and reliable means to generate random numbers 16is the $(D_PARAM Mt19937) generator, which derives its name from 17"$(LINK2 https://en.wikipedia.org/wiki/Mersenne_Twister, Mersenne Twister) 18with a period of 2 to the power of 1919937". In memory-constrained situations, 20$(LINK2 https://en.wikipedia.org/wiki/Linear_congruential_generator, 21linear congruential generators) such as $(D MinstdRand0) and $(D MinstdRand) might be 22useful. The standard library provides an alias $(D_PARAM Random) for 23whichever generator it considers the most fit for the target 24environment. 25 26In addition to random number generators, this module features 27distributions, which skew a generator's output statistical 28distribution in various ways. So far the uniform distribution for 29integers and real numbers have been implemented. 30 31Source: $(PHOBOSSRC std/_random.d) 32 33Macros: 34 35Copyright: Copyright Andrei Alexandrescu 2008 - 2009, Joseph Rushton Wakeling 2012. 36License: $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0). 37Authors: $(HTTP erdani.org, Andrei Alexandrescu) 38 Masahiro Nakagawa (Xorshift random generator) 39 $(HTTP braingam.es, Joseph Rushton Wakeling) (Algorithm D for random sampling) 40 Ilya Yaroshenko (Mersenne Twister implementation, adapted from $(HTTPS github.com/libmir/mir-_random, mir-_random)) 41Credits: The entire random number library architecture is derived from the 42 excellent $(HTTP open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2461.pdf, C++0X) 43 random number facility proposed by Jens Maurer and contributed to by 44 researchers at the Fermi laboratory (excluding Xorshift). 45*/ 46/* 47 Copyright Andrei Alexandrescu 2008 - 2009. 48Distributed under the Boost Software License, Version 1.0. 49 (See accompanying file LICENSE_1_0.txt or copy at 50 http://www.boost.org/LICENSE_1_0.txt) 51*/ 52module std.random; 53 54 55import std.range.primitives; 56import std.traits; 57 58/// 59@safe unittest 60{ 61 // seed a random generator with a constant 62 auto rnd = Random(42); 63 64 // Generate a uniformly-distributed integer in the range [0, 14] 65 // If no random generator is passed, the global `rndGen` would be used 66 auto i = uniform(0, 15, rnd); 67 assert(i >= 0 && i < 15); 68 69 // Generate a uniformly-distributed real in the range [0, 100) 70 auto r = uniform(0.0L, 100.0L, rnd); 71 assert(r >= 0 && r < 100); 72 73 // Generate a 32-bit random number 74 auto u = uniform!uint(rnd); 75 static assert(is(typeof(u) == uint)); 76} 77 78version (unittest) 79{ 80 static import std.meta; 81 package alias PseudoRngTypes = std.meta.AliasSeq!(MinstdRand0, MinstdRand, Mt19937, Xorshift32, Xorshift64, 82 Xorshift96, Xorshift128, Xorshift160, Xorshift192); 83} 84 85// Segments of the code in this file Copyright (c) 1997 by Rick Booth 86// From "Inner Loops" by Rick Booth, Addison-Wesley 87 88// Work derived from: 89 90/* 91 A C-program for MT19937, with initialization improved 2002/1/26. 92 Coded by Takuji Nishimura and Makoto Matsumoto. 93 94 Before using, initialize the state by using init_genrand(seed) 95 or init_by_array(init_key, key_length). 96 97 Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, 98 All rights reserved. 99 100 Redistribution and use in source and binary forms, with or without 101 modification, are permitted provided that the following conditions 102 are met: 103 104 1. Redistributions of source code must retain the above copyright 105 notice, this list of conditions and the following disclaimer. 106 107 2. Redistributions in binary form must reproduce the above copyright 108 notice, this list of conditions and the following disclaimer in the 109 documentation and/or other materials provided with the distribution. 110 111 3. The names of its contributors may not be used to endorse or promote 112 products derived from this software without specific prior written 113 permission. 114 115 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 116 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 117 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 118 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 119 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 120 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 121 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 122 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 123 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 124 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 125 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 126 127 128 Any feedback is very welcome. 129 http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html 130 email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) 131*/ 132 133/** 134 * Test if Rng is a random-number generator. The overload 135 * taking a ElementType also makes sure that the Rng generates 136 * values of that type. 137 * 138 * A random-number generator has at least the following features: 139 * $(UL 140 * $(LI it's an InputRange) 141 * $(LI it has a 'bool isUniformRandom' field readable in CTFE) 142 * ) 143 */ 144template isUniformRNG(Rng, ElementType) 145{ 146 enum bool isUniformRNG = isInputRange!Rng && 147 is(typeof(Rng.front) == ElementType) && 148 is(typeof( 149 { 150 static assert(Rng.isUniformRandom); //tag 151 })); 152} 153 154/** 155 * ditto 156 */ 157template isUniformRNG(Rng) 158{ 159 enum bool isUniformRNG = isInputRange!Rng && 160 is(typeof( 161 { 162 static assert(Rng.isUniformRandom); //tag 163 })); 164} 165 166/** 167 * Test if Rng is seedable. The overload 168 * taking a SeedType also makes sure that the Rng can be seeded with SeedType. 169 * 170 * A seedable random-number generator has the following additional features: 171 * $(UL 172 * $(LI it has a 'seed(ElementType)' function) 173 * ) 174 */ 175template isSeedable(Rng, SeedType) 176{ 177 enum bool isSeedable = isUniformRNG!(Rng) && 178 is(typeof( 179 { 180 Rng r = void; // can define a Rng object 181 r.seed(SeedType.init); // can seed a Rng 182 })); 183} 184 185///ditto 186template isSeedable(Rng) 187{ 188 enum bool isSeedable = isUniformRNG!Rng && 189 is(typeof( 190 { 191 Rng r = void; // can define a Rng object 192 r.seed(typeof(r.front).init); // can seed a Rng 193 })); 194} 195 196@safe pure nothrow unittest 197{ 198 struct NoRng 199 { 200 @property uint front() {return 0;} 201 @property bool empty() {return false;} 202 void popFront() {} 203 } 204 static assert(!isUniformRNG!(NoRng, uint)); 205 static assert(!isUniformRNG!(NoRng)); 206 static assert(!isSeedable!(NoRng, uint)); 207 static assert(!isSeedable!(NoRng)); 208 209 struct NoRng2 210 { 211 @property uint front() {return 0;} 212 @property bool empty() {return false;} 213 void popFront() {} 214 215 enum isUniformRandom = false; 216 } 217 static assert(!isUniformRNG!(NoRng2, uint)); 218 static assert(!isUniformRNG!(NoRng2)); 219 static assert(!isSeedable!(NoRng2, uint)); 220 static assert(!isSeedable!(NoRng2)); 221 222 struct NoRng3 223 { 224 @property bool empty() {return false;} 225 void popFront() {} 226 227 enum isUniformRandom = true; 228 } 229 static assert(!isUniformRNG!(NoRng3, uint)); 230 static assert(!isUniformRNG!(NoRng3)); 231 static assert(!isSeedable!(NoRng3, uint)); 232 static assert(!isSeedable!(NoRng3)); 233 234 struct validRng 235 { 236 @property uint front() {return 0;} 237 @property bool empty() {return false;} 238 void popFront() {} 239 240 enum isUniformRandom = true; 241 } 242 static assert(isUniformRNG!(validRng, uint)); 243 static assert(isUniformRNG!(validRng)); 244 static assert(!isSeedable!(validRng, uint)); 245 static assert(!isSeedable!(validRng)); 246 247 struct seedRng 248 { 249 @property uint front() {return 0;} 250 @property bool empty() {return false;} 251 void popFront() {} 252 void seed(uint val){} 253 enum isUniformRandom = true; 254 } 255 static assert(isUniformRNG!(seedRng, uint)); 256 static assert(isUniformRNG!(seedRng)); 257 static assert(isSeedable!(seedRng, uint)); 258 static assert(isSeedable!(seedRng)); 259} 260 261/** 262Linear Congruential generator. 263 */ 264struct LinearCongruentialEngine(UIntType, UIntType a, UIntType c, UIntType m) 265if (isUnsigned!UIntType) 266{ 267 ///Mark this as a Rng 268 enum bool isUniformRandom = true; 269 /// Does this generator have a fixed range? ($(D_PARAM true)). 270 enum bool hasFixedRange = true; 271 /// Lowest generated value ($(D 1) if $(D c == 0), $(D 0) otherwise). 272 enum UIntType min = ( c == 0 ? 1 : 0 ); 273 /// Highest generated value ($(D modulus - 1)). 274 enum UIntType max = m - 1; 275/** 276The parameters of this distribution. The random number is $(D_PARAM x 277= (x * multipler + increment) % modulus). 278 */ 279 enum UIntType multiplier = a; 280 ///ditto 281 enum UIntType increment = c; 282 ///ditto 283 enum UIntType modulus = m; 284 285 static assert(isIntegral!(UIntType)); 286 static assert(m == 0 || a < m); 287 static assert(m == 0 || c < m); 288 static assert(m == 0 || 289 (cast(ulong) a * (m-1) + c) % m == (c < a ? c - a + m : c - a)); 290 291 // Check for maximum range 292 private static ulong gcd(ulong a, ulong b) @safe pure nothrow @nogc 293 { 294 while (b) 295 { 296 auto t = b; 297 b = a % b; 298 a = t; 299 } 300 return a; 301 } 302 303 private static ulong primeFactorsOnly(ulong n) @safe pure nothrow @nogc 304 { 305 ulong result = 1; 306 ulong iter = 2; 307 for (; n >= iter * iter; iter += 2 - (iter == 2)) 308 { 309 if (n % iter) continue; 310 result *= iter; 311 do 312 { 313 n /= iter; 314 } while (n % iter == 0); 315 } 316 return result * n; 317 } 318 319 @safe pure nothrow unittest 320 { 321 static assert(primeFactorsOnly(100) == 10); 322 //writeln(primeFactorsOnly(11)); 323 static assert(primeFactorsOnly(11) == 11); 324 static assert(primeFactorsOnly(7 * 7 * 7 * 11 * 15 * 11) == 7 * 11 * 15); 325 static assert(primeFactorsOnly(129 * 2) == 129 * 2); 326 // enum x = primeFactorsOnly(7 * 7 * 7 * 11 * 15); 327 // static assert(x == 7 * 11 * 15); 328 } 329 330 private static bool properLinearCongruentialParameters(ulong m, 331 ulong a, ulong c) @safe pure nothrow @nogc 332 { 333 if (m == 0) 334 { 335 static if (is(UIntType == uint)) 336 { 337 // Assume m is uint.max + 1 338 m = (1uL << 32); 339 } 340 else 341 { 342 return false; 343 } 344 } 345 // Bounds checking 346 if (a == 0 || a >= m || c >= m) return false; 347 // c and m are relatively prime 348 if (c > 0 && gcd(c, m) != 1) return false; 349 // a - 1 is divisible by all prime factors of m 350 if ((a - 1) % primeFactorsOnly(m)) return false; 351 // if a - 1 is multiple of 4, then m is a multiple of 4 too. 352 if ((a - 1) % 4 == 0 && m % 4) return false; 353 // Passed all tests 354 return true; 355 } 356 357 // check here 358 static assert(c == 0 || properLinearCongruentialParameters(m, a, c), 359 "Incorrect instantiation of LinearCongruentialEngine"); 360 361/** 362Constructs a $(D_PARAM LinearCongruentialEngine) generator seeded with 363$(D x0). 364 */ 365 this(UIntType x0) @safe pure 366 { 367 seed(x0); 368 } 369 370/** 371 (Re)seeds the generator. 372*/ 373 void seed(UIntType x0 = 1) @safe pure 374 { 375 static if (c == 0) 376 { 377 import std.exception : enforce; 378 enforce(x0, "Invalid (zero) seed for " 379 ~ LinearCongruentialEngine.stringof); 380 } 381 _x = modulus ? (x0 % modulus) : x0; 382 popFront(); 383 } 384 385/** 386 Advances the random sequence. 387*/ 388 void popFront() @safe pure nothrow @nogc 389 { 390 static if (m) 391 { 392 static if (is(UIntType == uint) && m == uint.max) 393 { 394 immutable ulong 395 x = (cast(ulong) a * _x + c), 396 v = x >> 32, 397 w = x & uint.max; 398 immutable y = cast(uint)(v + w); 399 _x = (y < v || y == uint.max) ? (y + 1) : y; 400 } 401 else static if (is(UIntType == uint) && m == int.max) 402 { 403 immutable ulong 404 x = (cast(ulong) a * _x + c), 405 v = x >> 31, 406 w = x & int.max; 407 immutable uint y = cast(uint)(v + w); 408 _x = (y >= int.max) ? (y - int.max) : y; 409 } 410 else 411 { 412 _x = cast(UIntType) ((cast(ulong) a * _x + c) % m); 413 } 414 } 415 else 416 { 417 _x = a * _x + c; 418 } 419 } 420 421/** 422 Returns the current number in the random sequence. 423*/ 424 @property UIntType front() const @safe pure nothrow @nogc 425 { 426 return _x; 427 } 428 429/// 430 @property typeof(this) save() @safe pure nothrow @nogc 431 { 432 return this; 433 } 434 435/** 436Always $(D false) (random generators are infinite ranges). 437 */ 438 enum bool empty = false; 439 440/** 441 Compares against $(D_PARAM rhs) for equality. 442 */ 443 bool opEquals(ref const LinearCongruentialEngine rhs) const @safe pure nothrow @nogc 444 { 445 return _x == rhs._x; 446 } 447 448 private UIntType _x = m ? (a + c) % m : (a + c); 449} 450 451/** 452Define $(D_PARAM LinearCongruentialEngine) generators with well-chosen 453parameters. $(D MinstdRand0) implements Park and Miller's "minimal 454standard" $(HTTP 455wikipedia.org/wiki/Park%E2%80%93Miller_random_number_generator, 456generator) that uses 16807 for the multiplier. $(D MinstdRand) 457implements a variant that has slightly better spectral behavior by 458using the multiplier 48271. Both generators are rather simplistic. 459 */ 460alias MinstdRand0 = LinearCongruentialEngine!(uint, 16_807, 0, 2_147_483_647); 461/// ditto 462alias MinstdRand = LinearCongruentialEngine!(uint, 48_271, 0, 2_147_483_647); 463 464/// 465@safe unittest 466{ 467 // seed with a constant 468 auto rnd0 = MinstdRand0(1); 469 auto n = rnd0.front; // same for each run 470 // Seed with an unpredictable value 471 rnd0.seed(unpredictableSeed); 472 n = rnd0.front; // different across runs 473} 474 475@safe unittest 476{ 477 import std.range; 478 static assert(isForwardRange!MinstdRand); 479 static assert(isUniformRNG!MinstdRand); 480 static assert(isUniformRNG!MinstdRand0); 481 static assert(isUniformRNG!(MinstdRand, uint)); 482 static assert(isUniformRNG!(MinstdRand0, uint)); 483 static assert(isSeedable!MinstdRand); 484 static assert(isSeedable!MinstdRand0); 485 static assert(isSeedable!(MinstdRand, uint)); 486 static assert(isSeedable!(MinstdRand0, uint)); 487 488 // The correct numbers are taken from The Database of Integer Sequences 489 // http://www.research.att.com/~njas/sequences/eisBTfry00128.txt 490 auto checking0 = [ 491 16807UL,282475249,1622650073,984943658,1144108930,470211272, 492 101027544,1457850878,1458777923,2007237709,823564440,1115438165, 493 1784484492,74243042,114807987,1137522503,1441282327,16531729, 494 823378840,143542612 ]; 495 //auto rnd0 = MinstdRand0(1); 496 MinstdRand0 rnd0; 497 498 foreach (e; checking0) 499 { 500 assert(rnd0.front == e); 501 rnd0.popFront(); 502 } 503 // Test the 10000th invocation 504 // Correct value taken from: 505 // http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2461.pdf 506 rnd0.seed(); 507 popFrontN(rnd0, 9999); 508 assert(rnd0.front == 1043618065); 509 510 // Test MinstdRand 511 auto checking = [48271UL,182605794,1291394886,1914720637,2078669041, 512 407355683]; 513 //auto rnd = MinstdRand(1); 514 MinstdRand rnd; 515 foreach (e; checking) 516 { 517 assert(rnd.front == e); 518 rnd.popFront(); 519 } 520 521 // Test the 10000th invocation 522 // Correct value taken from: 523 // http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2461.pdf 524 rnd.seed(); 525 popFrontN(rnd, 9999); 526 assert(rnd.front == 399268537); 527 528 // Check .save works 529 foreach (Type; std.meta.AliasSeq!(MinstdRand0, MinstdRand)) 530 { 531 auto rnd1 = Type(unpredictableSeed); 532 auto rnd2 = rnd1.save; 533 assert(rnd1 == rnd2); 534 // Enable next test when RNGs are reference types 535 version (none) { assert(rnd1 !is rnd2); } 536 assert(rnd1.take(100).array() == rnd2.take(100).array()); 537 } 538} 539 540/** 541The $(LINK2 https://en.wikipedia.org/wiki/Mersenne_Twister, Mersenne Twister) generator. 542 */ 543struct MersenneTwisterEngine(UIntType, size_t w, size_t n, size_t m, size_t r, 544 UIntType a, size_t u, UIntType d, size_t s, 545 UIntType b, size_t t, 546 UIntType c, size_t l, UIntType f) 547if (isUnsigned!UIntType) 548{ 549 static assert(0 < w && w <= UIntType.sizeof * 8); 550 static assert(1 <= m && m <= n); 551 static assert(0 <= r && 0 <= u && 0 <= s && 0 <= t && 0 <= l); 552 static assert(r <= w && u <= w && s <= w && t <= w && l <= w); 553 static assert(0 <= a && 0 <= b && 0 <= c); 554 static assert(n <= sizediff_t.max); 555 556 ///Mark this as a Rng 557 enum bool isUniformRandom = true; 558 559/** 560Parameters for the generator. 561*/ 562 enum size_t wordSize = w; 563 enum size_t stateSize = n; /// ditto 564 enum size_t shiftSize = m; /// ditto 565 enum size_t maskBits = r; /// ditto 566 enum UIntType xorMask = a; /// ditto 567 enum size_t temperingU = u; /// ditto 568 enum UIntType temperingD = d; /// ditto 569 enum size_t temperingS = s; /// ditto 570 enum UIntType temperingB = b; /// ditto 571 enum size_t temperingT = t; /// ditto 572 enum UIntType temperingC = c; /// ditto 573 enum size_t temperingL = l; /// ditto 574 enum UIntType initializationMultiplier = f; /// ditto 575 576 /// Smallest generated value (0). 577 enum UIntType min = 0; 578 /// Largest generated value. 579 enum UIntType max = UIntType.max >> (UIntType.sizeof * 8u - w); 580 // note, `max` also serves as a bitmask for the lowest `w` bits 581 static assert(a <= max && b <= max && c <= max && f <= max); 582 583 /// The default seed value. 584 enum UIntType defaultSeed = 5489u; 585 586 // Bitmasks used in the 'twist' part of the algorithm 587 private enum UIntType lowerMask = (cast(UIntType) 1u << r) - 1, 588 upperMask = (~lowerMask) & this.max; 589 590 /* 591 Collection of all state variables 592 used by the generator 593 */ 594 private struct State 595 { 596 /* 597 State array of the generator. This 598 is iterated through backwards (from 599 last element to first), providing a 600 few extra compiler optimizations by 601 comparison to the forward iteration 602 used in most implementations. 603 */ 604 UIntType[n] data; 605 606 /* 607 Cached copy of most recently updated 608 element of `data` state array, ready 609 to be tempered to generate next 610 `front` value 611 */ 612 UIntType z; 613 614 /* 615 Most recently generated random variate 616 */ 617 UIntType front; 618 619 /* 620 Index of the entry in the `data` 621 state array that will be twisted 622 in the next `popFront()` call 623 */ 624 size_t index; 625 } 626 627 /* 628 State variables used by the generator; 629 initialized to values equivalent to 630 explicitly seeding the generator with 631 `defaultSeed` 632 */ 633 private State state = defaultState(); 634 /* NOTE: the above is a workaround to ensure 635 backwards compatibility with the original 636 implementation, which permitted implicit 637 construction. With `@disable this();` 638 it would not be necessary. */ 639 640/** 641 Constructs a MersenneTwisterEngine object. 642*/ 643 this(UIntType value) @safe pure nothrow @nogc 644 { 645 seed(value); 646 } 647 648 /** 649 Generates the default initial state for a Mersenne 650 Twister; equivalent to the internal state obtained 651 by calling `seed(defaultSeed)` 652 */ 653 private static State defaultState() @safe pure nothrow @nogc 654 { 655 if (!__ctfe) assert(false); 656 State mtState; 657 seedImpl(defaultSeed, mtState); 658 return mtState; 659 } 660 661/** 662 Seeds a MersenneTwisterEngine object. 663 Note: 664 This seed function gives 2^w starting points (the lowest w bits of 665 the value provided will be used). To allow the RNG to be started 666 in any one of its internal states use the seed overload taking an 667 InputRange. 668*/ 669 void seed()(UIntType value = defaultSeed) @safe pure nothrow @nogc 670 { 671 this.seedImpl(value, this.state); 672 } 673 674 /** 675 Implementation of the seeding mechanism, which 676 can be used with an arbitrary `State` instance 677 */ 678 private static void seedImpl(UIntType value, ref State mtState) 679 { 680 mtState.data[$ - 1] = value; 681 static if (this.max != UIntType.max) 682 { 683 mtState.data[$ - 1] &= this.max; 684 } 685 686 foreach_reverse (size_t i, ref e; mtState.data[0 .. $ - 1]) 687 { 688 e = f * (mtState.data[i + 1] ^ (mtState.data[i + 1] >> (w - 2))) + cast(UIntType)(n - (i + 1)); 689 static if (this.max != UIntType.max) 690 { 691 e &= this.max; 692 } 693 } 694 695 mtState.index = n - 1; 696 697 /* double popFront() to guarantee both `mtState.z` 698 and `mtState.front` are derived from the newly 699 set values in `mtState.data` */ 700 MersenneTwisterEngine.popFrontImpl(mtState); 701 MersenneTwisterEngine.popFrontImpl(mtState); 702 } 703 704/** 705 Seeds a MersenneTwisterEngine object using an InputRange. 706 707 Throws: 708 $(D Exception) if the InputRange didn't provide enough elements to seed the generator. 709 The number of elements required is the 'n' template parameter of the MersenneTwisterEngine struct. 710 */ 711 void seed(T)(T range) if (isInputRange!T && is(Unqual!(ElementType!T) == UIntType)) 712 { 713 this.seedImpl(range, this.state); 714 } 715 716 /** 717 Implementation of the range-based seeding mechanism, 718 which can be used with an arbitrary `State` instance 719 */ 720 private static void seedImpl(T)(T range, ref State mtState) 721 if (isInputRange!T && is(Unqual!(ElementType!T) == UIntType)) 722 { 723 size_t j; 724 for (j = 0; j < n && !range.empty; ++j, range.popFront()) 725 { 726 sizediff_t idx = n - j - 1; 727 mtState.data[idx] = range.front; 728 } 729 730 mtState.index = n - 1; 731 732 if (range.empty && j < n) 733 { 734 import core.internal.string : UnsignedStringBuf, unsignedToTempString; 735 736 UnsignedStringBuf buf = void; 737 string s = "MersenneTwisterEngine.seed: Input range didn't provide enough elements: Need "; 738 s ~= unsignedToTempString(n, buf, 10) ~ " elements."; 739 throw new Exception(s); 740 } 741 742 /* double popFront() to guarantee both `mtState.z` 743 and `mtState.front` are derived from the newly 744 set values in `mtState.data` */ 745 MersenneTwisterEngine.popFrontImpl(mtState); 746 MersenneTwisterEngine.popFrontImpl(mtState); 747 } 748 749/** 750 Advances the generator. 751*/ 752 void popFront() @safe pure nothrow @nogc 753 { 754 this.popFrontImpl(this.state); 755 } 756 757 /* 758 Internal implementation of `popFront()`, which 759 can be used with an arbitrary `State` instance 760 */ 761 private static void popFrontImpl(ref State mtState) 762 { 763 /* This function blends two nominally independent 764 processes: (i) calculation of the next random 765 variate `mtState.front` from the cached previous 766 `data` entry `z`, and (ii) updating the value 767 of `data[index]` and `mtState.z` and advancing 768 the `index` value to the next in sequence. 769 770 By interweaving the steps involved in these 771 procedures, rather than performing each of 772 them separately in sequence, the variables 773 are kept 'hot' in CPU registers, allowing 774 for significantly faster performance. */ 775 sizediff_t index = mtState.index; 776 sizediff_t next = index - 1; 777 if (next < 0) 778 next = n - 1; 779 auto z = mtState.z; 780 sizediff_t conj = index - m; 781 if (conj < 0) 782 conj = index - m + n; 783 784 static if (d == UIntType.max) 785 { 786 z ^= (z >> u); 787 } 788 else 789 { 790 z ^= (z >> u) & d; 791 } 792 793 auto q = mtState.data[index] & upperMask; 794 auto p = mtState.data[next] & lowerMask; 795 z ^= (z << s) & b; 796 auto y = q | p; 797 auto x = y >> 1; 798 z ^= (z << t) & c; 799 if (y & 1) 800 x ^= a; 801 auto e = mtState.data[conj] ^ x; 802 z ^= (z >> l); 803 mtState.z = mtState.data[index] = e; 804 mtState.index = next; 805 806 /* technically we should take the lowest `w` 807 bits here, but if the tempering bitmasks 808 `b` and `c` are set correctly, this should 809 be unnecessary */ 810 mtState.front = z; 811 } 812 813/** 814 Returns the current random value. 815 */ 816 @property UIntType front() @safe const pure nothrow @nogc 817 { 818 return this.state.front; 819 } 820 821/// 822 @property typeof(this) save() @safe pure nothrow @nogc 823 { 824 return this; 825 } 826 827/** 828Always $(D false). 829 */ 830 enum bool empty = false; 831} 832 833/** 834A $(D MersenneTwisterEngine) instantiated with the parameters of the 835original engine $(HTTP math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html, 836MT19937), generating uniformly-distributed 32-bit numbers with a 837period of 2 to the power of 19937. Recommended for random number 838generation unless memory is severely restricted, in which case a $(D 839LinearCongruentialEngine) would be the generator of choice. 840 */ 841alias Mt19937 = MersenneTwisterEngine!(uint, 32, 624, 397, 31, 842 0x9908b0df, 11, 0xffffffff, 7, 843 0x9d2c5680, 15, 844 0xefc60000, 18, 1_812_433_253); 845 846/// 847@safe unittest 848{ 849 // seed with a constant 850 Mt19937 gen; 851 auto n = gen.front; // same for each run 852 // Seed with an unpredictable value 853 gen.seed(unpredictableSeed); 854 n = gen.front; // different across runs 855} 856 857@safe nothrow unittest 858{ 859 import std.algorithm; 860 import std.range; 861 static assert(isUniformRNG!Mt19937); 862 static assert(isUniformRNG!(Mt19937, uint)); 863 static assert(isSeedable!Mt19937); 864 static assert(isSeedable!(Mt19937, uint)); 865 static assert(isSeedable!(Mt19937, typeof(map!((a) => unpredictableSeed)(repeat(0))))); 866 Mt19937 gen; 867 assert(gen.front == 3499211612); 868 popFrontN(gen, 9999); 869 assert(gen.front == 4123659995); 870 try { gen.seed(iota(624u)); } catch (Exception) { assert(false); } 871 assert(gen.front == 3708921088u); 872 popFrontN(gen, 9999); 873 assert(gen.front == 165737292u); 874} 875 876/** 877A $(D MersenneTwisterEngine) instantiated with the parameters of the 878original engine $(HTTP en.wikipedia.org/wiki/Mersenne_Twister, 879MT19937-64), generating uniformly-distributed 64-bit numbers with a 880period of 2 to the power of 19937. 881*/ 882alias Mt19937_64 = MersenneTwisterEngine!(ulong, 64, 312, 156, 31, 883 0xb5026f5aa96619e9, 29, 0x5555555555555555, 17, 884 0x71d67fffeda60000, 37, 885 0xfff7eee000000000, 43, 6_364_136_223_846_793_005); 886 887/// 888@safe unittest 889{ 890 // Seed with a constant 891 auto gen = Mt19937_64(12345); 892 auto n = gen.front; // same for each run 893 // Seed with an unpredictable value 894 gen.seed(unpredictableSeed); 895 n = gen.front; // different across runs 896} 897 898@safe nothrow unittest 899{ 900 import std.algorithm; 901 import std.range; 902 static assert(isUniformRNG!Mt19937_64); 903 static assert(isUniformRNG!(Mt19937_64, ulong)); 904 static assert(isSeedable!Mt19937_64); 905 static assert(isSeedable!(Mt19937_64, ulong)); 906 // FIXME: this test demonstrates viably that Mt19937_64 907 // is seedable with an infinite range of `ulong` values 908 // but it's a poor example of how to actually seed the 909 // generator, since it can't cover the full range of 910 // possible seed values. Ideally we need a 64-bit 911 // unpredictable seed to complement the 32-bit one! 912 static assert(isSeedable!(Mt19937_64, typeof(map!((a) => (cast(ulong) unpredictableSeed))(repeat(0))))); 913 Mt19937_64 gen; 914 assert(gen.front == 14514284786278117030uL); 915 popFrontN(gen, 9999); 916 assert(gen.front == 9981545732273789042uL); 917 try { gen.seed(iota(312uL)); } catch (Exception) { assert(false); } 918 assert(gen.front == 14660652410669508483uL); 919 popFrontN(gen, 9999); 920 assert(gen.front == 15956361063660440239uL); 921} 922 923@safe unittest 924{ 925 import std.algorithm; 926 import std.exception; 927 import std.range; 928 929 Mt19937 gen; 930 931 assertThrown(gen.seed(map!((a) => unpredictableSeed)(repeat(0, 623)))); 932 933 gen.seed(map!((a) => unpredictableSeed)(repeat(0, 624))); 934 //infinite Range 935 gen.seed(map!((a) => unpredictableSeed)(repeat(0))); 936} 937 938@safe pure nothrow unittest 939{ 940 uint a, b; 941 { 942 Mt19937 gen; 943 a = gen.front; 944 } 945 { 946 Mt19937 gen; 947 gen.popFront(); 948 //popFrontN(gen, 1); // skip 1 element 949 b = gen.front; 950 } 951 assert(a != b); 952} 953 954@safe unittest 955{ 956 import std.range; 957 // Check .save works 958 foreach (Type; std.meta.AliasSeq!(Mt19937, Mt19937_64)) 959 { 960 auto gen1 = Type(unpredictableSeed); 961 auto gen2 = gen1.save; 962 assert(gen1 == gen2); // Danger, Will Robinson -- no opEquals for MT 963 // Enable next test when RNGs are reference types 964 version (none) { assert(gen1 !is gen2); } 965 assert(gen1.take(100).array() == gen2.take(100).array()); 966 } 967} 968 969@safe pure nothrow unittest //11690 970{ 971 alias MT(UIntType, uint w) = MersenneTwisterEngine!(UIntType, w, 624, 397, 31, 972 0x9908b0df, 11, 0xffffffff, 7, 973 0x9d2c5680, 15, 974 0xefc60000, 18, 1812433253); 975 976 ulong[] expectedFirstValue = [3499211612uL, 3499211612uL, 977 171143175841277uL, 1145028863177033374uL]; 978 979 ulong[] expected10kValue = [4123659995uL, 4123659995uL, 980 51991688252792uL, 3031481165133029945uL]; 981 982 foreach (i, R; std.meta.AliasSeq!(MT!(uint, 32), MT!(ulong, 32), MT!(ulong, 48), MT!(ulong, 64))) 983 { 984 auto a = R(); 985 a.seed(a.defaultSeed); // checks that some alternative paths in `seed` are utilized 986 assert(a.front == expectedFirstValue[i]); 987 a.popFrontN(9999); 988 assert(a.front == expected10kValue[i]); 989 } 990} 991 992 993/** 994 * Xorshift generator using 32bit algorithm. 995 * 996 * Implemented according to $(HTTP www.jstatsoft.org/v08/i14/paper, Xorshift RNGs). 997 * Supporting bits are below, $(D bits) means second parameter of XorshiftEngine. 998 * 999 * $(BOOKTABLE , 1000 * $(TR $(TH bits) $(TH period)) 1001 * $(TR $(TD 32) $(TD 2^32 - 1)) 1002 * $(TR $(TD 64) $(TD 2^64 - 1)) 1003 * $(TR $(TD 96) $(TD 2^96 - 1)) 1004 * $(TR $(TD 128) $(TD 2^128 - 1)) 1005 * $(TR $(TD 160) $(TD 2^160 - 1)) 1006 * $(TR $(TD 192) $(TD 2^192 - 2^32)) 1007 * ) 1008 */ 1009struct XorshiftEngine(UIntType, UIntType bits, UIntType a, UIntType b, UIntType c) 1010if (isUnsigned!UIntType) 1011{ 1012 static assert(bits == 32 || bits == 64 || bits == 96 || bits == 128 || bits == 160 || bits == 192, 1013 "Xorshift supports only 32, 64, 96, 128, 160 and 192 bit versions. " 1014 ~ to!string(bits) ~ " is not supported."); 1015 1016 public: 1017 ///Mark this as a Rng 1018 enum bool isUniformRandom = true; 1019 /// Always $(D false) (random generators are infinite ranges). 1020 enum empty = false; 1021 /// Smallest generated value. 1022 enum UIntType min = 0; 1023 /// Largest generated value. 1024 enum UIntType max = UIntType.max; 1025 1026 1027 private: 1028 enum size = bits / 32; 1029 1030 static if (bits == 32) 1031 UIntType[size] seeds_ = [2_463_534_242]; 1032 else static if (bits == 64) 1033 UIntType[size] seeds_ = [123_456_789, 362_436_069]; 1034 else static if (bits == 96) 1035 UIntType[size] seeds_ = [123_456_789, 362_436_069, 521_288_629]; 1036 else static if (bits == 128) 1037 UIntType[size] seeds_ = [123_456_789, 362_436_069, 521_288_629, 88_675_123]; 1038 else static if (bits == 160) 1039 UIntType[size] seeds_ = [123_456_789, 362_436_069, 521_288_629, 88_675_123, 5_783_321]; 1040 else static if (bits == 192) 1041 { 1042 UIntType[size] seeds_ = [123_456_789, 362_436_069, 521_288_629, 88_675_123, 5_783_321, 6_615_241]; 1043 UIntType value_; 1044 } 1045 else 1046 { 1047 static assert(false, "Phobos Error: Xorshift has no instantiation rule for " 1048 ~ to!string(bits) ~ " bits."); 1049 } 1050 1051 1052 public: 1053 /** 1054 * Constructs a $(D XorshiftEngine) generator seeded with $(D_PARAM x0). 1055 */ 1056 this(UIntType x0) @safe pure nothrow @nogc 1057 { 1058 seed(x0); 1059 } 1060 1061 1062 /** 1063 * (Re)seeds the generator. 1064 */ 1065 void seed(UIntType x0) @safe pure nothrow @nogc 1066 { 1067 // Initialization routine from MersenneTwisterEngine. 1068 foreach (i, e; seeds_) 1069 seeds_[i] = x0 = cast(UIntType)(1_812_433_253U * (x0 ^ (x0 >> 30)) + i + 1); 1070 1071 // All seeds must not be 0. 1072 sanitizeSeeds(seeds_); 1073 1074 popFront(); 1075 } 1076 1077 1078 /** 1079 * Returns the current number in the random sequence. 1080 */ 1081 @property 1082 UIntType front() const @safe pure nothrow @nogc 1083 { 1084 static if (bits == 192) 1085 return value_; 1086 else 1087 return seeds_[size - 1]; 1088 } 1089 1090 1091 /** 1092 * Advances the random sequence. 1093 */ 1094 void popFront() @safe pure nothrow @nogc 1095 { 1096 UIntType temp; 1097 1098 static if (bits == 32) 1099 { 1100 temp = seeds_[0] ^ (seeds_[0] << a); 1101 temp = temp ^ (temp >> b); 1102 seeds_[0] = temp ^ (temp << c); 1103 } 1104 else static if (bits == 64) 1105 { 1106 temp = seeds_[0] ^ (seeds_[0] << a); 1107 seeds_[0] = seeds_[1]; 1108 seeds_[1] = seeds_[1] ^ (seeds_[1] >> c) ^ temp ^ (temp >> b); 1109 } 1110 else static if (bits == 96) 1111 { 1112 temp = seeds_[0] ^ (seeds_[0] << a); 1113 seeds_[0] = seeds_[1]; 1114 seeds_[1] = seeds_[2]; 1115 seeds_[2] = seeds_[2] ^ (seeds_[2] >> c) ^ temp ^ (temp >> b); 1116 } 1117 else static if (bits == 128) 1118 { 1119 temp = seeds_[0] ^ (seeds_[0] << a); 1120 seeds_[0] = seeds_[1]; 1121 seeds_[1] = seeds_[2]; 1122 seeds_[2] = seeds_[3]; 1123 seeds_[3] = seeds_[3] ^ (seeds_[3] >> c) ^ temp ^ (temp >> b); 1124 } 1125 else static if (bits == 160) 1126 { 1127 temp = seeds_[0] ^ (seeds_[0] << a); 1128 seeds_[0] = seeds_[1]; 1129 seeds_[1] = seeds_[2]; 1130 seeds_[2] = seeds_[3]; 1131 seeds_[3] = seeds_[4]; 1132 seeds_[4] = seeds_[4] ^ (seeds_[4] >> c) ^ temp ^ (temp >> b); 1133 } 1134 else static if (bits == 192) 1135 { 1136 temp = seeds_[0] ^ (seeds_[0] >> a); 1137 seeds_[0] = seeds_[1]; 1138 seeds_[1] = seeds_[2]; 1139 seeds_[2] = seeds_[3]; 1140 seeds_[3] = seeds_[4]; 1141 seeds_[4] = seeds_[4] ^ (seeds_[4] << c) ^ temp ^ (temp << b); 1142 value_ = seeds_[4] + (seeds_[5] += 362_437); 1143 } 1144 else 1145 { 1146 static assert(false, "Phobos Error: Xorshift has no popFront() update for " 1147 ~ to!string(bits) ~ " bits."); 1148 } 1149 } 1150 1151 1152 /** 1153 * Captures a range state. 1154 */ 1155 @property 1156 typeof(this) save() @safe pure nothrow @nogc 1157 { 1158 return this; 1159 } 1160 1161 1162 /** 1163 * Compares against $(D_PARAM rhs) for equality. 1164 */ 1165 bool opEquals(ref const XorshiftEngine rhs) const @safe pure nothrow @nogc 1166 { 1167 return seeds_ == rhs.seeds_; 1168 } 1169 1170 1171 private: 1172 static void sanitizeSeeds(ref UIntType[size] seeds) @safe pure nothrow @nogc 1173 { 1174 for (uint i; i < seeds.length; i++) 1175 { 1176 if (seeds[i] == 0) 1177 seeds[i] = i + 1; 1178 } 1179 } 1180 1181 1182 @safe pure nothrow unittest 1183 { 1184 static if (size == 4) // Other bits too 1185 { 1186 UIntType[size] seeds = [1, 0, 0, 4]; 1187 1188 sanitizeSeeds(seeds); 1189 1190 assert(seeds == [1, 2, 3, 4]); 1191 } 1192 } 1193} 1194 1195 1196/** 1197 * Define $(D XorshiftEngine) generators with well-chosen parameters. See each bits examples of "Xorshift RNGs". 1198 * $(D Xorshift) is a Xorshift128's alias because 128bits implementation is mostly used. 1199 */ 1200alias Xorshift32 = XorshiftEngine!(uint, 32, 13, 17, 15) ; 1201alias Xorshift64 = XorshiftEngine!(uint, 64, 10, 13, 10); /// ditto 1202alias Xorshift96 = XorshiftEngine!(uint, 96, 10, 5, 26); /// ditto 1203alias Xorshift128 = XorshiftEngine!(uint, 128, 11, 8, 19); /// ditto 1204alias Xorshift160 = XorshiftEngine!(uint, 160, 2, 1, 4); /// ditto 1205alias Xorshift192 = XorshiftEngine!(uint, 192, 2, 1, 4); /// ditto 1206alias Xorshift = Xorshift128; /// ditto 1207 1208/// 1209@safe unittest 1210{ 1211 // Seed with a constant 1212 auto rnd = Xorshift(1); 1213 auto num = rnd.front; // same for each run 1214 1215 // Seed with an unpredictable value 1216 rnd.seed(unpredictableSeed); 1217 num = rnd.front; // different across rnd 1218} 1219 1220@safe unittest 1221{ 1222 import std.range; 1223 static assert(isForwardRange!Xorshift); 1224 static assert(isUniformRNG!Xorshift); 1225 static assert(isUniformRNG!(Xorshift, uint)); 1226 static assert(isSeedable!Xorshift); 1227 static assert(isSeedable!(Xorshift, uint)); 1228 1229 // Result from reference implementation. 1230 auto checking = [ 1231 [2463534242UL, 901999875, 3371835698, 2675058524, 1053936272, 3811264849, 1232 472493137, 3856898176, 2131710969, 2312157505], 1233 [362436069UL, 2113136921, 19051112, 3010520417, 951284840, 1213972223, 1234 3173832558, 2611145638, 2515869689, 2245824891], 1235 [521288629UL, 1950277231, 185954712, 1582725458, 3580567609, 2303633688, 1236 2394948066, 4108622809, 1116800180, 3357585673], 1237 [88675123UL, 3701687786, 458299110, 2500872618, 3633119408, 516391518, 1238 2377269574, 2599949379, 717229868, 137866584], 1239 [5783321UL, 393427209, 1947109840, 565829276, 1006220149, 971147905, 1240 1436324242, 2800460115, 1484058076, 3823330032], 1241 [0UL, 246875399, 3690007200, 1264581005, 3906711041, 1866187943, 2481925219, 1242 2464530826, 1604040631, 3653403911] 1243 ]; 1244 1245 alias XorshiftTypes = std.meta.AliasSeq!(Xorshift32, Xorshift64, Xorshift96, Xorshift128, Xorshift160, Xorshift192); 1246 1247 foreach (I, Type; XorshiftTypes) 1248 { 1249 Type rnd; 1250 1251 foreach (e; checking[I]) 1252 { 1253 assert(rnd.front == e); 1254 rnd.popFront(); 1255 } 1256 } 1257 1258 // Check .save works 1259 foreach (Type; XorshiftTypes) 1260 { 1261 auto rnd1 = Type(unpredictableSeed); 1262 auto rnd2 = rnd1.save; 1263 assert(rnd1 == rnd2); 1264 // Enable next test when RNGs are reference types 1265 version (none) { assert(rnd1 !is rnd2); } 1266 assert(rnd1.take(100).array() == rnd2.take(100).array()); 1267 } 1268} 1269 1270 1271/* A complete list of all pseudo-random number generators implemented in 1272 * std.random. This can be used to confirm that a given function or 1273 * object is compatible with all the pseudo-random number generators 1274 * available. It is enabled only in unittest mode. 1275 */ 1276@safe unittest 1277{ 1278 foreach (Rng; PseudoRngTypes) 1279 { 1280 static assert(isUniformRNG!Rng); 1281 auto rng = Rng(unpredictableSeed); 1282 } 1283} 1284 1285 1286/** 1287A "good" seed for initializing random number engines. Initializing 1288with $(D_PARAM unpredictableSeed) makes engines generate different 1289random number sequences every run. 1290 1291Returns: 1292A single unsigned integer seed value, different on each successive call 1293*/ 1294@property uint unpredictableSeed() @trusted 1295{ 1296 import core.thread : Thread, getpid, MonoTime; 1297 static bool seeded; 1298 static MinstdRand0 rand; 1299 if (!seeded) 1300 { 1301 uint threadID = cast(uint) cast(void*) Thread.getThis(); 1302 rand.seed((getpid() + threadID) ^ cast(uint) MonoTime.currTime.ticks); 1303 seeded = true; 1304 } 1305 rand.popFront(); 1306 return cast(uint) (MonoTime.currTime.ticks ^ rand.front); 1307} 1308 1309/// 1310@safe unittest 1311{ 1312 auto rnd = Random(unpredictableSeed); 1313 auto n = rnd.front; 1314 static assert(is(typeof(n) == uint)); 1315} 1316 1317/** 1318The "default", "favorite", "suggested" random number generator type on 1319the current platform. It is an alias for one of the previously-defined 1320generators. You may want to use it if (1) you need to generate some 1321nice random numbers, and (2) you don't care for the minutiae of the 1322method being used. 1323 */ 1324 1325alias Random = Mt19937; 1326 1327@safe unittest 1328{ 1329 static assert(isUniformRNG!Random); 1330 static assert(isUniformRNG!(Random, uint)); 1331 static assert(isSeedable!Random); 1332 static assert(isSeedable!(Random, uint)); 1333} 1334 1335/** 1336Global random number generator used by various functions in this 1337module whenever no generator is specified. It is allocated per-thread 1338and initialized to an unpredictable value for each thread. 1339 1340Returns: 1341A singleton instance of the default random number generator 1342 */ 1343@property ref Random rndGen() @safe 1344{ 1345 import std.algorithm.iteration : map; 1346 import std.range : repeat; 1347 1348 static Random result; 1349 static bool initialized; 1350 if (!initialized) 1351 { 1352 static if (isSeedable!(Random, typeof(map!((a) => unpredictableSeed)(repeat(0))))) 1353 result.seed(map!((a) => unpredictableSeed)(repeat(0))); 1354 else 1355 result = Random(unpredictableSeed); 1356 initialized = true; 1357 } 1358 return result; 1359} 1360 1361/** 1362Generates a number between $(D a) and $(D b). The $(D boundaries) 1363parameter controls the shape of the interval (open vs. closed on 1364either side). Valid values for $(D boundaries) are $(D "[]"), $(D 1365"$(LPAREN)]"), $(D "[$(RPAREN)"), and $(D "()"). The default interval 1366is closed to the left and open to the right. The version that does not 1367take $(D urng) uses the default generator $(D rndGen). 1368 1369Params: 1370 a = lower bound of the _uniform distribution 1371 b = upper bound of the _uniform distribution 1372 urng = (optional) random number generator to use; 1373 if not specified, defaults to $(D rndGen) 1374 1375Returns: 1376 A single random variate drawn from the _uniform distribution 1377 between $(D a) and $(D b), whose type is the common type of 1378 these parameters 1379 */ 1380auto uniform(string boundaries = "[)", T1, T2) 1381(T1 a, T2 b) 1382if (!is(CommonType!(T1, T2) == void)) 1383{ 1384 return uniform!(boundaries, T1, T2, Random)(a, b, rndGen); 1385} 1386 1387/// 1388@safe unittest 1389{ 1390 auto gen = Random(unpredictableSeed); 1391 // Generate an integer in [0, 1023] 1392 auto a = uniform(0, 1024, gen); 1393 // Generate a float in [0, 1) 1394 auto b = uniform(0.0f, 1.0f, gen); 1395} 1396 1397@safe unittest 1398{ 1399 MinstdRand0 gen; 1400 foreach (i; 0 .. 20) 1401 { 1402 auto x = uniform(0.0, 15.0, gen); 1403 assert(0 <= x && x < 15); 1404 } 1405 foreach (i; 0 .. 20) 1406 { 1407 auto x = uniform!"[]"('a', 'z', gen); 1408 assert('a' <= x && x <= 'z'); 1409 } 1410 1411 foreach (i; 0 .. 20) 1412 { 1413 auto x = uniform('a', 'z', gen); 1414 assert('a' <= x && x < 'z'); 1415 } 1416 1417 foreach (i; 0 .. 20) 1418 { 1419 immutable ubyte a = 0; 1420 immutable ubyte b = 15; 1421 auto x = uniform(a, b, gen); 1422 assert(a <= x && x < b); 1423 } 1424} 1425 1426// Implementation of uniform for floating-point types 1427/// ditto 1428auto uniform(string boundaries = "[)", 1429 T1, T2, UniformRandomNumberGenerator) 1430(T1 a, T2 b, ref UniformRandomNumberGenerator urng) 1431if (isFloatingPoint!(CommonType!(T1, T2)) && isUniformRNG!UniformRandomNumberGenerator) 1432{ 1433 import std.conv : text; 1434 import std.exception : enforce; 1435 alias NumberType = Unqual!(CommonType!(T1, T2)); 1436 static if (boundaries[0] == '(') 1437 { 1438 import std.math : nextafter; 1439 NumberType _a = nextafter(cast(NumberType) a, NumberType.infinity); 1440 } 1441 else 1442 { 1443 NumberType _a = a; 1444 } 1445 static if (boundaries[1] == ')') 1446 { 1447 import std.math : nextafter; 1448 NumberType _b = nextafter(cast(NumberType) b, -NumberType.infinity); 1449 } 1450 else 1451 { 1452 NumberType _b = b; 1453 } 1454 enforce(_a <= _b, 1455 text("std.random.uniform(): invalid bounding interval ", 1456 boundaries[0], a, ", ", b, boundaries[1])); 1457 NumberType result = 1458 _a + (_b - _a) * cast(NumberType) (urng.front - urng.min) 1459 / (urng.max - urng.min); 1460 urng.popFront(); 1461 return result; 1462} 1463 1464// Implementation of uniform for integral types 1465/+ Description of algorithm and suggestion of correctness: 1466 1467The modulus operator maps an integer to a small, finite space. For instance, `x 1468% 3` will map whatever x is into the range [0 .. 3). 0 maps to 0, 1 maps to 1, 2 1469maps to 2, 3 maps to 0, and so on infinitely. As long as the integer is 1470uniformly chosen from the infinite space of all non-negative integers then `x % 14713` will uniformly fall into that range. 1472 1473(Non-negative is important in this case because some definitions of modulus, 1474namely the one used in computers generally, map negative numbers differently to 1475(-3 .. 0]. `uniform` does not use negative number modulus, thus we can safely 1476ignore that fact.) 1477 1478The issue with computers is that integers have a finite space they must fit in, 1479and our uniformly chosen random number is picked in that finite space. So, that 1480method is not sufficient. You can look at it as the integer space being divided 1481into "buckets" and every bucket after the first bucket maps directly into that 1482first bucket. `[0, 1, 2]`, `[3, 4, 5]`, ... When integers are finite, then the 1483last bucket has the chance to be "incomplete": `[uint.max - 3, uint.max - 2, 1484uint.max - 1]`, `[uint.max]` ... (the last bucket only has 1!). The issue here 1485is that _every_ bucket maps _completely_ to the first bucket except for that 1486last one. The last one doesn't have corresponding mappings to 1 or 2, in this 1487case, which makes it unfair. 1488 1489So, the answer is to simply "reroll" if you're in that last bucket, since it's 1490the only unfair one. Eventually you'll roll into a fair bucket. Simply, instead 1491of the meaning of the last bucket being "maps to `[0]`", it changes to "maps to 1492`[0, 1, 2]`", which is precisely what we want. 1493 1494To generalize, `upperDist` represents the size of our buckets (and, thus, the 1495exclusive upper bound for our desired uniform number). `rnum` is a uniformly 1496random number picked from the space of integers that a computer can hold (we'll 1497say `UpperType` represents that type). 1498 1499We'll first try to do the mapping into the first bucket by doing `offset = rnum 1500% upperDist`. We can figure out the position of the front of the bucket we're in 1501by `bucketFront = rnum - offset`. 1502 1503If we start at `UpperType.max` and walk backwards `upperDist - 1` spaces, then 1504the space we land on is the last acceptable position where a full bucket can 1505fit: 1506 1507``` 1508 bucketFront UpperType.max 1509 v v 1510[..., 0, 1, 2, ..., upperDist - 1] 1511 ^~~ upperDist - 1 ~~^ 1512``` 1513 1514If the bucket starts any later, then it must have lost at least one number and 1515at least that number won't be represented fairly. 1516 1517``` 1518 bucketFront UpperType.max 1519 v v 1520[..., upperDist - 1, 0, 1, 2, ..., upperDist - 2] 1521 ^~~~~~~~ upperDist - 1 ~~~~~~~^ 1522``` 1523 1524Hence, our condition to reroll is 1525`bucketFront > (UpperType.max - (upperDist - 1))` 1526+/ 1527auto uniform(string boundaries = "[)", T1, T2, RandomGen) 1528(T1 a, T2 b, ref RandomGen rng) 1529if ((isIntegral!(CommonType!(T1, T2)) || isSomeChar!(CommonType!(T1, T2))) && 1530 isUniformRNG!RandomGen) 1531{ 1532 import std.conv : text, unsigned; 1533 import std.exception : enforce; 1534 alias ResultType = Unqual!(CommonType!(T1, T2)); 1535 static if (boundaries[0] == '(') 1536 { 1537 enforce(a < ResultType.max, 1538 text("std.random.uniform(): invalid left bound ", a)); 1539 ResultType lower = cast(ResultType) (a + 1); 1540 } 1541 else 1542 { 1543 ResultType lower = a; 1544 } 1545 1546 static if (boundaries[1] == ']') 1547 { 1548 enforce(lower <= b, 1549 text("std.random.uniform(): invalid bounding interval ", 1550 boundaries[0], a, ", ", b, boundaries[1])); 1551 /* Cannot use this next optimization with dchar, as dchar 1552 * only partially uses its full bit range 1553 */ 1554 static if (!is(ResultType == dchar)) 1555 { 1556 if (b == ResultType.max && lower == ResultType.min) 1557 { 1558 // Special case - all bits are occupied 1559 return std.random.uniform!ResultType(rng); 1560 } 1561 } 1562 auto upperDist = unsigned(b - lower) + 1u; 1563 } 1564 else 1565 { 1566 enforce(lower < b, 1567 text("std.random.uniform(): invalid bounding interval ", 1568 boundaries[0], a, ", ", b, boundaries[1])); 1569 auto upperDist = unsigned(b - lower); 1570 } 1571 1572 assert(upperDist != 0); 1573 1574 alias UpperType = typeof(upperDist); 1575 static assert(UpperType.min == 0); 1576 1577 UpperType offset, rnum, bucketFront; 1578 do 1579 { 1580 rnum = uniform!UpperType(rng); 1581 offset = rnum % upperDist; 1582 bucketFront = rnum - offset; 1583 } // while we're in an unfair bucket... 1584 while (bucketFront > (UpperType.max - (upperDist - 1))); 1585 1586 return cast(ResultType)(lower + offset); 1587} 1588 1589@safe unittest 1590{ 1591 import std.conv : to; 1592 auto gen = Mt19937(unpredictableSeed); 1593 static assert(isForwardRange!(typeof(gen))); 1594 1595 auto a = uniform(0, 1024, gen); 1596 assert(0 <= a && a <= 1024); 1597 auto b = uniform(0.0f, 1.0f, gen); 1598 assert(0 <= b && b < 1, to!string(b)); 1599 auto c = uniform(0.0, 1.0); 1600 assert(0 <= c && c < 1); 1601 1602 foreach (T; std.meta.AliasSeq!(char, wchar, dchar, byte, ubyte, short, ushort, 1603 int, uint, long, ulong, float, double, real)) 1604 { 1605 T lo = 0, hi = 100; 1606 1607 // Try tests with each of the possible bounds 1608 { 1609 T init = uniform(lo, hi); 1610 size_t i = 50; 1611 while (--i && uniform(lo, hi) == init) {} 1612 assert(i > 0); 1613 } 1614 { 1615 T init = uniform!"[)"(lo, hi); 1616 size_t i = 50; 1617 while (--i && uniform(lo, hi) == init) {} 1618 assert(i > 0); 1619 } 1620 { 1621 T init = uniform!"(]"(lo, hi); 1622 size_t i = 50; 1623 while (--i && uniform(lo, hi) == init) {} 1624 assert(i > 0); 1625 } 1626 { 1627 T init = uniform!"()"(lo, hi); 1628 size_t i = 50; 1629 while (--i && uniform(lo, hi) == init) {} 1630 assert(i > 0); 1631 } 1632 { 1633 T init = uniform!"[]"(lo, hi); 1634 size_t i = 50; 1635 while (--i && uniform(lo, hi) == init) {} 1636 assert(i > 0); 1637 } 1638 1639 /* Test case with closed boundaries covering whole range 1640 * of integral type 1641 */ 1642 static if (isIntegral!T || isSomeChar!T) 1643 { 1644 foreach (immutable _; 0 .. 100) 1645 { 1646 auto u = uniform!"[]"(T.min, T.max); 1647 static assert(is(typeof(u) == T)); 1648 assert(T.min <= u, "Lower bound violation for uniform!\"[]\" with " ~ T.stringof); 1649 assert(u <= T.max, "Upper bound violation for uniform!\"[]\" with " ~ T.stringof); 1650 } 1651 } 1652 } 1653 1654 auto reproRng = Xorshift(239842); 1655 1656 foreach (T; std.meta.AliasSeq!(char, wchar, dchar, byte, ubyte, short, 1657 ushort, int, uint, long, ulong)) 1658 { 1659 T lo = T.min + 10, hi = T.max - 10; 1660 T init = uniform(lo, hi, reproRng); 1661 size_t i = 50; 1662 while (--i && uniform(lo, hi, reproRng) == init) {} 1663 assert(i > 0); 1664 } 1665 1666 { 1667 bool sawLB = false, sawUB = false; 1668 foreach (i; 0 .. 50) 1669 { 1670 auto x = uniform!"[]"('a', 'd', reproRng); 1671 if (x == 'a') sawLB = true; 1672 if (x == 'd') sawUB = true; 1673 assert('a' <= x && x <= 'd'); 1674 } 1675 assert(sawLB && sawUB); 1676 } 1677 1678 { 1679 bool sawLB = false, sawUB = false; 1680 foreach (i; 0 .. 50) 1681 { 1682 auto x = uniform('a', 'd', reproRng); 1683 if (x == 'a') sawLB = true; 1684 if (x == 'c') sawUB = true; 1685 assert('a' <= x && x < 'd'); 1686 } 1687 assert(sawLB && sawUB); 1688 } 1689 1690 { 1691 bool sawLB = false, sawUB = false; 1692 foreach (i; 0 .. 50) 1693 { 1694 immutable int lo = -2, hi = 2; 1695 auto x = uniform!"()"(lo, hi, reproRng); 1696 if (x == (lo+1)) sawLB = true; 1697 if (x == (hi-1)) sawUB = true; 1698 assert(lo < x && x < hi); 1699 } 1700 assert(sawLB && sawUB); 1701 } 1702 1703 { 1704 bool sawLB = false, sawUB = false; 1705 foreach (i; 0 .. 50) 1706 { 1707 immutable ubyte lo = 0, hi = 5; 1708 auto x = uniform(lo, hi, reproRng); 1709 if (x == lo) sawLB = true; 1710 if (x == (hi-1)) sawUB = true; 1711 assert(lo <= x && x < hi); 1712 } 1713 assert(sawLB && sawUB); 1714 } 1715 1716 { 1717 foreach (i; 0 .. 30) 1718 { 1719 assert(i == uniform(i, i+1, reproRng)); 1720 } 1721 } 1722} 1723 1724/** 1725Generates a uniformly-distributed number in the range $(D [T.min, 1726T.max]) for any integral or character type $(D T). If no random 1727number generator is passed, uses the default $(D rndGen). 1728 1729Params: 1730 urng = (optional) random number generator to use; 1731 if not specified, defaults to $(D rndGen) 1732 1733Returns: 1734 Random variate drawn from the _uniform distribution across all 1735 possible values of the integral or character type $(D T). 1736 */ 1737auto uniform(T, UniformRandomNumberGenerator) 1738(ref UniformRandomNumberGenerator urng) 1739if (!is(T == enum) && (isIntegral!T || isSomeChar!T) && isUniformRNG!UniformRandomNumberGenerator) 1740{ 1741 /* dchar does not use its full bit range, so we must 1742 * revert to the uniform with specified bounds 1743 */ 1744 static if (is(T == dchar)) 1745 { 1746 return uniform!"[]"(T.min, T.max); 1747 } 1748 else 1749 { 1750 auto r = urng.front; 1751 urng.popFront(); 1752 static if (T.sizeof <= r.sizeof) 1753 { 1754 return cast(T) r; 1755 } 1756 else 1757 { 1758 static assert(T.sizeof == 8 && r.sizeof == 4); 1759 T r1 = urng.front | (cast(T) r << 32); 1760 urng.popFront(); 1761 return r1; 1762 } 1763 } 1764} 1765 1766/// Ditto 1767auto uniform(T)() 1768if (!is(T == enum) && (isIntegral!T || isSomeChar!T)) 1769{ 1770 return uniform!T(rndGen); 1771} 1772 1773@safe unittest 1774{ 1775 foreach (T; std.meta.AliasSeq!(char, wchar, dchar, byte, ubyte, short, ushort, 1776 int, uint, long, ulong)) 1777 { 1778 T init = uniform!T(); 1779 size_t i = 50; 1780 while (--i && uniform!T() == init) {} 1781 assert(i > 0); 1782 1783 foreach (immutable _; 0 .. 100) 1784 { 1785 auto u = uniform!T(); 1786 static assert(is(typeof(u) == T)); 1787 assert(T.min <= u, "Lower bound violation for uniform!" ~ T.stringof); 1788 assert(u <= T.max, "Upper bound violation for uniform!" ~ T.stringof); 1789 } 1790 } 1791} 1792 1793/** 1794Returns a uniformly selected member of enum $(D E). If no random number 1795generator is passed, uses the default $(D rndGen). 1796 1797Params: 1798 urng = (optional) random number generator to use; 1799 if not specified, defaults to $(D rndGen) 1800 1801Returns: 1802 Random variate drawn with equal probability from any 1803 of the possible values of the enum $(D E). 1804 */ 1805auto uniform(E, UniformRandomNumberGenerator) 1806(ref UniformRandomNumberGenerator urng) 1807if (is(E == enum) && isUniformRNG!UniformRandomNumberGenerator) 1808{ 1809 static immutable E[EnumMembers!E.length] members = [EnumMembers!E]; 1810 return members[std.random.uniform(0, members.length, urng)]; 1811} 1812 1813/// Ditto 1814auto uniform(E)() 1815if (is(E == enum)) 1816{ 1817 return uniform!E(rndGen); 1818} 1819 1820/// 1821@safe unittest 1822{ 1823 enum Fruit { apple, mango, pear } 1824 auto randFruit = uniform!Fruit(); 1825} 1826 1827@safe unittest 1828{ 1829 enum Fruit { Apple = 12, Mango = 29, Pear = 72 } 1830 foreach (_; 0 .. 100) 1831 { 1832 foreach (f; [uniform!Fruit(), rndGen.uniform!Fruit()]) 1833 { 1834 assert(f == Fruit.Apple || f == Fruit.Mango || f == Fruit.Pear); 1835 } 1836 } 1837} 1838 1839/** 1840 * Generates a uniformly-distributed floating point number of type 1841 * $(D T) in the range [0, 1$(RPAREN). If no random number generator is 1842 * specified, the default RNG $(D rndGen) will be used as the source 1843 * of randomness. 1844 * 1845 * $(D uniform01) offers a faster generation of random variates than 1846 * the equivalent $(D uniform!"[$(RPAREN)"(0.0, 1.0)) and so may be preferred 1847 * for some applications. 1848 * 1849 * Params: 1850 * rng = (optional) random number generator to use; 1851 * if not specified, defaults to $(D rndGen) 1852 * 1853 * Returns: 1854 * Floating-point random variate of type $(D T) drawn from the _uniform 1855 * distribution across the half-open interval [0, 1$(RPAREN). 1856 * 1857 */ 1858T uniform01(T = double)() 1859if (isFloatingPoint!T) 1860{ 1861 return uniform01!T(rndGen); 1862} 1863 1864/// ditto 1865T uniform01(T = double, UniformRNG)(ref UniformRNG rng) 1866if (isFloatingPoint!T && isUniformRNG!UniformRNG) 1867out (result) 1868{ 1869 assert(0 <= result); 1870 assert(result < 1); 1871} 1872body 1873{ 1874 alias R = typeof(rng.front); 1875 static if (isIntegral!R) 1876 { 1877 enum T factor = 1 / (T(1) + rng.max - rng.min); 1878 } 1879 else static if (isFloatingPoint!R) 1880 { 1881 enum T factor = 1 / (rng.max - rng.min); 1882 } 1883 else 1884 { 1885 static assert(false); 1886 } 1887 1888 while (true) 1889 { 1890 immutable T u = (rng.front - rng.min) * factor; 1891 rng.popFront(); 1892 1893 import core.stdc.limits : CHAR_BIT; // CHAR_BIT is always 8 1894 static if (isIntegral!R && T.mant_dig >= (CHAR_BIT * R.sizeof)) 1895 { 1896 /* If RNG variates are integral and T has enough precision to hold 1897 * R without loss, we're guaranteed by the definition of factor 1898 * that precisely u < 1. 1899 */ 1900 return u; 1901 } 1902 else 1903 { 1904 /* Otherwise we have to check whether u is beyond the assumed range 1905 * because of the loss of precision, or for another reason, a 1906 * floating-point RNG can return a variate that is exactly equal to 1907 * its maximum. 1908 */ 1909 if (u < 1) 1910 { 1911 return u; 1912 } 1913 } 1914 } 1915 1916 // Shouldn't ever get here. 1917 assert(false); 1918} 1919 1920@safe unittest 1921{ 1922 import std.meta; 1923 foreach (UniformRNG; PseudoRngTypes) 1924 { 1925 1926 foreach (T; std.meta.AliasSeq!(float, double, real)) 1927 (){ // avoid slow optimizations for large functions @@@BUG@@@ 2396 1928 UniformRNG rng = UniformRNG(unpredictableSeed); 1929 1930 auto a = uniform01(); 1931 assert(is(typeof(a) == double)); 1932 assert(0 <= a && a < 1); 1933 1934 auto b = uniform01(rng); 1935 assert(is(typeof(a) == double)); 1936 assert(0 <= b && b < 1); 1937 1938 auto c = uniform01!T(); 1939 assert(is(typeof(c) == T)); 1940 assert(0 <= c && c < 1); 1941 1942 auto d = uniform01!T(rng); 1943 assert(is(typeof(d) == T)); 1944 assert(0 <= d && d < 1); 1945 1946 T init = uniform01!T(rng); 1947 size_t i = 50; 1948 while (--i && uniform01!T(rng) == init) {} 1949 assert(i > 0); 1950 assert(i < 50); 1951 }(); 1952 } 1953} 1954 1955/** 1956Generates a uniform probability distribution of size $(D n), i.e., an 1957array of size $(D n) of positive numbers of type $(D F) that sum to 1958$(D 1). If $(D useThis) is provided, it is used as storage. 1959 */ 1960F[] uniformDistribution(F = double)(size_t n, F[] useThis = null) 1961if (isFloatingPoint!F) 1962{ 1963 import std.numeric : normalize; 1964 useThis.length = n; 1965 foreach (ref e; useThis) 1966 { 1967 e = uniform(0.0, 1); 1968 } 1969 normalize(useThis); 1970 return useThis; 1971} 1972 1973@safe unittest 1974{ 1975 import std.algorithm; 1976 import std.math; 1977 static assert(is(CommonType!(double, int) == double)); 1978 auto a = uniformDistribution(5); 1979 assert(a.length == 5); 1980 assert(approxEqual(reduce!"a + b"(a), 1)); 1981 a = uniformDistribution(10, a); 1982 assert(a.length == 10); 1983 assert(approxEqual(reduce!"a + b"(a), 1)); 1984} 1985 1986/** 1987Returns a random, uniformly chosen, element `e` from the supplied 1988$(D Range range). If no random number generator is passed, the default 1989`rndGen` is used. 1990 1991Params: 1992 range = a random access range that has the `length` property defined 1993 urng = (optional) random number generator to use; 1994 if not specified, defaults to `rndGen` 1995 1996Returns: 1997 A single random element drawn from the `range`. If it can, it will 1998 return a `ref` to the $(D range element), otherwise it will return 1999 a copy. 2000 */ 2001auto ref choice(Range, RandomGen = Random)(auto ref Range range, 2002 ref RandomGen urng = rndGen) 2003if (isRandomAccessRange!Range && hasLength!Range && isUniformRNG!RandomGen) 2004{ 2005 assert(range.length > 0, 2006 __PRETTY_FUNCTION__ ~ ": invalid Range supplied. Range cannot be empty"); 2007 2008 return range[uniform(size_t(0), $, urng)]; 2009} 2010 2011/// 2012@safe unittest 2013{ 2014 import std.algorithm.searching : canFind; 2015 2016 auto array = [1, 2, 3, 4, 5]; 2017 auto elem = choice(array); 2018 2019 assert(canFind(array, elem), 2020 "Choice did not return a valid element from the given Range"); 2021 2022 auto urng = Random(unpredictableSeed); 2023 elem = choice(array, urng); 2024 2025 assert(canFind(array, elem), 2026 "Choice did not return a valid element from the given Range"); 2027} 2028 2029@safe unittest 2030{ 2031 import std.algorithm.searching : canFind; 2032 2033 class MyTestClass 2034 { 2035 int x; 2036 2037 this(int x) 2038 { 2039 this.x = x; 2040 } 2041 } 2042 2043 MyTestClass[] testClass; 2044 foreach (i; 0 .. 5) 2045 { 2046 testClass ~= new MyTestClass(i); 2047 } 2048 2049 auto elem = choice(testClass); 2050 2051 assert(canFind!((ref MyTestClass a, ref MyTestClass b) => a.x == b.x)(testClass, elem), 2052 "Choice did not return a valid element from the given Range"); 2053} 2054 2055@system unittest 2056{ 2057 import std.algorithm.iteration : map; 2058 import std.algorithm.searching : canFind; 2059 2060 auto array = [1, 2, 3, 4, 5]; 2061 auto elemAddr = &choice(array); 2062 2063 assert(array.map!((ref e) => &e).canFind(elemAddr), 2064 "Choice did not return a ref to an element from the given Range"); 2065 assert(array.canFind(*(cast(int *)(elemAddr))), 2066 "Choice did not return a valid element from the given Range"); 2067} 2068 2069/** 2070Shuffles elements of $(D r) using $(D gen) as a shuffler. $(D r) must be 2071a random-access range with length. If no RNG is specified, $(D rndGen) 2072will be used. 2073 2074Params: 2075 r = random-access range whose elements are to be shuffled 2076 gen = (optional) random number generator to use; if not 2077 specified, defaults to $(D rndGen) 2078 */ 2079 2080void randomShuffle(Range, RandomGen)(Range r, ref RandomGen gen) 2081if (isRandomAccessRange!Range && isUniformRNG!RandomGen) 2082{ 2083 return partialShuffle!(Range, RandomGen)(r, r.length, gen); 2084} 2085 2086/// ditto 2087void randomShuffle(Range)(Range r) 2088if (isRandomAccessRange!Range) 2089{ 2090 return randomShuffle(r, rndGen); 2091} 2092 2093@safe unittest 2094{ 2095 import std.algorithm.sorting : sort; 2096 foreach (RandomGen; PseudoRngTypes) 2097 { 2098 // Also tests partialShuffle indirectly. 2099 auto a = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]; 2100 auto b = a.dup; 2101 auto gen = RandomGen(unpredictableSeed); 2102 randomShuffle(a, gen); 2103 sort(a); 2104 assert(a == b); 2105 randomShuffle(a); 2106 sort(a); 2107 assert(a == b); 2108 } 2109} 2110 2111/** 2112Partially shuffles the elements of $(D r) such that upon returning $(D r[0 .. n]) 2113is a random subset of $(D r) and is randomly ordered. $(D r[n .. r.length]) 2114will contain the elements not in $(D r[0 .. n]). These will be in an undefined 2115order, but will not be random in the sense that their order after 2116$(D partialShuffle) returns will not be independent of their order before 2117$(D partialShuffle) was called. 2118 2119$(D r) must be a random-access range with length. $(D n) must be less than 2120or equal to $(D r.length). If no RNG is specified, $(D rndGen) will be used. 2121 2122Params: 2123 r = random-access range whose elements are to be shuffled 2124 n = number of elements of $(D r) to shuffle (counting from the beginning); 2125 must be less than $(D r.length) 2126 gen = (optional) random number generator to use; if not 2127 specified, defaults to $(D rndGen) 2128*/ 2129void partialShuffle(Range, RandomGen)(Range r, in size_t n, ref RandomGen gen) 2130if (isRandomAccessRange!Range && isUniformRNG!RandomGen) 2131{ 2132 import std.algorithm.mutation : swapAt; 2133 import std.exception : enforce; 2134 enforce(n <= r.length, "n must be <= r.length for partialShuffle."); 2135 foreach (i; 0 .. n) 2136 { 2137 r.swapAt(i, uniform(i, r.length, gen)); 2138 } 2139} 2140 2141/// ditto 2142void partialShuffle(Range)(Range r, in size_t n) 2143if (isRandomAccessRange!Range) 2144{ 2145 return partialShuffle(r, n, rndGen); 2146} 2147 2148@safe unittest 2149{ 2150 import std.algorithm; 2151 foreach (RandomGen; PseudoRngTypes) 2152 { 2153 auto a = [0, 1, 1, 2, 3]; 2154 auto b = a.dup; 2155 2156 // Pick a fixed seed so that the outcome of the statistical 2157 // test below is deterministic. 2158 auto gen = RandomGen(12345); 2159 2160 // NUM times, pick LEN elements from the array at random. 2161 immutable int LEN = 2; 2162 immutable int NUM = 750; 2163 int[][] chk; 2164 foreach (step; 0 .. NUM) 2165 { 2166 partialShuffle(a, LEN, gen); 2167 chk ~= a[0 .. LEN].dup; 2168 } 2169 2170 // Check that each possible a[0 .. LEN] was produced at least once. 2171 // For a perfectly random RandomGen, the probability that each 2172 // particular combination failed to appear would be at most 2173 // 0.95 ^^ NUM which is approximately 1,962e-17. 2174 // As long as hardware failure (e.g. bit flip) probability 2175 // is higher, we are fine with this unittest. 2176 sort(chk); 2177 assert(equal(uniq(chk), [ [0,1], [0,2], [0,3], 2178 [1,0], [1,1], [1,2], [1,3], 2179 [2,0], [2,1], [2,3], 2180 [3,0], [3,1], [3,2], ])); 2181 2182 // Check that all the elements are still there. 2183 sort(a); 2184 assert(equal(a, b)); 2185 } 2186} 2187 2188/** 2189Rolls a dice with relative probabilities stored in $(D 2190proportions). Returns the index in $(D proportions) that was chosen. 2191 2192Params: 2193 rnd = (optional) random number generator to use; if not 2194 specified, defaults to $(D rndGen) 2195 proportions = forward range or list of individual values 2196 whose elements correspond to the probabilities 2197 with which to choose the corresponding index 2198 value 2199 2200Returns: 2201 Random variate drawn from the index values 2202 [0, ... $(D proportions.length) - 1], with the probability 2203 of getting an individual index value $(D i) being proportional to 2204 $(D proportions[i]). 2205*/ 2206size_t dice(Rng, Num)(ref Rng rnd, Num[] proportions...) 2207if (isNumeric!Num && isForwardRange!Rng) 2208{ 2209 return diceImpl(rnd, proportions); 2210} 2211 2212/// Ditto 2213size_t dice(R, Range)(ref R rnd, Range proportions) 2214if (isForwardRange!Range && isNumeric!(ElementType!Range) && !isArray!Range) 2215{ 2216 return diceImpl(rnd, proportions); 2217} 2218 2219/// Ditto 2220size_t dice(Range)(Range proportions) 2221if (isForwardRange!Range && isNumeric!(ElementType!Range) && !isArray!Range) 2222{ 2223 return diceImpl(rndGen, proportions); 2224} 2225 2226/// Ditto 2227size_t dice(Num)(Num[] proportions...) 2228if (isNumeric!Num) 2229{ 2230 return diceImpl(rndGen, proportions); 2231} 2232 2233/// 2234@safe unittest 2235{ 2236 auto x = dice(0.5, 0.5); // x is 0 or 1 in equal proportions 2237 auto y = dice(50, 50); // y is 0 or 1 in equal proportions 2238 auto z = dice(70, 20, 10); // z is 0 70% of the time, 1 20% of the time, 2239 // and 2 10% of the time 2240} 2241 2242private size_t diceImpl(Rng, Range)(ref Rng rng, scope Range proportions) 2243if (isForwardRange!Range && isNumeric!(ElementType!Range) && isForwardRange!Rng) 2244in 2245{ 2246 import std.algorithm.searching : all; 2247 assert(proportions.save.all!"a >= 0"); 2248} 2249body 2250{ 2251 import std.algorithm.iteration : reduce; 2252 import std.exception : enforce; 2253 double sum = reduce!"a + b"(0.0, proportions.save); 2254 enforce(sum > 0, "Proportions in a dice cannot sum to zero"); 2255 immutable point = uniform(0.0, sum, rng); 2256 assert(point < sum); 2257 auto mass = 0.0; 2258 2259 size_t i = 0; 2260 foreach (e; proportions) 2261 { 2262 mass += e; 2263 if (point < mass) return i; 2264 i++; 2265 } 2266 // this point should not be reached 2267 assert(false); 2268} 2269 2270@safe unittest 2271{ 2272 auto rnd = Random(unpredictableSeed); 2273 auto i = dice(rnd, 0.0, 100.0); 2274 assert(i == 1); 2275 i = dice(rnd, 100.0, 0.0); 2276 assert(i == 0); 2277 2278 i = dice(100U, 0U); 2279 assert(i == 0); 2280} 2281 2282/** 2283Covers a given range $(D r) in a random manner, i.e. goes through each 2284element of $(D r) once and only once, just in a random order. $(D r) 2285must be a random-access range with length. 2286 2287If no random number generator is passed to $(D randomCover), the 2288thread-global RNG rndGen will be used internally. 2289 2290Params: 2291 r = random-access range to cover 2292 rng = (optional) random number generator to use; 2293 if not specified, defaults to $(D rndGen) 2294 2295Returns: 2296 Range whose elements consist of the elements of $(D r), 2297 in random order. Will be a forward range if both $(D r) and 2298 $(D rng) are forward ranges, an input range otherwise. 2299 2300Example: 2301---- 2302int[] a = [ 0, 1, 2, 3, 4, 5, 6, 7, 8 ]; 2303foreach (e; randomCover(a)) 2304{ 2305 writeln(e); 2306} 2307---- 2308 2309$(B WARNING:) If an alternative RNG is desired, it is essential for this 2310to be a $(I new) RNG seeded in an unpredictable manner. Passing it a RNG 2311used elsewhere in the program will result in unintended correlations, 2312due to the current implementation of RNGs as value types. 2313 2314Example: 2315---- 2316int[] a = [ 0, 1, 2, 3, 4, 5, 6, 7, 8 ]; 2317foreach (e; randomCover(a, Random(unpredictableSeed))) // correct! 2318{ 2319 writeln(e); 2320} 2321 2322foreach (e; randomCover(a, rndGen)) // DANGEROUS!! rndGen gets copied by value 2323{ 2324 writeln(e); 2325} 2326 2327foreach (e; randomCover(a, rndGen)) // ... so this second random cover 2328{ // will output the same sequence as 2329 writeln(e); // the previous one. 2330} 2331---- 2332 */ 2333struct RandomCover(Range, UniformRNG = void) 2334if (isRandomAccessRange!Range && (isUniformRNG!UniformRNG || is(UniformRNG == void))) 2335{ 2336 private Range _input; 2337 private bool[] _chosen; 2338 private size_t _current; 2339 private size_t _alreadyChosen = 0; 2340 private bool _isEmpty = false; 2341 2342 static if (is(UniformRNG == void)) 2343 { 2344 this(Range input) 2345 { 2346 _input = input; 2347 _chosen.length = _input.length; 2348 if (_input.empty) 2349 { 2350 _isEmpty = true; 2351 } 2352 else 2353 { 2354 _current = uniform(0, _chosen.length); 2355 } 2356 } 2357 } 2358 else 2359 { 2360 private UniformRNG _rng; 2361 2362 this(Range input, ref UniformRNG rng) 2363 { 2364 _input = input; 2365 _rng = rng; 2366 _chosen.length = _input.length; 2367 if (_input.empty) 2368 { 2369 _isEmpty = true; 2370 } 2371 else 2372 { 2373 _current = uniform(0, _chosen.length, rng); 2374 } 2375 } 2376 2377 this(Range input, UniformRNG rng) 2378 { 2379 this(input, rng); 2380 } 2381 } 2382 2383 static if (hasLength!Range) 2384 { 2385 @property size_t length() 2386 { 2387 return _input.length - _alreadyChosen; 2388 } 2389 } 2390 2391 @property auto ref front() 2392 { 2393 assert(!_isEmpty); 2394 return _input[_current]; 2395 } 2396 2397 void popFront() 2398 { 2399 assert(!_isEmpty); 2400 2401 size_t k = _input.length - _alreadyChosen - 1; 2402 if (k == 0) 2403 { 2404 _isEmpty = true; 2405 ++_alreadyChosen; 2406 return; 2407 } 2408 2409 size_t i; 2410 foreach (e; _input) 2411 { 2412 if (_chosen[i] || i == _current) { ++i; continue; } 2413 // Roll a dice with k faces 2414 static if (is(UniformRNG == void)) 2415 { 2416 auto chooseMe = uniform(0, k) == 0; 2417 } 2418 else 2419 { 2420 auto chooseMe = uniform(0, k, _rng) == 0; 2421 } 2422 assert(k > 1 || chooseMe); 2423 if (chooseMe) 2424 { 2425 _chosen[_current] = true; 2426 _current = i; 2427 ++_alreadyChosen; 2428 return; 2429 } 2430 --k; 2431 ++i; 2432 } 2433 } 2434 2435 static if (isForwardRange!UniformRNG) 2436 { 2437 @property typeof(this) save() 2438 { 2439 auto ret = this; 2440 ret._input = _input.save; 2441 ret._rng = _rng.save; 2442 return ret; 2443 } 2444 } 2445 2446 @property bool empty() { return _isEmpty; } 2447} 2448 2449/// Ditto 2450auto randomCover(Range, UniformRNG)(Range r, auto ref UniformRNG rng) 2451if (isRandomAccessRange!Range && isUniformRNG!UniformRNG) 2452{ 2453 return RandomCover!(Range, UniformRNG)(r, rng); 2454} 2455 2456/// Ditto 2457auto randomCover(Range)(Range r) 2458if (isRandomAccessRange!Range) 2459{ 2460 return RandomCover!(Range, void)(r); 2461} 2462 2463@safe unittest 2464{ 2465 import std.algorithm; 2466 import std.conv; 2467 int[] a = [ 0, 1, 2, 3, 4, 5, 6, 7, 8 ]; 2468 int[] c; 2469 foreach (UniformRNG; std.meta.AliasSeq!(void, PseudoRngTypes)) 2470 { 2471 static if (is(UniformRNG == void)) 2472 { 2473 auto rc = randomCover(a); 2474 static assert(isInputRange!(typeof(rc))); 2475 static assert(!isForwardRange!(typeof(rc))); 2476 } 2477 else 2478 { 2479 auto rng = UniformRNG(unpredictableSeed); 2480 auto rc = randomCover(a, rng); 2481 static assert(isForwardRange!(typeof(rc))); 2482 // check for constructor passed a value-type RNG 2483 auto rc2 = RandomCover!(int[], UniformRNG)(a, UniformRNG(unpredictableSeed)); 2484 static assert(isForwardRange!(typeof(rc2))); 2485 auto rcEmpty = randomCover(c, rng); 2486 assert(rcEmpty.length == 0); 2487 } 2488 2489 int[] b = new int[9]; 2490 uint i; 2491 foreach (e; rc) 2492 { 2493 //writeln(e); 2494 b[i++] = e; 2495 } 2496 sort(b); 2497 assert(a == b, text(b)); 2498 } 2499} 2500 2501@safe unittest 2502{ 2503 // Bugzilla 12589 2504 int[] r = []; 2505 auto rc = randomCover(r); 2506 assert(rc.length == 0); 2507 assert(rc.empty); 2508 2509 // Bugzilla 16724 2510 import std.range : iota; 2511 auto range = iota(10); 2512 auto randy = range.randomCover; 2513 2514 for (int i=1; i <= range.length; i++) 2515 { 2516 randy.popFront; 2517 assert(randy.length == range.length - i); 2518 } 2519} 2520 2521// RandomSample 2522/** 2523Selects a random subsample out of $(D r), containing exactly $(D n) 2524elements. The order of elements is the same as in the original 2525range. The total length of $(D r) must be known. If $(D total) is 2526passed in, the total number of sample is considered to be $(D 2527total). Otherwise, $(D RandomSample) uses $(D r.length). 2528 2529Params: 2530 r = range to sample from 2531 n = number of elements to include in the sample; 2532 must be less than or equal to the total number 2533 of elements in $(D r) and/or the parameter 2534 $(D total) (if provided) 2535 total = (semi-optional) number of elements of $(D r) 2536 from which to select the sample (counting from 2537 the beginning); must be less than or equal to 2538 the total number of elements in $(D r) itself. 2539 May be omitted if $(D r) has the $(D .length) 2540 property and the sample is to be drawn from 2541 all elements of $(D r). 2542 rng = (optional) random number generator to use; 2543 if not specified, defaults to $(D rndGen) 2544 2545Returns: 2546 Range whose elements consist of a randomly selected subset of 2547 the elements of $(D r), in the same order as these elements 2548 appear in $(D r) itself. Will be a forward range if both $(D r) 2549 and $(D rng) are forward ranges, an input range otherwise. 2550 2551$(D RandomSample) implements Jeffrey Scott Vitter's Algorithm D 2552(see Vitter $(HTTP dx.doi.org/10.1145/358105.893, 1984), $(HTTP 2553dx.doi.org/10.1145/23002.23003, 1987)), which selects a sample 2554of size $(D n) in O(n) steps and requiring O(n) random variates, 2555regardless of the size of the data being sampled. The exception 2556to this is if traversing k elements on the input range is itself 2557an O(k) operation (e.g. when sampling lines from an input file), 2558in which case the sampling calculation will inevitably be of 2559O(total). 2560 2561RandomSample will throw an exception if $(D total) is verifiably 2562less than the total number of elements available in the input, 2563or if $(D n > total). 2564 2565If no random number generator is passed to $(D randomSample), the 2566thread-global RNG rndGen will be used internally. 2567 2568Example: 2569---- 2570int[] a = [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 ]; 2571// Print 5 random elements picked off from a 2572foreach (e; randomSample(a, 5)) 2573{ 2574 writeln(e); 2575} 2576---- 2577 2578$(B WARNING:) If an alternative RNG is desired, it is essential for this 2579to be a $(I new) RNG seeded in an unpredictable manner. Passing it a RNG 2580used elsewhere in the program will result in unintended correlations, 2581due to the current implementation of RNGs as value types. 2582 2583Example: 2584---- 2585int[] a = [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 ]; 2586foreach (e; randomSample(a, 5, Random(unpredictableSeed))) // correct! 2587{ 2588 writeln(e); 2589} 2590 2591foreach (e; randomSample(a, 5, rndGen)) // DANGEROUS!! rndGen gets 2592{ // copied by value 2593 writeln(e); 2594} 2595 2596foreach (e; randomSample(a, 5, rndGen)) // ... so this second random 2597{ // sample will select the same 2598 writeln(e); // values as the previous one. 2599} 2600---- 2601*/ 2602struct RandomSample(Range, UniformRNG = void) 2603if (isInputRange!Range && (isUniformRNG!UniformRNG || is(UniformRNG == void))) 2604{ 2605 private size_t _available, _toSelect; 2606 private enum ushort _alphaInverse = 13; // Vitter's recommended value. 2607 private double _Vprime; 2608 private Range _input; 2609 private size_t _index; 2610 private enum Skip { None, A, D } 2611 private Skip _skip = Skip.None; 2612 2613 // If we're using the default thread-local random number generator then 2614 // we shouldn't store a copy of it here. UniformRNG == void is a sentinel 2615 // for this. If we're using a user-specified generator then we have no 2616 // choice but to store a copy. 2617 static if (is(UniformRNG == void)) 2618 { 2619 static if (hasLength!Range) 2620 { 2621 this(Range input, size_t howMany) 2622 { 2623 _input = input; 2624 initialize(howMany, input.length); 2625 } 2626 } 2627 2628 this(Range input, size_t howMany, size_t total) 2629 { 2630 _input = input; 2631 initialize(howMany, total); 2632 } 2633 } 2634 else 2635 { 2636 UniformRNG _rng; 2637 2638 static if (hasLength!Range) 2639 { 2640 this(Range input, size_t howMany, ref UniformRNG rng) 2641 { 2642 _rng = rng; 2643 _input = input; 2644 initialize(howMany, input.length); 2645 } 2646 2647 this(Range input, size_t howMany, UniformRNG rng) 2648 { 2649 this(input, howMany, rng); 2650 } 2651 } 2652 2653 this(Range input, size_t howMany, size_t total, ref UniformRNG rng) 2654 { 2655 _rng = rng; 2656 _input = input; 2657 initialize(howMany, total); 2658 } 2659 2660 this(Range input, size_t howMany, size_t total, UniformRNG rng) 2661 { 2662 this(input, howMany, total, rng); 2663 } 2664 } 2665 2666 private void initialize(size_t howMany, size_t total) 2667 { 2668 import std.conv : text; 2669 import std.exception : enforce; 2670 _available = total; 2671 _toSelect = howMany; 2672 enforce(_toSelect <= _available, 2673 text("RandomSample: cannot sample ", _toSelect, 2674 " items when only ", _available, " are available")); 2675 static if (hasLength!Range) 2676 { 2677 enforce(_available <= _input.length, 2678 text("RandomSample: specified ", _available, 2679 " items as available when input contains only ", 2680 _input.length)); 2681 } 2682 } 2683 2684 private void initializeFront() 2685 { 2686 assert(_skip == Skip.None); 2687 // We can save ourselves a random variate by checking right 2688 // at the beginning if we should use Algorithm A. 2689 if ((_alphaInverse * _toSelect) > _available) 2690 { 2691 _skip = Skip.A; 2692 } 2693 else 2694 { 2695 _skip = Skip.D; 2696 _Vprime = newVprime(_toSelect); 2697 } 2698 prime(); 2699 } 2700 2701/** 2702 Range primitives. 2703*/ 2704 @property bool empty() const 2705 { 2706 return _toSelect == 0; 2707 } 2708 2709/// Ditto 2710 @property auto ref front() 2711 { 2712 assert(!empty); 2713 // The first sample point must be determined here to avoid 2714 // having it always correspond to the first element of the 2715 // input. The rest of the sample points are determined each 2716 // time we call popFront(). 2717 if (_skip == Skip.None) 2718 { 2719 initializeFront(); 2720 } 2721 return _input.front; 2722 } 2723 2724/// Ditto 2725 void popFront() 2726 { 2727 // First we need to check if the sample has 2728 // been initialized in the first place. 2729 if (_skip == Skip.None) 2730 { 2731 initializeFront(); 2732 } 2733 2734 _input.popFront(); 2735 --_available; 2736 --_toSelect; 2737 ++_index; 2738 prime(); 2739 } 2740 2741/// Ditto 2742 static if (isForwardRange!Range && isForwardRange!UniformRNG) 2743 { 2744 @property typeof(this) save() 2745 { 2746 auto ret = this; 2747 ret._input = _input.save; 2748 ret._rng = _rng.save; 2749 return ret; 2750 } 2751 } 2752 2753/// Ditto 2754 @property size_t length() 2755 { 2756 return _toSelect; 2757 } 2758 2759/** 2760Returns the index of the visited record. 2761 */ 2762 @property size_t index() 2763 { 2764 if (_skip == Skip.None) 2765 { 2766 initializeFront(); 2767 } 2768 return _index; 2769 } 2770 2771 private size_t skip() 2772 { 2773 assert(_skip != Skip.None); 2774 2775 // Step D1: if the number of points still to select is greater 2776 // than a certain proportion of the remaining data points, i.e. 2777 // if n >= alpha * N where alpha = 1/13, we carry out the 2778 // sampling with Algorithm A. 2779 if (_skip == Skip.A) 2780 { 2781 return skipA(); 2782 } 2783 else if ((_alphaInverse * _toSelect) > _available) 2784 { 2785 // We shouldn't get here unless the current selected 2786 // algorithm is D. 2787 assert(_skip == Skip.D); 2788 _skip = Skip.A; 2789 return skipA(); 2790 } 2791 else 2792 { 2793 assert(_skip == Skip.D); 2794 return skipD(); 2795 } 2796 } 2797 2798/* 2799Vitter's Algorithm A, used when the ratio of needed sample values 2800to remaining data values is sufficiently large. 2801*/ 2802 private size_t skipA() 2803 { 2804 size_t s; 2805 double v, quot, top; 2806 2807 if (_toSelect == 1) 2808 { 2809 static if (is(UniformRNG == void)) 2810 { 2811 s = uniform(0, _available); 2812 } 2813 else 2814 { 2815 s = uniform(0, _available, _rng); 2816 } 2817 } 2818 else 2819 { 2820 v = 0; 2821 top = _available - _toSelect; 2822 quot = top / _available; 2823 2824 static if (is(UniformRNG == void)) 2825 { 2826 v = uniform!"()"(0.0, 1.0); 2827 } 2828 else 2829 { 2830 v = uniform!"()"(0.0, 1.0, _rng); 2831 } 2832 2833 while (quot > v) 2834 { 2835 ++s; 2836 quot *= (top - s) / (_available - s); 2837 } 2838 } 2839 2840 return s; 2841 } 2842 2843/* 2844Randomly reset the value of _Vprime. 2845*/ 2846 private double newVprime(size_t remaining) 2847 { 2848 static if (is(UniformRNG == void)) 2849 { 2850 double r = uniform!"()"(0.0, 1.0); 2851 } 2852 else 2853 { 2854 double r = uniform!"()"(0.0, 1.0, _rng); 2855 } 2856 2857 return r ^^ (1.0 / remaining); 2858 } 2859 2860/* 2861Vitter's Algorithm D. For an extensive description of the algorithm 2862and its rationale, see: 2863 2864 * Vitter, J.S. (1984), "Faster methods for random sampling", 2865 Commun. ACM 27(7): 703--718 2866 2867 * Vitter, J.S. (1987) "An efficient algorithm for sequential random 2868 sampling", ACM Trans. Math. Softw. 13(1): 58-67. 2869 2870Variable names are chosen to match those in Vitter's paper. 2871*/ 2872 private size_t skipD() 2873 { 2874 import std.math : isNaN, trunc; 2875 // Confirm that the check in Step D1 is valid and we 2876 // haven't been sent here by mistake 2877 assert((_alphaInverse * _toSelect) <= _available); 2878 2879 // Now it's safe to use the standard Algorithm D mechanism. 2880 if (_toSelect > 1) 2881 { 2882 size_t s; 2883 size_t qu1 = 1 + _available - _toSelect; 2884 double x, y1; 2885 2886 assert(!_Vprime.isNaN()); 2887 2888 while (true) 2889 { 2890 // Step D2: set values of x and u. 2891 while (1) 2892 { 2893 x = _available * (1-_Vprime); 2894 s = cast(size_t) trunc(x); 2895 if (s < qu1) 2896 break; 2897 _Vprime = newVprime(_toSelect); 2898 } 2899 2900 static if (is(UniformRNG == void)) 2901 { 2902 double u = uniform!"()"(0.0, 1.0); 2903 } 2904 else 2905 { 2906 double u = uniform!"()"(0.0, 1.0, _rng); 2907 } 2908 2909 y1 = (u * (cast(double) _available) / qu1) ^^ (1.0/(_toSelect - 1)); 2910 2911 _Vprime = y1 * ((-x/_available)+1.0) * ( qu1/( (cast(double) qu1) - s ) ); 2912 2913 // Step D3: if _Vprime <= 1.0 our work is done and we return S. 2914 // Otherwise ... 2915 if (_Vprime > 1.0) 2916 { 2917 size_t top = _available - 1, limit; 2918 double y2 = 1.0, bottom; 2919 2920 if (_toSelect > (s+1)) 2921 { 2922 bottom = _available - _toSelect; 2923 limit = _available - s; 2924 } 2925 else 2926 { 2927 bottom = _available - (s+1); 2928 limit = qu1; 2929 } 2930 2931 foreach (size_t t; limit .. _available) 2932 { 2933 y2 *= top/bottom; 2934 top--; 2935 bottom--; 2936 } 2937 2938 // Step D4: decide whether or not to accept the current value of S. 2939 if (_available/(_available-x) < y1 * (y2 ^^ (1.0/(_toSelect-1)))) 2940 { 2941 // If it's not acceptable, we generate a new value of _Vprime 2942 // and go back to the start of the for (;;) loop. 2943 _Vprime = newVprime(_toSelect); 2944 } 2945 else 2946 { 2947 // If it's acceptable we generate a new value of _Vprime 2948 // based on the remaining number of sample points needed, 2949 // and return S. 2950 _Vprime = newVprime(_toSelect-1); 2951 return s; 2952 } 2953 } 2954 else 2955 { 2956 // Return if condition D3 satisfied. 2957 return s; 2958 } 2959 } 2960 } 2961 else 2962 { 2963 // If only one sample point remains to be taken ... 2964 return cast(size_t) trunc(_available * _Vprime); 2965 } 2966 } 2967 2968 private void prime() 2969 { 2970 if (empty) 2971 { 2972 return; 2973 } 2974 assert(_available && _available >= _toSelect); 2975 immutable size_t s = skip(); 2976 assert(s + _toSelect <= _available); 2977 static if (hasLength!Range) 2978 { 2979 assert(s + _toSelect <= _input.length); 2980 } 2981 assert(!_input.empty); 2982 _input.popFrontExactly(s); 2983 _index += s; 2984 _available -= s; 2985 assert(_available > 0); 2986 } 2987} 2988 2989/// Ditto 2990auto randomSample(Range)(Range r, size_t n, size_t total) 2991if (isInputRange!Range) 2992{ 2993 return RandomSample!(Range, void)(r, n, total); 2994} 2995 2996/// Ditto 2997auto randomSample(Range)(Range r, size_t n) 2998if (isInputRange!Range && hasLength!Range) 2999{ 3000 return RandomSample!(Range, void)(r, n, r.length); 3001} 3002 3003/// Ditto 3004auto randomSample(Range, UniformRNG)(Range r, size_t n, size_t total, auto ref UniformRNG rng) 3005if (isInputRange!Range && isUniformRNG!UniformRNG) 3006{ 3007 return RandomSample!(Range, UniformRNG)(r, n, total, rng); 3008} 3009 3010/// Ditto 3011auto randomSample(Range, UniformRNG)(Range r, size_t n, auto ref UniformRNG rng) 3012if (isInputRange!Range && hasLength!Range && isUniformRNG!UniformRNG) 3013{ 3014 return RandomSample!(Range, UniformRNG)(r, n, r.length, rng); 3015} 3016 3017@system unittest 3018{ 3019 // @system because it takes the address of a local 3020 import std.conv : text; 3021 import std.exception; 3022 import std.range; 3023 // For test purposes, an infinite input range 3024 struct TestInputRange 3025 { 3026 private auto r = recurrence!"a[n-1] + 1"(0); 3027 bool empty() @property const pure nothrow { return r.empty; } 3028 auto front() @property pure nothrow { return r.front; } 3029 void popFront() pure nothrow { r.popFront(); } 3030 } 3031 static assert(isInputRange!TestInputRange); 3032 static assert(!isForwardRange!TestInputRange); 3033 3034 int[] a = [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 ]; 3035 3036 foreach (UniformRNG; PseudoRngTypes) 3037 { 3038 auto rng = UniformRNG(1234); 3039 /* First test the most general case: randomSample of input range, with and 3040 * without a specified random number generator. 3041 */ 3042 static assert(isInputRange!(typeof(randomSample(TestInputRange(), 5, 10)))); 3043 static assert(isInputRange!(typeof(randomSample(TestInputRange(), 5, 10, rng)))); 3044 static assert(!isForwardRange!(typeof(randomSample(TestInputRange(), 5, 10)))); 3045 static assert(!isForwardRange!(typeof(randomSample(TestInputRange(), 5, 10, rng)))); 3046 // test case with range initialized by direct call to struct 3047 { 3048 auto sample = 3049 RandomSample!(TestInputRange, UniformRNG) 3050 (TestInputRange(), 5, 10, UniformRNG(unpredictableSeed)); 3051 static assert(isInputRange!(typeof(sample))); 3052 static assert(!isForwardRange!(typeof(sample))); 3053 } 3054 3055 /* Now test the case of an input range with length. We ignore the cases 3056 * already covered by the previous tests. 3057 */ 3058 static assert(isInputRange!(typeof(randomSample(TestInputRange().takeExactly(10), 5)))); 3059 static assert(isInputRange!(typeof(randomSample(TestInputRange().takeExactly(10), 5, rng)))); 3060 static assert(!isForwardRange!(typeof(randomSample(TestInputRange().takeExactly(10), 5)))); 3061 static assert(!isForwardRange!(typeof(randomSample(TestInputRange().takeExactly(10), 5, rng)))); 3062 // test case with range initialized by direct call to struct 3063 { 3064 auto sample = 3065 RandomSample!(typeof(TestInputRange().takeExactly(10)), UniformRNG) 3066 (TestInputRange().takeExactly(10), 5, 10, UniformRNG(unpredictableSeed)); 3067 static assert(isInputRange!(typeof(sample))); 3068 static assert(!isForwardRange!(typeof(sample))); 3069 } 3070 3071 // Now test the case of providing a forward range as input. 3072 static assert(!isForwardRange!(typeof(randomSample(a, 5)))); 3073 static if (isForwardRange!UniformRNG) 3074 { 3075 static assert(isForwardRange!(typeof(randomSample(a, 5, rng)))); 3076 // ... and test with range initialized directly 3077 { 3078 auto sample = 3079 RandomSample!(int[], UniformRNG) 3080 (a, 5, UniformRNG(unpredictableSeed)); 3081 static assert(isForwardRange!(typeof(sample))); 3082 } 3083 } 3084 else 3085 { 3086 static assert(isInputRange!(typeof(randomSample(a, 5, rng)))); 3087 static assert(!isForwardRange!(typeof(randomSample(a, 5, rng)))); 3088 // ... and test with range initialized directly 3089 { 3090 auto sample = 3091 RandomSample!(int[], UniformRNG) 3092 (a, 5, UniformRNG(unpredictableSeed)); 3093 static assert(isInputRange!(typeof(sample))); 3094 static assert(!isForwardRange!(typeof(sample))); 3095 } 3096 } 3097 3098 /* Check that randomSample will throw an error if we claim more 3099 * items are available than there actually are, or if we try to 3100 * sample more items than are available. */ 3101 assert(collectExceptionMsg( 3102 randomSample(a, 5, 15) 3103 ) == "RandomSample: specified 15 items as available when input contains only 10"); 3104 assert(collectExceptionMsg( 3105 randomSample(a, 15) 3106 ) == "RandomSample: cannot sample 15 items when only 10 are available"); 3107 assert(collectExceptionMsg( 3108 randomSample(a, 9, 8) 3109 ) == "RandomSample: cannot sample 9 items when only 8 are available"); 3110 assert(collectExceptionMsg( 3111 randomSample(TestInputRange(), 12, 11) 3112 ) == "RandomSample: cannot sample 12 items when only 11 are available"); 3113 3114 /* Check that sampling algorithm never accidentally overruns the end of 3115 * the input range. If input is an InputRange without .length, this 3116 * relies on the user specifying the total number of available items 3117 * correctly. 3118 */ 3119 { 3120 uint i = 0; 3121 foreach (e; randomSample(a, a.length)) 3122 { 3123 assert(e == i); 3124 ++i; 3125 } 3126 assert(i == a.length); 3127 3128 i = 0; 3129 foreach (e; randomSample(TestInputRange(), 17, 17)) 3130 { 3131 assert(e == i); 3132 ++i; 3133 } 3134 assert(i == 17); 3135 } 3136 3137 3138 // Check length properties of random samples. 3139 assert(randomSample(a, 5).length == 5); 3140 assert(randomSample(a, 5, 10).length == 5); 3141 assert(randomSample(a, 5, rng).length == 5); 3142 assert(randomSample(a, 5, 10, rng).length == 5); 3143 assert(randomSample(TestInputRange(), 5, 10).length == 5); 3144 assert(randomSample(TestInputRange(), 5, 10, rng).length == 5); 3145 3146 // ... and emptiness! 3147 assert(randomSample(a, 0).empty); 3148 assert(randomSample(a, 0, 5).empty); 3149 assert(randomSample(a, 0, rng).empty); 3150 assert(randomSample(a, 0, 5, rng).empty); 3151 assert(randomSample(TestInputRange(), 0, 10).empty); 3152 assert(randomSample(TestInputRange(), 0, 10, rng).empty); 3153 3154 /* Test that the (lazy) evaluation of random samples works correctly. 3155 * 3156 * We cover 2 different cases: a sample where the ratio of sample points 3157 * to total points is greater than the threshold for using Algorithm, and 3158 * one where the ratio is small enough (< 1/13) for Algorithm D to be used. 3159 * 3160 * For each, we also cover the case with and without a specified RNG. 3161 */ 3162 { 3163 // Small sample/source ratio, no specified RNG. 3164 uint i = 0; 3165 foreach (e; randomSample(randomCover(a), 5)) 3166 { 3167 ++i; 3168 } 3169 assert(i == 5); 3170 3171 // Small sample/source ratio, specified RNG. 3172 i = 0; 3173 foreach (e; randomSample(randomCover(a), 5, rng)) 3174 { 3175 ++i; 3176 } 3177 assert(i == 5); 3178 3179 // Large sample/source ratio, no specified RNG. 3180 i = 0; 3181 foreach (e; randomSample(TestInputRange(), 123, 123_456)) 3182 { 3183 ++i; 3184 } 3185 assert(i == 123); 3186 3187 // Large sample/source ratio, specified RNG. 3188 i = 0; 3189 foreach (e; randomSample(TestInputRange(), 123, 123_456, rng)) 3190 { 3191 ++i; 3192 } 3193 assert(i == 123); 3194 3195 /* Sample/source ratio large enough to start with Algorithm D, 3196 * small enough to switch to Algorithm A. 3197 */ 3198 i = 0; 3199 foreach (e; randomSample(TestInputRange(), 10, 131)) 3200 { 3201 ++i; 3202 } 3203 assert(i == 10); 3204 } 3205 3206 // Test that the .index property works correctly 3207 { 3208 auto sample1 = randomSample(TestInputRange(), 654, 654_321); 3209 for (; !sample1.empty; sample1.popFront()) 3210 { 3211 assert(sample1.front == sample1.index); 3212 } 3213 3214 auto sample2 = randomSample(TestInputRange(), 654, 654_321, rng); 3215 for (; !sample2.empty; sample2.popFront()) 3216 { 3217 assert(sample2.front == sample2.index); 3218 } 3219 3220 /* Check that it also works if .index is called before .front. 3221 * See: http://d.puremagic.com/issues/show_bug.cgi?id=10322 3222 */ 3223 auto sample3 = randomSample(TestInputRange(), 654, 654_321); 3224 for (; !sample3.empty; sample3.popFront()) 3225 { 3226 assert(sample3.index == sample3.front); 3227 } 3228 3229 auto sample4 = randomSample(TestInputRange(), 654, 654_321, rng); 3230 for (; !sample4.empty; sample4.popFront()) 3231 { 3232 assert(sample4.index == sample4.front); 3233 } 3234 } 3235 3236 /* Test behaviour if .popFront() is called before sample is read. 3237 * This is a rough-and-ready check that the statistical properties 3238 * are in the ballpark -- not a proper validation of statistical 3239 * quality! This incidentally also checks for reference-type 3240 * initialization bugs, as the foreach () loop will operate on a 3241 * copy of the popFronted (and hence initialized) sample. 3242 */ 3243 { 3244 size_t count0, count1, count99; 3245 foreach (_; 0 .. 100_000) 3246 { 3247 auto sample = randomSample(iota(100), 5, &rng); 3248 sample.popFront(); 3249 foreach (s; sample) 3250 { 3251 if (s == 0) 3252 { 3253 ++count0; 3254 } 3255 else if (s == 1) 3256 { 3257 ++count1; 3258 } 3259 else if (s == 99) 3260 { 3261 ++count99; 3262 } 3263 } 3264 } 3265 /* Statistical assumptions here: this is a sequential sampling process 3266 * so (i) 0 can only be the first sample point, so _can't_ be in the 3267 * remainder of the sample after .popFront() is called. (ii) By similar 3268 * token, 1 can only be in the remainder if it's the 2nd point of the 3269 * whole sample, and hence if 0 was the first; probability of 0 being 3270 * first and 1 second is 5/100 * 4/99 (thank you, Algorithm S:-) and 3271 * so the mean count of 1 should be about 202. Finally, 99 can only 3272 * be the _last_ sample point to be picked, so its probability of 3273 * inclusion should be independent of the .popFront() and it should 3274 * occur with frequency 5/100, hence its count should be about 5000. 3275 * Unfortunately we have to set quite a high tolerance because with 3276 * sample size small enough for unittests to run in reasonable time, 3277 * the variance can be quite high. 3278 */ 3279 assert(count0 == 0); 3280 assert(count1 < 300, text("1: ", count1, " > 300.")); 3281 assert(4_700 < count99, text("99: ", count99, " < 4700.")); 3282 assert(count99 < 5_300, text("99: ", count99, " > 5300.")); 3283 } 3284 3285 /* Odd corner-cases: RandomSample has 2 constructors that are not called 3286 * by the randomSample() helper functions, but that can be used if the 3287 * constructor is called directly. These cover the case of the user 3288 * specifying input but not input length. 3289 */ 3290 { 3291 auto input1 = TestInputRange().takeExactly(456_789); 3292 static assert(hasLength!(typeof(input1))); 3293 auto sample1 = RandomSample!(typeof(input1), void)(input1, 789); 3294 static assert(isInputRange!(typeof(sample1))); 3295 static assert(!isForwardRange!(typeof(sample1))); 3296 assert(sample1.length == 789); 3297 assert(sample1._available == 456_789); 3298 uint i = 0; 3299 for (; !sample1.empty; sample1.popFront()) 3300 { 3301 assert(sample1.front == sample1.index); 3302 ++i; 3303 } 3304 assert(i == 789); 3305 3306 auto input2 = TestInputRange().takeExactly(456_789); 3307 static assert(hasLength!(typeof(input2))); 3308 auto sample2 = RandomSample!(typeof(input2), typeof(rng))(input2, 789, rng); 3309 static assert(isInputRange!(typeof(sample2))); 3310 static assert(!isForwardRange!(typeof(sample2))); 3311 assert(sample2.length == 789); 3312 assert(sample2._available == 456_789); 3313 i = 0; 3314 for (; !sample2.empty; sample2.popFront()) 3315 { 3316 assert(sample2.front == sample2.index); 3317 ++i; 3318 } 3319 assert(i == 789); 3320 } 3321 3322 /* Test that the save property works where input is a forward range, 3323 * and RandomSample is using a (forward range) random number generator 3324 * that is not rndGen. 3325 */ 3326 static if (isForwardRange!UniformRNG) 3327 { 3328 auto sample1 = randomSample(a, 5, rng); 3329 auto sample2 = sample1.save; 3330 assert(sample1.array() == sample2.array()); 3331 } 3332 3333 // Bugzilla 8314 3334 { 3335 auto sample(RandomGen)(uint seed) { return randomSample(a, 1, RandomGen(seed)).front; } 3336 3337 // Start from 1 because not all RNGs accept 0 as seed. 3338 immutable fst = sample!UniformRNG(1); 3339 uint n = 1; 3340 while (sample!UniformRNG(++n) == fst && n < n.max) {} 3341 assert(n < n.max); 3342 } 3343 } 3344} 3345