d: Merge upstream dmd 3982604c5, druntime bc58b1e9, phobos 12329adb6.
[official-gcc.git] / libphobos / src / std / random.d
blob106e51ceedb193eb62dd9a17e7cc4a5a1edd6a6d
1 // Written in the D programming language.
3 /**
4 Facilities for random number generation.
6 $(SCRIPT inhibitQuickIndex = 1;)
7 $(DIVC quickindex,
8 $(BOOKTABLE,
9 $(TR $(TH Category) $(TH Functions))
10 $(TR $(TD Uniform sampling) $(TD
11 $(LREF uniform)
12 $(LREF uniform01)
13 $(LREF uniformDistribution)
15 $(TR $(TD Element sampling) $(TD
16 $(LREF choice)
17 $(LREF dice)
19 $(TR $(TD Range sampling) $(TD
20 $(LREF randomCover)
21 $(LREF randomSample)
23 $(TR $(TD Default Random Engines) $(TD
24 $(LREF rndGen)
25 $(LREF Random)
26 $(LREF unpredictableSeed)
28 $(TR $(TD Linear Congruential Engines) $(TD
29 $(LREF MinstdRand)
30 $(LREF MinstdRand0)
31 $(LREF LinearCongruentialEngine)
33 $(TR $(TD Mersenne Twister Engines) $(TD
34 $(LREF Mt19937)
35 $(LREF Mt19937_64)
36 $(LREF MersenneTwisterEngine)
38 $(TR $(TD Xorshift Engines) $(TD
39 $(LREF Xorshift)
40 $(LREF XorshiftEngine)
41 $(LREF Xorshift32)
42 $(LREF Xorshift64)
43 $(LREF Xorshift96)
44 $(LREF Xorshift128)
45 $(LREF Xorshift160)
46 $(LREF Xorshift192)
48 $(TR $(TD Shuffle) $(TD
49 $(LREF partialShuffle)
50 $(LREF randomShuffle)
52 $(TR $(TD Traits) $(TD
53 $(LREF isSeedable)
54 $(LREF isUniformRNG)
58 $(RED Disclaimer:) The random number generators and API provided in this
59 module are not designed to be cryptographically secure, and are therefore
60 unsuitable for cryptographic or security-related purposes such as generating
61 authentication tokens or network sequence numbers. For such needs, please use a
62 reputable cryptographic library instead.
64 The new-style generator objects hold their own state so they are
65 immune of threading issues. The generators feature a number of
66 well-known and well-documented methods of generating random
67 numbers. An overall fast and reliable means to generate random numbers
68 is the $(D_PARAM Mt19937) generator, which derives its name from
69 "$(LINK2 https://en.wikipedia.org/wiki/Mersenne_Twister, Mersenne Twister)
70 with a period of 2 to the power of
71 19937". In memory-constrained situations,
72 $(LINK2 https://en.wikipedia.org/wiki/Linear_congruential_generator,
73 linear congruential generators) such as `MinstdRand0` and `MinstdRand` might be
74 useful. The standard library provides an alias $(D_PARAM Random) for
75 whichever generator it considers the most fit for the target
76 environment.
78 In addition to random number generators, this module features
79 distributions, which skew a generator's output statistical
80 distribution in various ways. So far the uniform distribution for
81 integers and real numbers have been implemented.
83 Source: $(PHOBOSSRC std/random.d)
85 Macros:
87 Copyright: Copyright Andrei Alexandrescu 2008 - 2009, Joseph Rushton Wakeling 2012.
88 License: $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0).
89 Authors: $(HTTP erdani.org, Andrei Alexandrescu)
90 Masahiro Nakagawa (Xorshift random generator)
91 $(HTTP braingam.es, Joseph Rushton Wakeling) (Algorithm D for random sampling)
92 Ilya Yaroshenko (Mersenne Twister implementation, adapted from $(HTTPS github.com/libmir/mir-random, mir-random))
93 Credits: The entire random number library architecture is derived from the
94 excellent $(HTTP open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2461.pdf, C++0X)
95 random number facility proposed by Jens Maurer and contributed to by
96 researchers at the Fermi laboratory (excluding Xorshift).
99 Copyright Andrei Alexandrescu 2008 - 2009.
100 Distributed under the Boost Software License, Version 1.0.
101 (See accompanying file LICENSE_1_0.txt or copy at
102 http://www.boost.org/LICENSE_1_0.txt)
104 module std.random;
107 import std.range.primitives;
108 import std.traits;
110 version (OSX)
111 version = Darwin;
112 else version (iOS)
113 version = Darwin;
114 else version (TVOS)
115 version = Darwin;
116 else version (WatchOS)
117 version = Darwin;
119 version (D_InlineAsm_X86) version = InlineAsm_X86_Any;
120 version (D_InlineAsm_X86_64) version = InlineAsm_X86_Any;
123 @safe unittest
125 import std.algorithm.comparison : among, equal;
126 import std.range : iota;
128 // seed a random generator with a constant
129 auto rnd = Random(42);
131 // Generate a uniformly-distributed integer in the range [0, 14]
132 // If no random generator is passed, the global `rndGen` would be used
133 auto i = uniform(0, 15, rnd);
134 assert(i >= 0 && i < 15);
136 // Generate a uniformly-distributed real in the range [0, 100)
137 auto r = uniform(0.0L, 100.0L, rnd);
138 assert(r >= 0 && r < 100);
140 // Sample from a custom type
141 enum Fruit { apple, mango, pear }
142 auto f = rnd.uniform!Fruit;
143 with(Fruit)
144 assert(f.among(apple, mango, pear));
146 // Generate a 32-bit random number
147 auto u = uniform!uint(rnd);
148 static assert(is(typeof(u) == uint));
150 // Generate a random number in the range in the range [0, 1)
151 auto u2 = uniform01(rnd);
152 assert(u2 >= 0 && u2 < 1);
154 // Select an element randomly
155 auto el = 10.iota.choice(rnd);
156 assert(0 <= el && el < 10);
158 // Throw a dice with custom proportions
159 // 0: 20%, 1: 10%, 2: 60%
160 auto val = rnd.dice(0.2, 0.1, 0.6);
161 assert(0 <= val && val <= 2);
163 auto rnd2 = MinstdRand0(42);
165 // Select a random subsample from a range
166 assert(10.iota.randomSample(3, rnd2).equal([7, 8, 9]));
168 // Cover all elements in an array in random order
169 version (D_LP64) // https://issues.dlang.org/show_bug.cgi?id=15147
170 assert(10.iota.randomCover(rnd2).equal([7, 4, 2, 0, 1, 6, 8, 3, 9, 5]));
171 else
172 assert(10.iota.randomCover(rnd2).equal([4, 8, 7, 3, 5, 9, 2, 6, 0, 1]));
174 // Shuffle an array
175 version (D_LP64) // https://issues.dlang.org/show_bug.cgi?id=15147
176 assert([0, 1, 2, 4, 5].randomShuffle(rnd2).equal([2, 0, 4, 5, 1]));
177 else
178 assert([0, 1, 2, 4, 5].randomShuffle(rnd2).equal([4, 2, 5, 0, 1]));
181 version (StdUnittest)
183 static import std.meta;
184 package alias Xorshift64_64 = XorshiftEngine!(ulong, 64, -12, 25, -27);
185 package alias Xorshift128_64 = XorshiftEngine!(ulong, 128, 23, -18, -5);
186 package alias PseudoRngTypes = std.meta.AliasSeq!(MinstdRand0, MinstdRand, Mt19937, Xorshift32, Xorshift64,
187 Xorshift96, Xorshift128, Xorshift160, Xorshift192,
188 Xorshift64_64, Xorshift128_64);
191 // Segments of the code in this file Copyright (c) 1997 by Rick Booth
192 // From "Inner Loops" by Rick Booth, Addison-Wesley
194 // Work derived from:
197 A C-program for MT19937, with initialization improved 2002/1/26.
198 Coded by Takuji Nishimura and Makoto Matsumoto.
200 Before using, initialize the state by using init_genrand(seed)
201 or init_by_array(init_key, key_length).
203 Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
204 All rights reserved.
206 Redistribution and use in source and binary forms, with or without
207 modification, are permitted provided that the following conditions
208 are met:
210 1. Redistributions of source code must retain the above copyright
211 notice, this list of conditions and the following disclaimer.
213 2. Redistributions in binary form must reproduce the above copyright
214 notice, this list of conditions and the following disclaimer in the
215 documentation and/or other materials provided with the distribution.
217 3. The names of its contributors may not be used to endorse or promote
218 products derived from this software without specific prior written
219 permission.
221 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
222 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
223 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
224 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
225 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
226 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
227 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
228 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
229 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
230 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
231 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
234 Any feedback is very welcome.
235 http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
236 email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
240 * Test if Rng is a random-number generator. The overload
241 * taking a ElementType also makes sure that the Rng generates
242 * values of that type.
244 * A random-number generator has at least the following features:
245 * $(UL
246 * $(LI it's an InputRange)
247 * $(LI it has a 'bool isUniformRandom' field readable in CTFE)
250 template isUniformRNG(Rng, ElementType)
252 enum bool isUniformRNG = .isUniformRNG!Rng &&
253 is(std.range.primitives.ElementType!Rng == ElementType);
257 * ditto
259 template isUniformRNG(Rng)
261 enum bool isUniformRNG = isInputRange!Rng &&
262 is(typeof(
264 static assert(Rng.isUniformRandom); //tag
265 }));
269 @safe unittest
271 struct NoRng
273 @property uint front() {return 0;}
274 @property bool empty() {return false;}
275 void popFront() {}
277 static assert(!isUniformRNG!(NoRng));
279 struct validRng
281 @property uint front() {return 0;}
282 @property bool empty() {return false;}
283 void popFront() {}
285 enum isUniformRandom = true;
287 static assert(isUniformRNG!(validRng, uint));
288 static assert(isUniformRNG!(validRng));
291 @safe unittest
293 // two-argument predicate should not require @property on `front`
294 // https://issues.dlang.org/show_bug.cgi?id=19837
295 struct Rng
297 float front() {return 0;}
298 void popFront() {}
299 enum empty = false;
300 enum isUniformRandom = true;
302 static assert(isUniformRNG!(Rng, float));
306 * Test if Rng is seedable. The overload
307 * taking a SeedType also makes sure that the Rng can be seeded with SeedType.
309 * A seedable random-number generator has the following additional features:
310 * $(UL
311 * $(LI it has a 'seed(ElementType)' function)
314 template isSeedable(Rng, SeedType)
316 enum bool isSeedable = isUniformRNG!(Rng) &&
317 is(typeof(
319 Rng r = void; // can define a Rng object
320 SeedType s = void;
321 r.seed(s); // can seed a Rng
322 }));
325 ///ditto
326 template isSeedable(Rng)
328 enum bool isSeedable = isUniformRNG!Rng &&
329 is(typeof(
331 Rng r = void; // can define a Rng object
332 alias SeedType = typeof(r.front);
333 SeedType s = void;
334 r.seed(s); // can seed a Rng
335 }));
339 @safe unittest
341 struct validRng
343 @property uint front() {return 0;}
344 @property bool empty() {return false;}
345 void popFront() {}
347 enum isUniformRandom = true;
349 static assert(!isSeedable!(validRng, uint));
350 static assert(!isSeedable!(validRng));
352 struct seedRng
354 @property uint front() {return 0;}
355 @property bool empty() {return false;}
356 void popFront() {}
357 void seed(uint val){}
358 enum isUniformRandom = true;
360 static assert(isSeedable!(seedRng, uint));
361 static assert(!isSeedable!(seedRng, ulong));
362 static assert(isSeedable!(seedRng));
365 @safe @nogc pure nothrow unittest
367 struct NoRng
369 @property uint front() {return 0;}
370 @property bool empty() {return false;}
371 void popFront() {}
373 static assert(!isUniformRNG!(NoRng, uint));
374 static assert(!isUniformRNG!(NoRng));
375 static assert(!isSeedable!(NoRng, uint));
376 static assert(!isSeedable!(NoRng));
378 struct NoRng2
380 @property uint front() {return 0;}
381 @property bool empty() {return false;}
382 void popFront() {}
384 enum isUniformRandom = false;
386 static assert(!isUniformRNG!(NoRng2, uint));
387 static assert(!isUniformRNG!(NoRng2));
388 static assert(!isSeedable!(NoRng2, uint));
389 static assert(!isSeedable!(NoRng2));
391 struct NoRng3
393 @property bool empty() {return false;}
394 void popFront() {}
396 enum isUniformRandom = true;
398 static assert(!isUniformRNG!(NoRng3, uint));
399 static assert(!isUniformRNG!(NoRng3));
400 static assert(!isSeedable!(NoRng3, uint));
401 static assert(!isSeedable!(NoRng3));
403 struct validRng
405 @property uint front() {return 0;}
406 @property bool empty() {return false;}
407 void popFront() {}
409 enum isUniformRandom = true;
411 static assert(isUniformRNG!(validRng, uint));
412 static assert(isUniformRNG!(validRng));
413 static assert(!isSeedable!(validRng, uint));
414 static assert(!isSeedable!(validRng));
416 struct seedRng
418 @property uint front() {return 0;}
419 @property bool empty() {return false;}
420 void popFront() {}
421 void seed(uint val){}
422 enum isUniformRandom = true;
424 static assert(isUniformRNG!(seedRng, uint));
425 static assert(isUniformRNG!(seedRng));
426 static assert(isSeedable!(seedRng, uint));
427 static assert(!isSeedable!(seedRng, ulong));
428 static assert(isSeedable!(seedRng));
432 Linear Congruential generator. When m = 0, no modulus is used.
434 struct LinearCongruentialEngine(UIntType, UIntType a, UIntType c, UIntType m)
435 if (isUnsigned!UIntType)
437 ///Mark this as a Rng
438 enum bool isUniformRandom = true;
439 /// Does this generator have a fixed range? ($(D_PARAM true)).
440 enum bool hasFixedRange = true;
441 /// Lowest generated value (`1` if $(D c == 0), `0` otherwise).
442 enum UIntType min = ( c == 0 ? 1 : 0 );
443 /// Highest generated value ($(D modulus - 1)).
444 enum UIntType max = m - 1;
446 The parameters of this distribution. The random number is $(D_PARAM x
447 = (x * multipler + increment) % modulus).
449 enum UIntType multiplier = a;
450 ///ditto
451 enum UIntType increment = c;
452 ///ditto
453 enum UIntType modulus = m;
455 static assert(isIntegral!(UIntType));
456 static assert(m == 0 || a < m);
457 static assert(m == 0 || c < m);
458 static assert(m == 0 ||
459 (cast(ulong) a * (m-1) + c) % m == (c < a ? c - a + m : c - a));
461 // Check for maximum range
462 private static ulong gcd(ulong a, ulong b) @safe pure nothrow @nogc
464 while (b)
466 auto t = b;
467 b = a % b;
468 a = t;
470 return a;
473 private static ulong primeFactorsOnly(ulong n) @safe pure nothrow @nogc
475 ulong result = 1;
476 ulong iter = 2;
477 for (; n >= iter * iter; iter += 2 - (iter == 2))
479 if (n % iter) continue;
480 result *= iter;
483 n /= iter;
484 } while (n % iter == 0);
486 return result * n;
489 @safe pure nothrow unittest
491 static assert(primeFactorsOnly(100) == 10);
492 //writeln(primeFactorsOnly(11));
493 static assert(primeFactorsOnly(11) == 11);
494 static assert(primeFactorsOnly(7 * 7 * 7 * 11 * 15 * 11) == 7 * 11 * 15);
495 static assert(primeFactorsOnly(129 * 2) == 129 * 2);
496 // enum x = primeFactorsOnly(7 * 7 * 7 * 11 * 15);
497 // static assert(x == 7 * 11 * 15);
500 private static bool properLinearCongruentialParameters(ulong m,
501 ulong a, ulong c) @safe pure nothrow @nogc
503 if (m == 0)
505 static if (is(UIntType == uint))
507 // Assume m is uint.max + 1
508 m = (1uL << 32);
510 else
512 return false;
515 // Bounds checking
516 if (a == 0 || a >= m || c >= m) return false;
517 // c and m are relatively prime
518 if (c > 0 && gcd(c, m) != 1) return false;
519 // a - 1 is divisible by all prime factors of m
520 if ((a - 1) % primeFactorsOnly(m)) return false;
521 // if a - 1 is multiple of 4, then m is a multiple of 4 too.
522 if ((a - 1) % 4 == 0 && m % 4) return false;
523 // Passed all tests
524 return true;
527 // check here
528 static assert(c == 0 || properLinearCongruentialParameters(m, a, c),
529 "Incorrect instantiation of LinearCongruentialEngine");
532 Constructs a $(D_PARAM LinearCongruentialEngine) generator seeded with
533 `x0`.
535 this(UIntType x0) @safe pure nothrow @nogc
537 seed(x0);
541 (Re)seeds the generator.
543 void seed(UIntType x0 = 1) @safe pure nothrow @nogc
545 _x = modulus ? (x0 % modulus) : x0;
546 static if (c == 0)
548 //Necessary to prevent generator from outputting an endless series of zeroes.
549 if (_x == 0)
550 _x = max;
552 popFront();
556 Advances the random sequence.
558 void popFront() @safe pure nothrow @nogc
560 static if (m)
562 static if (is(UIntType == uint) && m == uint.max)
564 immutable ulong
565 x = (cast(ulong) a * _x + c),
566 v = x >> 32,
567 w = x & uint.max;
568 immutable y = cast(uint)(v + w);
569 _x = (y < v || y == uint.max) ? (y + 1) : y;
571 else static if (is(UIntType == uint) && m == int.max)
573 immutable ulong
574 x = (cast(ulong) a * _x + c),
575 v = x >> 31,
576 w = x & int.max;
577 immutable uint y = cast(uint)(v + w);
578 _x = (y >= int.max) ? (y - int.max) : y;
580 else
582 _x = cast(UIntType) ((cast(ulong) a * _x + c) % m);
585 else
587 _x = a * _x + c;
592 Returns the current number in the random sequence.
594 @property UIntType front() const @safe pure nothrow @nogc
596 return _x;
600 @property typeof(this) save() const @safe pure nothrow @nogc
602 return this;
606 Always `false` (random generators are infinite ranges).
608 enum bool empty = false;
610 // https://issues.dlang.org/show_bug.cgi?id=21610
611 static if (m)
613 private UIntType _x = (a + c) % m;
615 else
617 private UIntType _x = a + c;
621 /// Declare your own linear congruential engine
622 @safe unittest
624 alias CPP11LCG = LinearCongruentialEngine!(uint, 48271, 0, 2_147_483_647);
626 // seed with a constant
627 auto rnd = CPP11LCG(42);
628 auto n = rnd.front; // same for each run
629 assert(n == 2027382);
632 /// Declare your own linear congruential engine
633 @safe unittest
635 // glibc's LCG
636 alias GLibcLCG = LinearCongruentialEngine!(uint, 1103515245, 12345, 2_147_483_648);
638 // Seed with an unpredictable value
639 auto rnd = GLibcLCG(unpredictableSeed);
640 auto n = rnd.front; // different across runs
643 /// Declare your own linear congruential engine
644 @safe unittest
646 // Visual C++'s LCG
647 alias MSVCLCG = LinearCongruentialEngine!(uint, 214013, 2531011, 0);
649 // seed with a constant
650 auto rnd = MSVCLCG(1);
651 auto n = rnd.front; // same for each run
652 assert(n == 2745024);
655 // Ensure that unseeded LCGs produce correct values
656 @safe unittest
658 alias LGE = LinearCongruentialEngine!(uint, 10000, 19682, 19683);
660 LGE rnd;
661 assert(rnd.front == 9999);
665 Define $(D_PARAM LinearCongruentialEngine) generators with well-chosen
666 parameters. `MinstdRand0` implements Park and Miller's "minimal
667 standard" $(HTTP
668 wikipedia.org/wiki/Park%E2%80%93Miller_random_number_generator,
669 generator) that uses 16807 for the multiplier. `MinstdRand`
670 implements a variant that has slightly better spectral behavior by
671 using the multiplier 48271. Both generators are rather simplistic.
673 alias MinstdRand0 = LinearCongruentialEngine!(uint, 16_807, 0, 2_147_483_647);
674 /// ditto
675 alias MinstdRand = LinearCongruentialEngine!(uint, 48_271, 0, 2_147_483_647);
678 @safe @nogc unittest
680 // seed with a constant
681 auto rnd0 = MinstdRand0(1);
682 auto n = rnd0.front;
683 // same for each run
684 assert(n == 16807);
686 // Seed with an unpredictable value
687 rnd0.seed(unpredictableSeed);
688 n = rnd0.front; // different across runs
691 @safe @nogc unittest
693 import std.range;
694 static assert(isForwardRange!MinstdRand);
695 static assert(isUniformRNG!MinstdRand);
696 static assert(isUniformRNG!MinstdRand0);
697 static assert(isUniformRNG!(MinstdRand, uint));
698 static assert(isUniformRNG!(MinstdRand0, uint));
699 static assert(isSeedable!MinstdRand);
700 static assert(isSeedable!MinstdRand0);
701 static assert(isSeedable!(MinstdRand, uint));
702 static assert(isSeedable!(MinstdRand0, uint));
704 // The correct numbers are taken from The Database of Integer Sequences
705 // http://www.research.att.com/~njas/sequences/eisBTfry00128.txt
706 enum ulong[20] checking0 = [
707 16807UL,282475249,1622650073,984943658,1144108930,470211272,
708 101027544,1457850878,1458777923,2007237709,823564440,1115438165,
709 1784484492,74243042,114807987,1137522503,1441282327,16531729,
710 823378840,143542612 ];
711 //auto rnd0 = MinstdRand0(1);
712 MinstdRand0 rnd0;
714 foreach (e; checking0)
716 assert(rnd0.front == e);
717 rnd0.popFront();
719 // Test the 10000th invocation
720 // Correct value taken from:
721 // http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2461.pdf
722 rnd0.seed();
723 popFrontN(rnd0, 9999);
724 assert(rnd0.front == 1043618065);
726 // Test MinstdRand
727 enum ulong[6] checking = [48271UL,182605794,1291394886,1914720637,2078669041,
728 407355683];
729 //auto rnd = MinstdRand(1);
730 MinstdRand rnd;
731 foreach (e; checking)
733 assert(rnd.front == e);
734 rnd.popFront();
737 // Test the 10000th invocation
738 // Correct value taken from:
739 // http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2461.pdf
740 rnd.seed();
741 popFrontN(rnd, 9999);
742 assert(rnd.front == 399268537);
744 // Check .save works
745 static foreach (Type; std.meta.AliasSeq!(MinstdRand0, MinstdRand))
747 auto rnd1 = Type(123_456_789);
748 rnd1.popFront();
749 // https://issues.dlang.org/show_bug.cgi?id=15853
750 auto rnd2 = ((const ref Type a) => a.save())(rnd1);
751 assert(rnd1 == rnd2);
752 // Enable next test when RNGs are reference types
753 version (none) { assert(rnd1 !is rnd2); }
754 for (auto i = 0; i < 3; i++, rnd1.popFront, rnd2.popFront)
755 assert(rnd1.front() == rnd2.front());
759 @safe @nogc unittest
761 auto rnd0 = MinstdRand0(MinstdRand0.modulus);
762 auto n = rnd0.front;
763 rnd0.popFront();
764 assert(n != rnd0.front);
768 The $(LINK2 https://en.wikipedia.org/wiki/Mersenne_Twister, Mersenne Twister) generator.
770 struct MersenneTwisterEngine(UIntType, size_t w, size_t n, size_t m, size_t r,
771 UIntType a, size_t u, UIntType d, size_t s,
772 UIntType b, size_t t,
773 UIntType c, size_t l, UIntType f)
774 if (isUnsigned!UIntType)
776 static assert(0 < w && w <= UIntType.sizeof * 8);
777 static assert(1 <= m && m <= n);
778 static assert(0 <= r && 0 <= u && 0 <= s && 0 <= t && 0 <= l);
779 static assert(r <= w && u <= w && s <= w && t <= w && l <= w);
780 static assert(0 <= a && 0 <= b && 0 <= c);
781 static assert(n <= ptrdiff_t.max);
783 ///Mark this as a Rng
784 enum bool isUniformRandom = true;
787 Parameters for the generator.
789 enum size_t wordSize = w;
790 enum size_t stateSize = n; /// ditto
791 enum size_t shiftSize = m; /// ditto
792 enum size_t maskBits = r; /// ditto
793 enum UIntType xorMask = a; /// ditto
794 enum size_t temperingU = u; /// ditto
795 enum UIntType temperingD = d; /// ditto
796 enum size_t temperingS = s; /// ditto
797 enum UIntType temperingB = b; /// ditto
798 enum size_t temperingT = t; /// ditto
799 enum UIntType temperingC = c; /// ditto
800 enum size_t temperingL = l; /// ditto
801 enum UIntType initializationMultiplier = f; /// ditto
803 /// Smallest generated value (0).
804 enum UIntType min = 0;
805 /// Largest generated value.
806 enum UIntType max = UIntType.max >> (UIntType.sizeof * 8u - w);
807 // note, `max` also serves as a bitmask for the lowest `w` bits
808 static assert(a <= max && b <= max && c <= max && f <= max);
810 /// The default seed value.
811 enum UIntType defaultSeed = 5489u;
813 // Bitmasks used in the 'twist' part of the algorithm
814 private enum UIntType lowerMask = (cast(UIntType) 1u << r) - 1,
815 upperMask = (~lowerMask) & this.max;
818 Collection of all state variables
819 used by the generator
821 private struct State
824 State array of the generator. This
825 is iterated through backwards (from
826 last element to first), providing a
827 few extra compiler optimizations by
828 comparison to the forward iteration
829 used in most implementations.
831 UIntType[n] data;
834 Cached copy of most recently updated
835 element of `data` state array, ready
836 to be tempered to generate next
837 `front` value
839 UIntType z;
842 Most recently generated random variate
844 UIntType front;
847 Index of the entry in the `data`
848 state array that will be twisted
849 in the next `popFront()` call
851 size_t index;
855 State variables used by the generator;
856 initialized to values equivalent to
857 explicitly seeding the generator with
858 `defaultSeed`
860 private State state = defaultState();
861 /* NOTE: the above is a workaround to ensure
862 backwards compatibility with the original
863 implementation, which permitted implicit
864 construction. With `@disable this();`
865 it would not be necessary. */
868 Constructs a MersenneTwisterEngine object.
870 this(UIntType value) @safe pure nothrow @nogc
872 seed(value);
876 Generates the default initial state for a Mersenne
877 Twister; equivalent to the internal state obtained
878 by calling `seed(defaultSeed)`
880 private static State defaultState() @safe pure nothrow @nogc
882 if (!__ctfe) assert(false);
883 State mtState;
884 seedImpl(defaultSeed, mtState);
885 return mtState;
889 Seeds a MersenneTwisterEngine object.
890 Note:
891 This seed function gives 2^w starting points (the lowest w bits of
892 the value provided will be used). To allow the RNG to be started
893 in any one of its internal states use the seed overload taking an
894 InputRange.
896 void seed()(UIntType value = defaultSeed) @safe pure nothrow @nogc
898 this.seedImpl(value, this.state);
902 Implementation of the seeding mechanism, which
903 can be used with an arbitrary `State` instance
905 private static void seedImpl(UIntType value, ref State mtState) @nogc
907 mtState.data[$ - 1] = value;
908 static if (this.max != UIntType.max)
910 mtState.data[$ - 1] &= this.max;
913 foreach_reverse (size_t i, ref e; mtState.data[0 .. $ - 1])
915 e = f * (mtState.data[i + 1] ^ (mtState.data[i + 1] >> (w - 2))) + cast(UIntType)(n - (i + 1));
916 static if (this.max != UIntType.max)
918 e &= this.max;
922 mtState.index = n - 1;
924 /* double popFront() to guarantee both `mtState.z`
925 and `mtState.front` are derived from the newly
926 set values in `mtState.data` */
927 MersenneTwisterEngine.popFrontImpl(mtState);
928 MersenneTwisterEngine.popFrontImpl(mtState);
932 Seeds a MersenneTwisterEngine object using an InputRange.
934 Throws:
935 `Exception` if the InputRange didn't provide enough elements to seed the generator.
936 The number of elements required is the 'n' template parameter of the MersenneTwisterEngine struct.
938 void seed(T)(T range) if (isInputRange!T && is(immutable ElementType!T == immutable UIntType))
940 this.seedImpl(range, this.state);
944 Implementation of the range-based seeding mechanism,
945 which can be used with an arbitrary `State` instance
947 private static void seedImpl(T)(T range, ref State mtState)
948 if (isInputRange!T && is(immutable ElementType!T == immutable UIntType))
950 size_t j;
951 for (j = 0; j < n && !range.empty; ++j, range.popFront())
953 ptrdiff_t idx = n - j - 1;
954 mtState.data[idx] = range.front;
957 mtState.index = n - 1;
959 if (range.empty && j < n)
961 import core.internal.string : UnsignedStringBuf, unsignedToTempString;
963 UnsignedStringBuf buf = void;
964 string s = "MersenneTwisterEngine.seed: Input range didn't provide enough elements: Need ";
965 s ~= unsignedToTempString(n, buf) ~ " elements.";
966 throw new Exception(s);
969 /* double popFront() to guarantee both `mtState.z`
970 and `mtState.front` are derived from the newly
971 set values in `mtState.data` */
972 MersenneTwisterEngine.popFrontImpl(mtState);
973 MersenneTwisterEngine.popFrontImpl(mtState);
977 Advances the generator.
979 void popFront() @safe pure nothrow @nogc
981 this.popFrontImpl(this.state);
985 Internal implementation of `popFront()`, which
986 can be used with an arbitrary `State` instance
988 private static void popFrontImpl(ref State mtState) @nogc
990 /* This function blends two nominally independent
991 processes: (i) calculation of the next random
992 variate `mtState.front` from the cached previous
993 `data` entry `z`, and (ii) updating the value
994 of `data[index]` and `mtState.z` and advancing
995 the `index` value to the next in sequence.
997 By interweaving the steps involved in these
998 procedures, rather than performing each of
999 them separately in sequence, the variables
1000 are kept 'hot' in CPU registers, allowing
1001 for significantly faster performance. */
1002 ptrdiff_t index = mtState.index;
1003 ptrdiff_t next = index - 1;
1004 if (next < 0)
1005 next = n - 1;
1006 auto z = mtState.z;
1007 ptrdiff_t conj = index - m;
1008 if (conj < 0)
1009 conj = index - m + n;
1011 static if (d == UIntType.max)
1013 z ^= (z >> u);
1015 else
1017 z ^= (z >> u) & d;
1020 auto q = mtState.data[index] & upperMask;
1021 auto p = mtState.data[next] & lowerMask;
1022 z ^= (z << s) & b;
1023 auto y = q | p;
1024 auto x = y >> 1;
1025 z ^= (z << t) & c;
1026 if (y & 1)
1027 x ^= a;
1028 auto e = mtState.data[conj] ^ x;
1029 z ^= (z >> l);
1030 mtState.z = mtState.data[index] = e;
1031 mtState.index = next;
1033 /* technically we should take the lowest `w`
1034 bits here, but if the tempering bitmasks
1035 `b` and `c` are set correctly, this should
1036 be unnecessary */
1037 mtState.front = z;
1041 Returns the current random value.
1043 @property UIntType front() @safe const pure nothrow @nogc
1045 return this.state.front;
1049 @property typeof(this) save() @safe const pure nothrow @nogc
1051 return this;
1055 Always `false`.
1057 enum bool empty = false;
1061 @safe unittest
1063 // seed with a constant
1064 Mt19937 gen;
1065 auto n = gen.front; // same for each run
1066 assert(n == 3499211612);
1068 // Seed with an unpredictable value
1069 gen.seed(unpredictableSeed);
1070 n = gen.front; // different across runs
1074 A `MersenneTwisterEngine` instantiated with the parameters of the
1075 original engine $(HTTP math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html,
1076 MT19937), generating uniformly-distributed 32-bit numbers with a
1077 period of 2 to the power of 19937. Recommended for random number
1078 generation unless memory is severely restricted, in which case a $(LREF
1079 LinearCongruentialEngine) would be the generator of choice.
1081 alias Mt19937 = MersenneTwisterEngine!(uint, 32, 624, 397, 31,
1082 0x9908b0df, 11, 0xffffffff, 7,
1083 0x9d2c5680, 15,
1084 0xefc60000, 18, 1_812_433_253);
1087 @safe @nogc unittest
1089 // seed with a constant
1090 Mt19937 gen;
1091 auto n = gen.front; // same for each run
1092 assert(n == 3499211612);
1094 // Seed with an unpredictable value
1095 gen.seed(unpredictableSeed);
1096 n = gen.front; // different across runs
1099 @safe nothrow unittest
1101 import std.algorithm;
1102 import std.range;
1103 static assert(isUniformRNG!Mt19937);
1104 static assert(isUniformRNG!(Mt19937, uint));
1105 static assert(isSeedable!Mt19937);
1106 static assert(isSeedable!(Mt19937, uint));
1107 static assert(isSeedable!(Mt19937, typeof(map!((a) => unpredictableSeed)(repeat(0)))));
1108 Mt19937 gen;
1109 assert(gen.front == 3499211612);
1110 popFrontN(gen, 9999);
1111 assert(gen.front == 4123659995);
1112 try { gen.seed(iota(624u)); } catch (Exception) { assert(false); }
1113 assert(gen.front == 3708921088u);
1114 popFrontN(gen, 9999);
1115 assert(gen.front == 165737292u);
1119 A `MersenneTwisterEngine` instantiated with the parameters of the
1120 original engine $(HTTP en.wikipedia.org/wiki/Mersenne_Twister,
1121 MT19937-64), generating uniformly-distributed 64-bit numbers with a
1122 period of 2 to the power of 19937.
1124 alias Mt19937_64 = MersenneTwisterEngine!(ulong, 64, 312, 156, 31,
1125 0xb5026f5aa96619e9, 29, 0x5555555555555555, 17,
1126 0x71d67fffeda60000, 37,
1127 0xfff7eee000000000, 43, 6_364_136_223_846_793_005);
1130 @safe @nogc unittest
1132 // Seed with a constant
1133 auto gen = Mt19937_64(12345);
1134 auto n = gen.front; // same for each run
1135 assert(n == 6597103971274460346);
1137 // Seed with an unpredictable value
1138 gen.seed(unpredictableSeed!ulong);
1139 n = gen.front; // different across runs
1142 @safe nothrow unittest
1144 import std.algorithm;
1145 import std.range;
1146 static assert(isUniformRNG!Mt19937_64);
1147 static assert(isUniformRNG!(Mt19937_64, ulong));
1148 static assert(isSeedable!Mt19937_64);
1149 static assert(isSeedable!(Mt19937_64, ulong));
1150 static assert(isSeedable!(Mt19937_64, typeof(map!((a) => unpredictableSeed!ulong)(repeat(0)))));
1151 Mt19937_64 gen;
1152 assert(gen.front == 14514284786278117030uL);
1153 popFrontN(gen, 9999);
1154 assert(gen.front == 9981545732273789042uL);
1155 try { gen.seed(iota(312uL)); } catch (Exception) { assert(false); }
1156 assert(gen.front == 14660652410669508483uL);
1157 popFrontN(gen, 9999);
1158 assert(gen.front == 15956361063660440239uL);
1161 @safe unittest
1163 import std.algorithm;
1164 import std.exception;
1165 import std.range;
1167 Mt19937 gen;
1169 assertThrown(gen.seed(map!((a) => 123_456_789U)(repeat(0, 623))));
1171 gen.seed(123_456_789U.repeat(624));
1172 //infinite Range
1173 gen.seed(123_456_789U.repeat);
1176 @safe @nogc pure nothrow unittest
1178 uint a, b;
1180 Mt19937 gen;
1181 a = gen.front;
1184 Mt19937 gen;
1185 gen.popFront();
1186 //popFrontN(gen, 1); // skip 1 element
1187 b = gen.front;
1189 assert(a != b);
1192 @safe @nogc unittest
1194 // Check .save works
1195 static foreach (Type; std.meta.AliasSeq!(Mt19937, Mt19937_64))
1197 auto gen1 = Type(123_456_789);
1198 gen1.popFront();
1199 // https://issues.dlang.org/show_bug.cgi?id=15853
1200 auto gen2 = ((const ref Type a) => a.save())(gen1);
1201 assert(gen1 == gen2); // Danger, Will Robinson -- no opEquals for MT
1202 // Enable next test when RNGs are reference types
1203 version (none) { assert(gen1 !is gen2); }
1204 for (auto i = 0; i < 100; i++, gen1.popFront, gen2.popFront)
1205 assert(gen1.front() == gen2.front());
1209 // https://issues.dlang.org/show_bug.cgi?id=11690
1210 @safe @nogc pure nothrow unittest
1212 alias MT(UIntType, uint w) = MersenneTwisterEngine!(UIntType, w, 624, 397, 31,
1213 0x9908b0df, 11, 0xffffffff, 7,
1214 0x9d2c5680, 15,
1215 0xefc60000, 18, 1812433253);
1217 static immutable ulong[] expectedFirstValue = [3499211612uL, 3499211612uL,
1218 171143175841277uL, 1145028863177033374uL];
1220 static immutable ulong[] expected10kValue = [4123659995uL, 4123659995uL,
1221 51991688252792uL, 3031481165133029945uL];
1223 static foreach (i, R; std.meta.AliasSeq!(MT!(uint, 32), MT!(ulong, 32), MT!(ulong, 48), MT!(ulong, 64)))
1225 auto a = R();
1226 a.seed(a.defaultSeed); // checks that some alternative paths in `seed` are utilized
1227 assert(a.front == expectedFirstValue[i]);
1228 a.popFrontN(9999);
1229 assert(a.front == expected10kValue[i]);
1234 Xorshift generator.
1235 Implemented according to $(HTTP www.jstatsoft.org/v08/i14/paper, Xorshift RNGs)
1236 (Marsaglia, 2003) when the size is small. For larger sizes the generator
1237 uses Sebastino Vigna's optimization of using an index to avoid needing
1238 to rotate the internal array.
1240 Period is `2 ^^ nbits - 1` except for a legacy 192-bit uint version (see
1241 note below).
1243 Params:
1244 UIntType = Word size of this xorshift generator and the return type
1245 of `opCall`.
1246 nbits = The number of bits of state of this generator. This must be
1247 a positive multiple of the size in bits of UIntType. If
1248 nbits is large this struct may occupy slightly more memory
1249 than this so it can use a circular counter instead of
1250 shifting the entire array.
1251 sa = The direction and magnitude of the 1st shift. Positive
1252 means left, negative means right.
1253 sb = The direction and magnitude of the 2nd shift. Positive
1254 means left, negative means right.
1255 sc = The direction and magnitude of the 3rd shift. Positive
1256 means left, negative means right.
1258 Note:
1259 For historical compatibility when `nbits == 192` and `UIntType` is `uint`
1260 a legacy hybrid PRNG is used consisting of a 160-bit xorshift combined
1261 with a 32-bit counter. This combined generator has period equal to the
1262 least common multiple of `2^^160 - 1` and `2^^32`.
1264 Previous versions of `XorshiftEngine` did not provide any mechanism to specify
1265 the directions of the shifts, taking each shift as an unsigned magnitude.
1266 For backwards compatibility, because three shifts in the same direction
1267 cannot result in a full-period XorshiftEngine, when all three of `sa`, `sb`,
1268 `sc, are positive `XorshiftEngine` treats them as unsigned magnitudes and
1269 uses shift directions to match the old behavior of `XorshiftEngine`.
1271 Not every set of shifts results in a full-period xorshift generator.
1272 The template does not currently at compile-time perform a full check
1273 for maximum period but in a future version might reject parameters
1274 resulting in shorter periods.
1276 struct XorshiftEngine(UIntType, uint nbits, int sa, int sb, int sc)
1277 if (isUnsigned!UIntType && !(sa > 0 && sb > 0 && sc > 0))
1279 static assert(nbits > 0 && nbits % (UIntType.sizeof * 8) == 0,
1280 "nbits must be an even multiple of "~UIntType.stringof
1281 ~".sizeof * 8, not "~nbits.stringof~".");
1283 static assert(!((sa >= 0) == (sb >= 0) && (sa >= 0) == (sc >= 0))
1284 && (sa * sb * sc != 0),
1285 "shifts cannot be zero and cannot all be in same direction: cannot be ["
1286 ~sa.stringof~", "~sb.stringof~", "~sc.stringof~"].");
1288 static assert(sa != sb && sb != sc,
1289 "consecutive shifts with the same magnitude and direction would partially or completely cancel!");
1291 static assert(UIntType.sizeof == uint.sizeof || UIntType.sizeof == ulong.sizeof,
1292 "XorshiftEngine currently does not support type " ~ UIntType.sizeof
1293 ~ " because it does not have a `seed` implementation for arrays "
1294 ~ " of element types other than uint and ulong.");
1296 public:
1297 ///Mark this as a Rng
1298 enum bool isUniformRandom = true;
1299 /// Always `false` (random generators are infinite ranges).
1300 enum empty = false;
1301 /// Smallest generated value.
1302 enum UIntType min = _state.length == 1 ? 1 : 0;
1303 /// Largest generated value.
1304 enum UIntType max = UIntType.max;
1307 private:
1308 // Legacy 192 bit uint hybrid counter/xorshift.
1309 enum bool isLegacy192Bit = UIntType.sizeof == uint.sizeof && nbits == 192;
1311 // Shift magnitudes.
1312 enum a = (sa < 0 ? -sa : sa);
1313 enum b = (sb < 0 ? -sb : sb);
1314 enum c = (sc < 0 ? -sc : sc);
1316 // Shift expressions to mix in.
1317 enum shiftA(string expr) = `((`~expr~`) `~(sa > 0 ? `<< a)` : ` >>> a)`);
1318 enum shiftB(string expr) = `((`~expr~`) `~(sb > 0 ? `<< b)` : ` >>> b)`);
1319 enum shiftC(string expr) = `((`~expr~`) `~(sc > 0 ? `<< c)` : ` >>> c)`);
1321 enum N = nbits / (UIntType.sizeof * 8);
1323 // For N <= 2 it is strictly worse to use an index.
1324 // Informal third-party benchmarks suggest that for `ulong` it is
1325 // faster to use an index when N == 4. For `uint` we err on the side
1326 // of not increasing the struct's size and only switch to the other
1327 // implementation when N > 4.
1328 enum useIndex = !isLegacy192Bit && (UIntType.sizeof >= ulong.sizeof ? N > 3 : N > 4);
1329 static if (useIndex)
1331 enum initialIndex = N - 1;
1332 uint _index = initialIndex;
1335 static if (N == 1 && UIntType.sizeof <= uint.sizeof)
1337 UIntType[N] _state = [cast(UIntType) 2_463_534_242];
1339 else static if (isLegacy192Bit)
1341 UIntType[N] _state = [123_456_789, 362_436_069, 521_288_629, 88_675_123, 5_783_321, 6_615_241];
1342 UIntType value_;
1344 else static if (N <= 5 && UIntType.sizeof <= uint.sizeof)
1346 UIntType[N] _state = [
1347 cast(UIntType) 123_456_789,
1348 cast(UIntType) 362_436_069,
1349 cast(UIntType) 521_288_629,
1350 cast(UIntType) 88_675_123,
1351 cast(UIntType) 5_783_321][0 .. N];
1353 else
1355 UIntType[N] _state = ()
1357 static if (UIntType.sizeof < ulong.sizeof)
1359 uint x0 = 123_456_789;
1360 enum uint m = 1_812_433_253U;
1362 else static if (UIntType.sizeof <= ulong.sizeof)
1364 ulong x0 = 123_456_789;
1365 enum ulong m = 6_364_136_223_846_793_005UL;
1367 else
1369 static assert(0, "Phobos Error: Xorshift has no instantiation rule for "
1370 ~ UIntType.stringof);
1372 enum uint rshift = typeof(x0).sizeof * 8 - 2;
1373 UIntType[N] result = void;
1374 foreach (i, ref e; result)
1376 e = cast(UIntType) (x0 = (m * (x0 ^ (x0 >>> rshift)) + i + 1));
1377 if (e == 0)
1378 e = cast(UIntType) (i + 1);
1380 return result;
1381 }();
1385 public:
1387 Constructs a `XorshiftEngine` generator seeded with $(D_PARAM x0).
1389 Params:
1390 x0 = value used to deterministically initialize internal state
1392 this()(UIntType x0) @safe pure nothrow @nogc
1394 seed(x0);
1399 (Re)seeds the generator.
1401 Params:
1402 x0 = value used to deterministically initialize internal state
1404 void seed()(UIntType x0) @safe pure nothrow @nogc
1406 static if (useIndex)
1407 _index = initialIndex;
1409 static if (UIntType.sizeof == uint.sizeof)
1411 // Initialization routine from MersenneTwisterEngine.
1412 foreach (uint i, ref e; _state)
1414 e = (x0 = (1_812_433_253U * (x0 ^ (x0 >> 30)) + i + 1));
1415 // Xorshift requires merely that not every word of the internal
1416 // array is 0. For historical compatibility the 32-bit word version
1417 // has the stronger requirement that not any word of the state
1418 // array is 0 after initial seeding.
1419 if (e == 0)
1420 e = (i + 1);
1423 else static if (UIntType.sizeof == ulong.sizeof)
1425 static if (N > 1)
1427 // Initialize array using splitmix64 as recommended by Sebastino Vigna.
1428 // By construction the array will not be all zeroes.
1429 // http://xoroshiro.di.unimi.it/splitmix64.c
1430 foreach (ref e; _state)
1432 ulong z = (x0 += 0x9e37_79b9_7f4a_7c15UL);
1433 z = (z ^ (z >>> 30)) * 0xbf58_476d_1ce4_e5b9UL;
1434 z = (z ^ (z >>> 27)) * 0x94d0_49bb_1331_11ebUL;
1435 e = z ^ (z >>> 31);
1438 else
1440 // Apply a transformation when N == 1 instead of just copying x0
1441 // directly because it's not unlikely that a user might initialize
1442 // a PRNG with small counting numbers (e.g. 1, 2, 3) that have the
1443 // statistically rare property of having only 1 or 2 non-zero bits.
1444 // Additionally we need to ensure that the internal state is not
1445 // entirely zero.
1446 if (x0 != 0)
1447 _state[0] = x0 * 6_364_136_223_846_793_005UL;
1448 else
1449 _state[0] = typeof(this).init._state[0];
1452 else
1454 static assert(0, "Phobos Error: Xorshift has no `seed` implementation for "
1455 ~ UIntType.stringof);
1458 popFront();
1463 * Returns the current number in the random sequence.
1465 @property
1466 UIntType front() const @safe pure nothrow @nogc
1468 static if (isLegacy192Bit)
1469 return value_;
1470 else static if (!useIndex)
1471 return _state[N-1];
1472 else
1473 return _state[_index];
1478 * Advances the random sequence.
1480 void popFront() @safe pure nothrow @nogc
1482 alias s = _state;
1483 static if (isLegacy192Bit)
1485 auto x = _state[0] ^ mixin(shiftA!`s[0]`);
1486 static foreach (i; 0 .. N-2)
1487 s[i] = s[i + 1];
1488 s[N-2] = s[N-2] ^ mixin(shiftC!`s[N-2]`) ^ x ^ mixin(shiftB!`x`);
1489 value_ = s[N-2] + (s[N-1] += 362_437);
1491 else static if (N == 1)
1493 s[0] ^= mixin(shiftA!`s[0]`);
1494 s[0] ^= mixin(shiftB!`s[0]`);
1495 s[0] ^= mixin(shiftC!`s[0]`);
1497 else static if (!useIndex)
1499 auto x = s[0] ^ mixin(shiftA!`s[0]`);
1500 static foreach (i; 0 .. N-1)
1501 s[i] = s[i + 1];
1502 s[N-1] = s[N-1] ^ mixin(shiftC!`s[N-1]`) ^ x ^ mixin(shiftB!`x`);
1504 else
1506 assert(_index < N); // Invariant.
1507 const sIndexMinus1 = s[_index];
1508 static if ((N & (N - 1)) == 0)
1509 _index = (_index + 1) & typeof(_index)(N - 1);
1510 else
1512 if (++_index >= N)
1513 _index = 0;
1515 auto x = s[_index];
1516 x ^= mixin(shiftA!`x`);
1517 s[_index] = sIndexMinus1 ^ mixin(shiftC!`sIndexMinus1`) ^ x ^ mixin(shiftB!`x`);
1523 * Captures a range state.
1525 @property
1526 typeof(this) save() const @safe pure nothrow @nogc
1528 return this;
1531 private:
1532 // Workaround for a DScanner bug. If we remove this `private:` DScanner
1533 // gives erroneous warnings about missing documentation for public symbols
1534 // later in the module.
1537 /// ditto
1538 template XorshiftEngine(UIntType, int bits, int a, int b, int c)
1539 if (isUnsigned!UIntType && a > 0 && b > 0 && c > 0)
1541 // Compatibility with old parameterizations without explicit shift directions.
1542 static if (bits == UIntType.sizeof * 8)
1543 alias XorshiftEngine = .XorshiftEngine!(UIntType, bits, a, -b, c);//left, right, left
1544 else static if (bits == 192 && UIntType.sizeof == uint.sizeof)
1545 alias XorshiftEngine = .XorshiftEngine!(UIntType, bits, -a, b, c);//right, left, left
1546 else
1547 alias XorshiftEngine = .XorshiftEngine!(UIntType, bits, a, -b, -c);//left, right, right
1551 @safe unittest
1553 alias Xorshift96 = XorshiftEngine!(uint, 96, 10, 5, 26);
1554 auto rnd = Xorshift96(42);
1555 auto num = rnd.front; // same for each run
1556 assert(num == 2704588748);
1561 * Define `XorshiftEngine` generators with well-chosen parameters. See each bits examples of "Xorshift RNGs".
1562 * `Xorshift` is a Xorshift128's alias because 128bits implementation is mostly used.
1564 alias Xorshift32 = XorshiftEngine!(uint, 32, 13, 17, 15) ;
1565 alias Xorshift64 = XorshiftEngine!(uint, 64, 10, 13, 10); /// ditto
1566 alias Xorshift96 = XorshiftEngine!(uint, 96, 10, 5, 26); /// ditto
1567 alias Xorshift128 = XorshiftEngine!(uint, 128, 11, 8, 19); /// ditto
1568 alias Xorshift160 = XorshiftEngine!(uint, 160, 2, 1, 4); /// ditto
1569 alias Xorshift192 = XorshiftEngine!(uint, 192, 2, 1, 4); /// ditto
1570 alias Xorshift = Xorshift128; /// ditto
1573 @safe @nogc unittest
1575 // Seed with a constant
1576 auto rnd = Xorshift(1);
1577 auto num = rnd.front; // same for each run
1578 assert(num == 1405313047);
1580 // Seed with an unpredictable value
1581 rnd.seed(unpredictableSeed);
1582 num = rnd.front; // different across rnd
1585 @safe @nogc unittest
1587 import std.range;
1588 static assert(isForwardRange!Xorshift);
1589 static assert(isUniformRNG!Xorshift);
1590 static assert(isUniformRNG!(Xorshift, uint));
1591 static assert(isSeedable!Xorshift);
1592 static assert(isSeedable!(Xorshift, uint));
1594 static assert(Xorshift32.min == 1);
1596 // Result from reference implementation.
1597 static ulong[][] checking = [
1598 [2463534242UL, 901999875, 3371835698, 2675058524, 1053936272, 3811264849,
1599 472493137, 3856898176, 2131710969, 2312157505],
1600 [362436069UL, 2113136921, 19051112, 3010520417, 951284840, 1213972223,
1601 3173832558, 2611145638, 2515869689, 2245824891],
1602 [521288629UL, 1950277231, 185954712, 1582725458, 3580567609, 2303633688,
1603 2394948066, 4108622809, 1116800180, 3357585673],
1604 [88675123UL, 3701687786, 458299110, 2500872618, 3633119408, 516391518,
1605 2377269574, 2599949379, 717229868, 137866584],
1606 [5783321UL, 393427209, 1947109840, 565829276, 1006220149, 971147905,
1607 1436324242, 2800460115, 1484058076, 3823330032],
1608 [0UL, 246875399, 3690007200, 1264581005, 3906711041, 1866187943, 2481925219,
1609 2464530826, 1604040631, 3653403911],
1610 [16749904790159980466UL, 14489774923612894650UL, 148813570191443371UL,
1611 6529783008780612886UL, 10182425759614080046UL, 16549659571055687844UL,
1612 542957868271744939UL, 9459127085596028450UL, 16001049981702441780UL,
1613 7351634712593111741],
1614 [14750058843113249055UL, 17731577314455387619UL, 1314705253499959044UL,
1615 3113030620614841056UL, 9468075444678629182UL, 13962152036600088141UL,
1616 9030205609946043947UL, 1856726150434672917UL, 8098922200110395314UL,
1617 2772699174618556175UL],
1620 alias XorshiftTypes = std.meta.AliasSeq!(Xorshift32, Xorshift64, Xorshift96,
1621 Xorshift128, Xorshift160, Xorshift192, Xorshift64_64, Xorshift128_64);
1623 foreach (I, Type; XorshiftTypes)
1625 Type rnd;
1627 foreach (e; checking[I])
1629 assert(rnd.front == e);
1630 rnd.popFront();
1634 // Check .save works
1635 foreach (Type; XorshiftTypes)
1637 auto rnd1 = Type(123_456_789);
1638 rnd1.popFront();
1639 // https://issues.dlang.org/show_bug.cgi?id=15853
1640 auto rnd2 = ((const ref Type a) => a.save())(rnd1);
1641 assert(rnd1 == rnd2);
1642 // Enable next test when RNGs are reference types
1643 version (none) { assert(rnd1 !is rnd2); }
1644 for (auto i = 0; i <= Type.sizeof / 4; i++, rnd1.popFront, rnd2.popFront)
1645 assert(rnd1.front() == rnd2.front());
1650 /* A complete list of all pseudo-random number generators implemented in
1651 * std.random. This can be used to confirm that a given function or
1652 * object is compatible with all the pseudo-random number generators
1653 * available. It is enabled only in unittest mode.
1655 @safe @nogc unittest
1657 foreach (Rng; PseudoRngTypes)
1659 static assert(isUniformRNG!Rng);
1660 auto rng = Rng(123_456_789);
1664 version (CRuntime_Bionic)
1665 version = SecureARC4Random; // ChaCha20
1666 version (Darwin)
1667 version = SecureARC4Random; // AES
1668 version (OpenBSD)
1669 version = SecureARC4Random; // ChaCha20
1670 version (NetBSD)
1671 version = SecureARC4Random; // ChaCha20
1673 version (CRuntime_UClibc)
1674 version = LegacyARC4Random; // ARC4
1675 version (FreeBSD)
1676 version = LegacyARC4Random; // ARC4
1677 version (DragonFlyBSD)
1678 version = LegacyARC4Random; // ARC4
1679 version (BSD)
1680 version = LegacyARC4Random; // Unknown implementation
1682 // For the current purpose of unpredictableSeed the difference between
1683 // a secure arc4random implementation and a legacy implementation is
1684 // unimportant. The source code documents this distinction in case in the
1685 // future Phobos is altered to require cryptographically secure sources
1686 // of randomness, and also so other people reading this source code (as
1687 // Phobos is often looked to as an example of good D programming practices)
1688 // do not mistakenly use insecure versions of arc4random in contexts where
1689 // cryptographically secure sources of randomness are needed.
1691 // Performance note: ChaCha20 is about 70% faster than ARC4, contrary to
1692 // what one might assume from it being more secure.
1694 version (SecureARC4Random)
1695 version = AnyARC4Random;
1696 version (LegacyARC4Random)
1697 version = AnyARC4Random;
1699 version (AnyARC4Random)
1701 extern(C) private @nogc nothrow
1703 uint arc4random() @safe;
1704 void arc4random_buf(scope void* buf, size_t nbytes) @system;
1707 else
1709 private ulong bootstrapSeed() @nogc nothrow
1711 // https://issues.dlang.org/show_bug.cgi?id=19580
1712 // previously used `ulong result = void` to start with an arbitary value
1713 // but using an uninitialized variable's value is undefined behavior
1714 // and enabled unwanted optimizations on non-DMD compilers.
1715 ulong result;
1716 enum ulong m = 0xc6a4_a793_5bd1_e995UL; // MurmurHash2_64A constant.
1717 void updateResult(ulong x)
1719 x *= m;
1720 x = (x ^ (x >>> 47)) * m;
1721 result = (result ^ x) * m;
1723 import core.thread : getpid, Thread;
1724 import core.time : MonoTime;
1726 updateResult(cast(ulong) cast(void*) Thread.getThis());
1727 updateResult(cast(ulong) getpid());
1728 updateResult(cast(ulong) MonoTime.currTime.ticks);
1729 result = (result ^ (result >>> 47)) * m;
1730 return result ^ (result >>> 47);
1733 // If we don't have arc4random and we don't have RDRAND fall back to this.
1734 private ulong fallbackSeed() @nogc nothrow
1736 // Bit avalanche function from MurmurHash3.
1737 static ulong fmix64(ulong k) @nogc nothrow pure @safe
1739 k = (k ^ (k >>> 33)) * 0xff51afd7ed558ccd;
1740 k = (k ^ (k >>> 33)) * 0xc4ceb9fe1a85ec53;
1741 return k ^ (k >>> 33);
1743 // Using SplitMix algorithm with constant gamma.
1744 // Chosen gamma is the odd number closest to 2^^64
1745 // divided by the silver ratio (1.0L + sqrt(2.0L)).
1746 enum gamma = 0x6a09e667f3bcc909UL;
1747 import core.atomic : has64BitCAS;
1748 static if (has64BitCAS)
1750 import core.atomic : MemoryOrder, atomicLoad, atomicOp, atomicStore, cas;
1751 shared static ulong seed;
1752 shared static bool initialized;
1753 if (0 == atomicLoad!(MemoryOrder.raw)(initialized))
1755 cas(&seed, 0UL, fmix64(bootstrapSeed()));
1756 atomicStore!(MemoryOrder.rel)(initialized, true);
1758 return fmix64(atomicOp!"+="(seed, gamma));
1760 else
1762 static ulong seed;
1763 static bool initialized;
1764 if (!initialized)
1766 seed = fmix64(bootstrapSeed());
1767 initialized = true;
1769 return fmix64(seed += gamma);
1775 A "good" seed for initializing random number engines. Initializing
1776 with $(D_PARAM unpredictableSeed) makes engines generate different
1777 random number sequences every run.
1779 Returns:
1780 A single unsigned integer seed value, different on each successive call
1781 Note:
1782 In general periodically 'reseeding' a PRNG does not improve its quality
1783 and in some cases may harm it. For an extreme example the Mersenne
1784 Twister has `2 ^^ 19937 - 1` distinct states but after `seed(uint)` is
1785 called it can only be in one of `2 ^^ 32` distinct states regardless of
1786 how excellent the source of entropy is.
1788 @property uint unpredictableSeed() @trusted nothrow @nogc
1790 version (AnyARC4Random)
1792 return arc4random();
1794 else
1796 version (InlineAsm_X86_Any)
1798 import core.cpuid : hasRdrand;
1799 if (hasRdrand)
1801 uint result;
1802 asm @nogc nothrow
1804 db 0x0f, 0xc7, 0xf0; // rdrand EAX
1805 jnc LnotUsingRdrand;
1806 mov result, EAX;
1807 // Some AMD CPUs shipped with bugs where RDRAND could fail
1808 // but still set the carry flag to 1. In those cases the
1809 // output will be -1.
1810 cmp EAX, 0xffff_ffff;
1811 jne LusingRdrand;
1812 // If result was -1 verify RDAND isn't constantly returning -1.
1813 db 0x0f, 0xc7, 0xf0;
1814 jnc LusingRdrand;
1815 cmp EAX, 0xffff_ffff;
1816 je LnotUsingRdrand;
1818 LusingRdrand:
1819 return result;
1821 LnotUsingRdrand:
1823 return cast(uint) fallbackSeed();
1827 /// ditto
1828 template unpredictableSeed(UIntType)
1829 if (isUnsigned!UIntType)
1831 static if (is(UIntType == uint))
1832 alias unpredictableSeed = .unpredictableSeed;
1833 else static if (!is(Unqual!UIntType == UIntType))
1834 alias unpredictableSeed = .unpredictableSeed!(Unqual!UIntType);
1835 else
1836 /// ditto
1837 @property UIntType unpredictableSeed() @nogc nothrow @trusted
1839 version (AnyARC4Random)
1841 static if (UIntType.sizeof <= uint.sizeof)
1843 return cast(UIntType) arc4random();
1845 else
1847 UIntType result = void;
1848 arc4random_buf(&result, UIntType.sizeof);
1849 return result;
1852 else
1854 version (InlineAsm_X86_Any)
1856 import core.cpuid : hasRdrand;
1857 if (hasRdrand)
1859 static if (UIntType.sizeof <= uint.sizeof)
1861 uint result;
1862 asm @nogc nothrow
1864 db 0x0f, 0xc7, 0xf0; // rdrand EAX
1865 jnc LnotUsingRdrand;
1866 mov result, EAX;
1867 // Some AMD CPUs shipped with bugs where RDRAND could fail
1868 // but still set the carry flag to 1. In those cases the
1869 // output will be -1.
1870 cmp EAX, 0xffff_ffff;
1871 jne LusingRdrand;
1872 // If result was -1 verify RDAND isn't constantly returning -1.
1873 db 0x0f, 0xc7, 0xf0;
1874 jnc LusingRdrand;
1875 cmp EAX, 0xffff_ffff;
1876 je LnotUsingRdrand;
1878 LusingRdrand:
1879 return cast(UIntType) result;
1881 else version (D_InlineAsm_X86_64)
1883 ulong result;
1884 asm @nogc nothrow
1886 db 0x48, 0x0f, 0xc7, 0xf0; // rdrand RAX
1887 jnc LnotUsingRdrand;
1888 mov result, RAX;
1889 // Some AMD CPUs shipped with bugs where RDRAND could fail
1890 // but still set the carry flag to 1. In those cases the
1891 // output will be -1.
1892 cmp RAX, 0xffff_ffff_ffff_ffff;
1893 jne LusingRdrand;
1894 // If result was -1 verify RDAND isn't constantly returning -1.
1895 db 0x48, 0x0f, 0xc7, 0xf0;
1896 jnc LusingRdrand;
1897 cmp RAX, 0xffff_ffff_ffff_ffff;
1898 je LnotUsingRdrand;
1900 LusingRdrand:
1901 return result;
1903 else
1905 uint resultLow, resultHigh;
1906 asm @nogc nothrow
1908 db 0x0f, 0xc7, 0xf0; // rdrand EAX
1909 jnc LnotUsingRdrand;
1910 mov resultLow, EAX;
1911 db 0x0f, 0xc7, 0xf0; // rdrand EAX
1912 jnc LnotUsingRdrand;
1913 mov resultHigh, EAX;
1915 if (resultLow != uint.max || resultHigh != uint.max) // Protect against AMD RDRAND bug.
1916 return ((cast(ulong) resultHigh) << 32) ^ resultLow;
1919 LnotUsingRdrand:
1921 return cast(UIntType) fallbackSeed();
1927 @safe @nogc unittest
1929 auto rnd = Random(unpredictableSeed);
1930 auto n = rnd.front;
1931 static assert(is(typeof(n) == uint));
1935 The "default", "favorite", "suggested" random number generator type on
1936 the current platform. It is an alias for one of the previously-defined
1937 generators. You may want to use it if (1) you need to generate some
1938 nice random numbers, and (2) you don't care for the minutiae of the
1939 method being used.
1942 alias Random = Mt19937;
1944 @safe @nogc unittest
1946 static assert(isUniformRNG!Random);
1947 static assert(isUniformRNG!(Random, uint));
1948 static assert(isSeedable!Random);
1949 static assert(isSeedable!(Random, uint));
1953 Global random number generator used by various functions in this
1954 module whenever no generator is specified. It is allocated per-thread
1955 and initialized to an unpredictable value for each thread.
1957 Returns:
1958 A singleton instance of the default random number generator
1960 @property ref Random rndGen() @safe nothrow @nogc
1962 static Random result;
1963 static bool initialized;
1964 if (!initialized)
1966 static if (isSeedable!(Random, ulong))
1967 result.seed(unpredictableSeed!ulong); // Avoid unnecessary copy.
1968 else static if (is(Random : MersenneTwisterEngine!Params, Params...))
1969 initMTEngine(result);
1970 else static if (isSeedable!(Random, uint))
1971 result.seed(unpredictableSeed!uint); // Avoid unnecessary copy.
1972 else
1973 result = Random(unpredictableSeed);
1974 initialized = true;
1976 return result;
1980 @safe nothrow @nogc unittest
1982 import std.algorithm.iteration : sum;
1983 import std.range : take;
1984 auto rnd = rndGen;
1985 assert(rnd.take(3).sum > 0);
1989 Initialize a 32-bit MersenneTwisterEngine from 64 bits of entropy.
1990 This is private and accepts no seed as a parameter, freeing the internal
1991 implementaton from any need for stability across releases.
1993 private void initMTEngine(MTEngine)(scope ref MTEngine mt)
1994 if (is(MTEngine : MersenneTwisterEngine!Params, Params...))
1996 pragma(inline, false); // Called no more than once per thread by rndGen.
1997 ulong seed = unpredictableSeed!ulong;
1998 static if (is(typeof(mt.seed(seed))))
2000 mt.seed(seed);
2002 else
2004 alias UIntType = typeof(mt.front());
2005 if (seed == 0) seed = -1; // Any number but 0 is fine.
2006 uint s0 = cast(uint) seed;
2007 uint s1 = cast(uint) (seed >> 32);
2008 foreach (ref e; mt.state.data)
2010 //http://xoshiro.di.unimi.it/xoroshiro64starstar.c
2011 const tmp = s0 * 0x9E3779BB;
2012 e = ((tmp << 5) | (tmp >> (32 - 5))) * 5;
2013 static if (MTEngine.max != UIntType.max) { e &= MTEngine.max; }
2015 const tmp1 = s0 ^ s1;
2016 s0 = ((s0 << 26) | (s0 >> (32 - 26))) ^ tmp1 ^ (tmp1 << 9);
2017 s1 = (tmp1 << 13) | (tmp1 >> (32 - 13));
2020 mt.state.index = mt.state.data.length - 1;
2021 // double popFront() to guarantee both `mt.state.z`
2022 // and `mt.state.front` are derived from the newly
2023 // set values in `mt.state.data`.
2024 mt.popFront();
2025 mt.popFront();
2030 Generates a number between `a` and `b`. The `boundaries`
2031 parameter controls the shape of the interval (open vs. closed on
2032 either side). Valid values for `boundaries` are `"[]"`, $(D
2033 "$(LPAREN)]"), `"[$(RPAREN)"`, and `"()"`. The default interval
2034 is closed to the left and open to the right. The version that does not
2035 take `urng` uses the default generator `rndGen`.
2037 Params:
2038 a = lower bound of the _uniform distribution
2039 b = upper bound of the _uniform distribution
2040 urng = (optional) random number generator to use;
2041 if not specified, defaults to `rndGen`
2043 Returns:
2044 A single random variate drawn from the _uniform distribution
2045 between `a` and `b`, whose type is the common type of
2046 these parameters
2048 auto uniform(string boundaries = "[)", T1, T2)
2049 (T1 a, T2 b)
2050 if (!is(CommonType!(T1, T2) == void))
2052 return uniform!(boundaries, T1, T2, Random)(a, b, rndGen);
2056 @safe unittest
2058 auto rnd = Random(unpredictableSeed);
2060 // Generate an integer in [0, 1023]
2061 auto a = uniform(0, 1024, rnd);
2062 assert(0 <= a && a < 1024);
2064 // Generate a float in [0, 1)
2065 auto b = uniform(0.0f, 1.0f, rnd);
2066 assert(0 <= b && b < 1);
2068 // Generate a float in [0, 1]
2069 b = uniform!"[]"(0.0f, 1.0f, rnd);
2070 assert(0 <= b && b <= 1);
2072 // Generate a float in (0, 1)
2073 b = uniform!"()"(0.0f, 1.0f, rnd);
2074 assert(0 < b && b < 1);
2077 /// Create an array of random numbers using range functions and UFCS
2078 @safe unittest
2080 import std.array : array;
2081 import std.range : generate, takeExactly;
2083 int[] arr = generate!(() => uniform(0, 100)).takeExactly(10).array;
2084 assert(arr.length == 10);
2085 assert(arr[0] >= 0 && arr[0] < 100);
2088 @safe unittest
2090 MinstdRand0 gen;
2091 foreach (i; 0 .. 20)
2093 auto x = uniform(0.0, 15.0, gen);
2094 assert(0 <= x && x < 15);
2096 foreach (i; 0 .. 20)
2098 auto x = uniform!"[]"('a', 'z', gen);
2099 assert('a' <= x && x <= 'z');
2102 foreach (i; 0 .. 20)
2104 auto x = uniform('a', 'z', gen);
2105 assert('a' <= x && x < 'z');
2108 foreach (i; 0 .. 20)
2110 immutable ubyte a = 0;
2111 immutable ubyte b = 15;
2112 auto x = uniform(a, b, gen);
2113 assert(a <= x && x < b);
2117 // Implementation of uniform for floating-point types
2118 /// ditto
2119 auto uniform(string boundaries = "[)",
2120 T1, T2, UniformRandomNumberGenerator)
2121 (T1 a, T2 b, ref UniformRandomNumberGenerator urng)
2122 if (isFloatingPoint!(CommonType!(T1, T2)) && isUniformRNG!UniformRandomNumberGenerator)
2124 import std.conv : text;
2125 import std.exception : enforce;
2126 alias NumberType = Unqual!(CommonType!(T1, T2));
2127 static if (boundaries[0] == '(')
2129 import std.math.operations : nextafter;
2130 NumberType _a = nextafter(cast(NumberType) a, NumberType.infinity);
2132 else
2134 NumberType _a = a;
2136 static if (boundaries[1] == ')')
2138 import std.math.operations : nextafter;
2139 NumberType _b = nextafter(cast(NumberType) b, -NumberType.infinity);
2141 else
2143 NumberType _b = b;
2145 enforce(_a <= _b,
2146 text("std.random.uniform(): invalid bounding interval ",
2147 boundaries[0], a, ", ", b, boundaries[1]));
2148 NumberType result =
2149 _a + (_b - _a) * cast(NumberType) (urng.front - urng.min)
2150 / (urng.max - urng.min);
2151 urng.popFront();
2152 return result;
2155 // Implementation of uniform for integral types
2156 /+ Description of algorithm and suggestion of correctness:
2158 The modulus operator maps an integer to a small, finite space. For instance, `x
2159 % 3` will map whatever x is into the range [0 .. 3). 0 maps to 0, 1 maps to 1, 2
2160 maps to 2, 3 maps to 0, and so on infinitely. As long as the integer is
2161 uniformly chosen from the infinite space of all non-negative integers then `x %
2162 3` will uniformly fall into that range.
2164 (Non-negative is important in this case because some definitions of modulus,
2165 namely the one used in computers generally, map negative numbers differently to
2166 (-3 .. 0]. `uniform` does not use negative number modulus, thus we can safely
2167 ignore that fact.)
2169 The issue with computers is that integers have a finite space they must fit in,
2170 and our uniformly chosen random number is picked in that finite space. So, that
2171 method is not sufficient. You can look at it as the integer space being divided
2172 into "buckets" and every bucket after the first bucket maps directly into that
2173 first bucket. `[0, 1, 2]`, `[3, 4, 5]`, ... When integers are finite, then the
2174 last bucket has the chance to be "incomplete": `[uint.max - 3, uint.max - 2,
2175 uint.max - 1]`, `[uint.max]` ... (the last bucket only has 1!). The issue here
2176 is that _every_ bucket maps _completely_ to the first bucket except for that
2177 last one. The last one doesn't have corresponding mappings to 1 or 2, in this
2178 case, which makes it unfair.
2180 So, the answer is to simply "reroll" if you're in that last bucket, since it's
2181 the only unfair one. Eventually you'll roll into a fair bucket. Simply, instead
2182 of the meaning of the last bucket being "maps to `[0]`", it changes to "maps to
2183 `[0, 1, 2]`", which is precisely what we want.
2185 To generalize, `upperDist` represents the size of our buckets (and, thus, the
2186 exclusive upper bound for our desired uniform number). `rnum` is a uniformly
2187 random number picked from the space of integers that a computer can hold (we'll
2188 say `UpperType` represents that type).
2190 We'll first try to do the mapping into the first bucket by doing `offset = rnum
2191 % upperDist`. We can figure out the position of the front of the bucket we're in
2192 by `bucketFront = rnum - offset`.
2194 If we start at `UpperType.max` and walk backwards `upperDist - 1` spaces, then
2195 the space we land on is the last acceptable position where a full bucket can
2196 fit:
2199 bucketFront UpperType.max
2201 [..., 0, 1, 2, ..., upperDist - 1]
2202 ^~~ upperDist - 1 ~~^
2205 If the bucket starts any later, then it must have lost at least one number and
2206 at least that number won't be represented fairly.
2209 bucketFront UpperType.max
2211 [..., upperDist - 1, 0, 1, 2, ..., upperDist - 2]
2212 ^~~~~~~~ upperDist - 1 ~~~~~~~^
2215 Hence, our condition to reroll is
2216 `bucketFront > (UpperType.max - (upperDist - 1))`
2218 auto uniform(string boundaries = "[)", T1, T2, RandomGen)
2219 (T1 a, T2 b, ref RandomGen rng)
2220 if ((isIntegral!(CommonType!(T1, T2)) || isSomeChar!(CommonType!(T1, T2))) &&
2221 isUniformRNG!RandomGen)
2223 import std.conv : text, unsigned;
2224 import std.exception : enforce;
2225 alias ResultType = Unqual!(CommonType!(T1, T2));
2226 static if (boundaries[0] == '(')
2228 enforce(a < ResultType.max,
2229 text("std.random.uniform(): invalid left bound ", a));
2230 ResultType lower = cast(ResultType) (a + 1);
2232 else
2234 ResultType lower = a;
2237 static if (boundaries[1] == ']')
2239 enforce(lower <= b,
2240 text("std.random.uniform(): invalid bounding interval ",
2241 boundaries[0], a, ", ", b, boundaries[1]));
2242 /* Cannot use this next optimization with dchar, as dchar
2243 * only partially uses its full bit range
2245 static if (!is(ResultType == dchar))
2247 if (b == ResultType.max && lower == ResultType.min)
2249 // Special case - all bits are occupied
2250 return std.random.uniform!ResultType(rng);
2253 auto upperDist = unsigned(b - lower) + 1u;
2255 else
2257 enforce(lower < b,
2258 text("std.random.uniform(): invalid bounding interval ",
2259 boundaries[0], a, ", ", b, boundaries[1]));
2260 auto upperDist = unsigned(b - lower);
2263 assert(upperDist != 0);
2265 alias UpperType = typeof(upperDist);
2266 static assert(UpperType.min == 0);
2268 UpperType offset, rnum, bucketFront;
2271 rnum = uniform!UpperType(rng);
2272 offset = rnum % upperDist;
2273 bucketFront = rnum - offset;
2274 } // while we're in an unfair bucket...
2275 while (bucketFront > (UpperType.max - (upperDist - 1)));
2277 return cast(ResultType)(lower + offset);
2280 @safe unittest
2282 import std.conv : to;
2283 auto gen = Mt19937(123_456_789);
2284 static assert(isForwardRange!(typeof(gen)));
2286 auto a = uniform(0, 1024, gen);
2287 assert(0 <= a && a <= 1024);
2288 auto b = uniform(0.0f, 1.0f, gen);
2289 assert(0 <= b && b < 1, to!string(b));
2290 auto c = uniform(0.0, 1.0);
2291 assert(0 <= c && c < 1);
2293 static foreach (T; std.meta.AliasSeq!(char, wchar, dchar, byte, ubyte, short, ushort,
2294 int, uint, long, ulong, float, double, real))
2296 T lo = 0, hi = 100;
2298 // Try tests with each of the possible bounds
2300 T init = uniform(lo, hi);
2301 size_t i = 50;
2302 while (--i && uniform(lo, hi) == init) {}
2303 assert(i > 0);
2306 T init = uniform!"[)"(lo, hi);
2307 size_t i = 50;
2308 while (--i && uniform(lo, hi) == init) {}
2309 assert(i > 0);
2312 T init = uniform!"(]"(lo, hi);
2313 size_t i = 50;
2314 while (--i && uniform(lo, hi) == init) {}
2315 assert(i > 0);
2318 T init = uniform!"()"(lo, hi);
2319 size_t i = 50;
2320 while (--i && uniform(lo, hi) == init) {}
2321 assert(i > 0);
2324 T init = uniform!"[]"(lo, hi);
2325 size_t i = 50;
2326 while (--i && uniform(lo, hi) == init) {}
2327 assert(i > 0);
2330 /* Test case with closed boundaries covering whole range
2331 * of integral type
2333 static if (isIntegral!T || isSomeChar!T)
2335 foreach (immutable _; 0 .. 100)
2337 auto u = uniform!"[]"(T.min, T.max);
2338 static assert(is(typeof(u) == T));
2339 assert(T.min <= u, "Lower bound violation for uniform!\"[]\" with " ~ T.stringof);
2340 assert(u <= T.max, "Upper bound violation for uniform!\"[]\" with " ~ T.stringof);
2345 auto reproRng = Xorshift(239842);
2347 static foreach (T; std.meta.AliasSeq!(char, wchar, dchar, byte, ubyte, short,
2348 ushort, int, uint, long, ulong))
2350 T lo = T.min + 10, hi = T.max - 10;
2351 T init = uniform(lo, hi, reproRng);
2352 size_t i = 50;
2353 while (--i && uniform(lo, hi, reproRng) == init) {}
2354 assert(i > 0);
2358 bool sawLB = false, sawUB = false;
2359 foreach (i; 0 .. 50)
2361 auto x = uniform!"[]"('a', 'd', reproRng);
2362 if (x == 'a') sawLB = true;
2363 if (x == 'd') sawUB = true;
2364 assert('a' <= x && x <= 'd');
2366 assert(sawLB && sawUB);
2370 bool sawLB = false, sawUB = false;
2371 foreach (i; 0 .. 50)
2373 auto x = uniform('a', 'd', reproRng);
2374 if (x == 'a') sawLB = true;
2375 if (x == 'c') sawUB = true;
2376 assert('a' <= x && x < 'd');
2378 assert(sawLB && sawUB);
2382 bool sawLB = false, sawUB = false;
2383 foreach (i; 0 .. 50)
2385 immutable int lo = -2, hi = 2;
2386 auto x = uniform!"()"(lo, hi, reproRng);
2387 if (x == (lo+1)) sawLB = true;
2388 if (x == (hi-1)) sawUB = true;
2389 assert(lo < x && x < hi);
2391 assert(sawLB && sawUB);
2395 bool sawLB = false, sawUB = false;
2396 foreach (i; 0 .. 50)
2398 immutable ubyte lo = 0, hi = 5;
2399 auto x = uniform(lo, hi, reproRng);
2400 if (x == lo) sawLB = true;
2401 if (x == (hi-1)) sawUB = true;
2402 assert(lo <= x && x < hi);
2404 assert(sawLB && sawUB);
2408 foreach (i; 0 .. 30)
2410 assert(i == uniform(i, i+1, reproRng));
2416 Generates an unsigned integer in the half-open range `[0, k)`.
2417 Non-public because we locally guarantee `k > 0`.
2419 Params:
2420 k = unsigned exclusive upper bound; caller guarantees this is non-zero
2421 rng = random number generator to use
2423 Returns:
2424 Pseudo-random unsigned integer strictly less than `k`.
2426 private UInt _uniformIndex(UniformRNG, UInt = size_t)(const UInt k, ref UniformRNG rng)
2427 if (isUnsigned!UInt && isUniformRNG!UniformRNG)
2429 alias ResultType = UInt;
2430 alias UpperType = Unsigned!(typeof(k - 0));
2431 alias upperDist = k;
2433 assert(upperDist != 0);
2435 // For backwards compatibility use same algorithm as uniform(0, k, rng).
2436 UpperType offset, rnum, bucketFront;
2439 rnum = uniform!UpperType(rng);
2440 offset = rnum % upperDist;
2441 bucketFront = rnum - offset;
2442 } // while we're in an unfair bucket...
2443 while (bucketFront > (UpperType.max - (upperDist - 1)));
2445 return cast(ResultType) offset;
2448 pure @safe unittest
2450 // For backwards compatibility check that _uniformIndex(k, rng)
2451 // has the same result as uniform(0, k, rng).
2452 auto rng1 = Xorshift(123_456_789);
2453 auto rng2 = rng1.save();
2454 const size_t k = (1U << 31) - 1;
2455 assert(_uniformIndex(k, rng1) == uniform(0, k, rng2));
2459 Generates a uniformly-distributed number in the range $(D [T.min,
2460 T.max]) for any integral or character type `T`. If no random
2461 number generator is passed, uses the default `rndGen`.
2463 If an `enum` is used as type, the random variate is drawn with
2464 equal probability from any of the possible values of the enum `E`.
2466 Params:
2467 urng = (optional) random number generator to use;
2468 if not specified, defaults to `rndGen`
2470 Returns:
2471 Random variate drawn from the _uniform distribution across all
2472 possible values of the integral, character or enum type `T`.
2474 auto uniform(T, UniformRandomNumberGenerator)
2475 (ref UniformRandomNumberGenerator urng)
2476 if (!is(T == enum) && (isIntegral!T || isSomeChar!T) && isUniformRNG!UniformRandomNumberGenerator)
2478 /* dchar does not use its full bit range, so we must
2479 * revert to the uniform with specified bounds
2481 static if (is(immutable T == immutable dchar))
2483 return uniform!"[]"(T.min, T.max, urng);
2485 else
2487 auto r = urng.front;
2488 urng.popFront();
2489 static if (T.sizeof <= r.sizeof)
2491 return cast(T) r;
2493 else
2495 static assert(T.sizeof == 8 && r.sizeof == 4);
2496 T r1 = urng.front | (cast(T) r << 32);
2497 urng.popFront();
2498 return r1;
2503 /// Ditto
2504 auto uniform(T)()
2505 if (!is(T == enum) && (isIntegral!T || isSomeChar!T))
2507 return uniform!T(rndGen);
2511 @safe unittest
2513 auto rnd = MinstdRand0(42);
2515 assert(rnd.uniform!ubyte == 102);
2516 assert(rnd.uniform!ulong == 4838462006927449017);
2518 enum Fruit { apple, mango, pear }
2519 version (X86_64) // https://issues.dlang.org/show_bug.cgi?id=15147
2520 assert(rnd.uniform!Fruit == Fruit.mango);
2523 @safe unittest
2525 // https://issues.dlang.org/show_bug.cgi?id=21383
2526 auto rng1 = Xorshift32(123456789);
2527 auto rng2 = rng1.save;
2528 assert(rng1.uniform!dchar == rng2.uniform!dchar);
2529 // https://issues.dlang.org/show_bug.cgi?id=21384
2530 assert(rng1.uniform!(const shared dchar) <= dchar.max);
2531 // https://issues.dlang.org/show_bug.cgi?id=8671
2532 double t8671 = 1.0 - uniform(0.0, 1.0);
2535 @safe unittest
2537 static foreach (T; std.meta.AliasSeq!(char, wchar, dchar, byte, ubyte, short, ushort,
2538 int, uint, long, ulong))
2540 T init = uniform!T();
2541 size_t i = 50;
2542 while (--i && uniform!T() == init) {}
2543 assert(i > 0);
2545 foreach (immutable _; 0 .. 100)
2547 auto u = uniform!T();
2548 static assert(is(typeof(u) == T));
2549 assert(T.min <= u, "Lower bound violation for uniform!" ~ T.stringof);
2550 assert(u <= T.max, "Upper bound violation for uniform!" ~ T.stringof);
2555 /// ditto
2556 auto uniform(E, UniformRandomNumberGenerator)
2557 (ref UniformRandomNumberGenerator urng)
2558 if (is(E == enum) && isUniformRNG!UniformRandomNumberGenerator)
2560 static immutable E[EnumMembers!E.length] members = [EnumMembers!E];
2561 return members[std.random.uniform(0, members.length, urng)];
2564 /// Ditto
2565 auto uniform(E)()
2566 if (is(E == enum))
2568 return uniform!E(rndGen);
2571 @safe unittest
2573 enum Fruit { Apple = 12, Mango = 29, Pear = 72 }
2574 foreach (_; 0 .. 100)
2576 foreach (f; [uniform!Fruit(), rndGen.uniform!Fruit()])
2578 assert(f == Fruit.Apple || f == Fruit.Mango || f == Fruit.Pear);
2584 * Generates a uniformly-distributed floating point number of type
2585 * `T` in the range [0, 1$(RPAREN). If no random number generator is
2586 * specified, the default RNG `rndGen` will be used as the source
2587 * of randomness.
2589 * `uniform01` offers a faster generation of random variates than
2590 * the equivalent $(D uniform!"[$(RPAREN)"(0.0, 1.0)) and so may be preferred
2591 * for some applications.
2593 * Params:
2594 * rng = (optional) random number generator to use;
2595 * if not specified, defaults to `rndGen`
2597 * Returns:
2598 * Floating-point random variate of type `T` drawn from the _uniform
2599 * distribution across the half-open interval [0, 1$(RPAREN).
2602 T uniform01(T = double)()
2603 if (isFloatingPoint!T)
2605 return uniform01!T(rndGen);
2608 /// ditto
2609 T uniform01(T = double, UniformRNG)(ref UniformRNG rng)
2610 if (isFloatingPoint!T && isUniformRNG!UniformRNG)
2611 out (result)
2613 assert(0 <= result);
2614 assert(result < 1);
2618 alias R = typeof(rng.front);
2619 static if (isIntegral!R)
2621 enum T factor = 1 / (T(1) + rng.max - rng.min);
2623 else static if (isFloatingPoint!R)
2625 enum T factor = 1 / (rng.max - rng.min);
2627 else
2629 static assert(false);
2632 while (true)
2634 immutable T u = (rng.front - rng.min) * factor;
2635 rng.popFront();
2637 static if (isIntegral!R && T.mant_dig >= (8 * R.sizeof))
2639 /* If RNG variates are integral and T has enough precision to hold
2640 * R without loss, we're guaranteed by the definition of factor
2641 * that precisely u < 1.
2643 return u;
2645 else
2647 /* Otherwise we have to check whether u is beyond the assumed range
2648 * because of the loss of precision, or for another reason, a
2649 * floating-point RNG can return a variate that is exactly equal to
2650 * its maximum.
2652 if (u < 1)
2654 return u;
2659 // Shouldn't ever get here.
2660 assert(false);
2664 @safe @nogc unittest
2666 import std.math.operations : feqrel;
2668 auto rnd = MinstdRand0(42);
2670 // Generate random numbers in the range in the range [0, 1)
2671 auto u1 = uniform01(rnd);
2672 assert(u1 >= 0 && u1 < 1);
2674 auto u2 = rnd.uniform01!float;
2675 assert(u2 >= 0 && u2 < 1);
2677 // Confirm that the random values with the initial seed 42 are 0.000328707 and 0.524587
2678 assert(u1.feqrel(0.000328707) > 20);
2679 assert(u2.feqrel(0.524587) > 20);
2682 @safe @nogc unittest
2684 import std.meta;
2685 static foreach (UniformRNG; PseudoRngTypes)
2688 static foreach (T; std.meta.AliasSeq!(float, double, real))
2690 UniformRNG rng = UniformRNG(123_456_789);
2692 auto a = uniform01();
2693 assert(is(typeof(a) == double));
2694 assert(0 <= a && a < 1);
2696 auto b = uniform01(rng);
2697 assert(is(typeof(a) == double));
2698 assert(0 <= b && b < 1);
2700 auto c = uniform01!T();
2701 assert(is(typeof(c) == T));
2702 assert(0 <= c && c < 1);
2704 auto d = uniform01!T(rng);
2705 assert(is(typeof(d) == T));
2706 assert(0 <= d && d < 1);
2708 T init = uniform01!T(rng);
2709 size_t i = 50;
2710 while (--i && uniform01!T(rng) == init) {}
2711 assert(i > 0);
2712 assert(i < 50);
2718 Generates a uniform probability distribution of size `n`, i.e., an
2719 array of size `n` of positive numbers of type `F` that sum to
2720 `1`. If `useThis` is provided, it is used as storage.
2722 F[] uniformDistribution(F = double)(size_t n, F[] useThis = null)
2723 if (isFloatingPoint!F)
2725 import std.numeric : normalize;
2726 useThis.length = n;
2727 foreach (ref e; useThis)
2729 e = uniform(0.0, 1);
2731 normalize(useThis);
2732 return useThis;
2736 @safe unittest
2738 import std.algorithm.iteration : reduce;
2739 import std.math.operations : isClose;
2741 auto a = uniformDistribution(5);
2742 assert(a.length == 5);
2743 assert(isClose(reduce!"a + b"(a), 1));
2745 a = uniformDistribution(10, a);
2746 assert(a.length == 10);
2747 assert(isClose(reduce!"a + b"(a), 1));
2751 Returns a random, uniformly chosen, element `e` from the supplied
2752 $(D Range range). If no random number generator is passed, the default
2753 `rndGen` is used.
2755 Params:
2756 range = a random access range that has the `length` property defined
2757 urng = (optional) random number generator to use;
2758 if not specified, defaults to `rndGen`
2760 Returns:
2761 A single random element drawn from the `range`. If it can, it will
2762 return a `ref` to the $(D range element), otherwise it will return
2763 a copy.
2765 auto ref choice(Range, RandomGen = Random)(auto ref Range range, ref RandomGen urng)
2766 if (isRandomAccessRange!Range && hasLength!Range && isUniformRNG!RandomGen)
2768 assert(range.length > 0,
2769 __PRETTY_FUNCTION__ ~ ": invalid Range supplied. Range cannot be empty");
2771 return range[uniform(size_t(0), $, urng)];
2774 /// ditto
2775 auto ref choice(Range)(auto ref Range range)
2777 return choice(range, rndGen);
2781 @safe unittest
2783 auto rnd = MinstdRand0(42);
2785 auto elem = [1, 2, 3, 4, 5].choice(rnd);
2786 version (X86_64) // https://issues.dlang.org/show_bug.cgi?id=15147
2787 assert(elem == 3);
2790 @safe unittest
2792 import std.algorithm.searching : canFind;
2794 class MyTestClass
2796 int x;
2798 this(int x)
2800 this.x = x;
2804 MyTestClass[] testClass;
2805 foreach (i; 0 .. 5)
2807 testClass ~= new MyTestClass(i);
2810 auto elem = choice(testClass);
2812 assert(canFind!((ref MyTestClass a, ref MyTestClass b) => a.x == b.x)(testClass, elem),
2813 "Choice did not return a valid element from the given Range");
2816 @system unittest
2818 import std.algorithm.iteration : map;
2819 import std.algorithm.searching : canFind;
2821 auto array = [1, 2, 3, 4, 5];
2822 auto elemAddr = &choice(array);
2824 assert(array.map!((ref e) => &e).canFind(elemAddr),
2825 "Choice did not return a ref to an element from the given Range");
2826 assert(array.canFind(*(cast(int *)(elemAddr))),
2827 "Choice did not return a valid element from the given Range");
2831 Shuffles elements of `r` using `gen` as a shuffler. `r` must be
2832 a random-access range with length. If no RNG is specified, `rndGen`
2833 will be used.
2835 Params:
2836 r = random-access range whose elements are to be shuffled
2837 gen = (optional) random number generator to use; if not
2838 specified, defaults to `rndGen`
2839 Returns:
2840 The shuffled random-access range.
2843 Range randomShuffle(Range, RandomGen)(Range r, ref RandomGen gen)
2844 if (isRandomAccessRange!Range && isUniformRNG!RandomGen)
2846 import std.algorithm.mutation : swapAt;
2847 const n = r.length;
2848 foreach (i; 0 .. n)
2850 r.swapAt(i, i + _uniformIndex(n - i, gen));
2852 return r;
2855 /// ditto
2856 Range randomShuffle(Range)(Range r)
2857 if (isRandomAccessRange!Range)
2859 return randomShuffle(r, rndGen);
2863 @safe unittest
2865 auto rnd = MinstdRand0(42);
2867 auto arr = [1, 2, 3, 4, 5].randomShuffle(rnd);
2868 version (X86_64) // https://issues.dlang.org/show_bug.cgi?id=15147
2869 assert(arr == [3, 5, 2, 4, 1]);
2872 @safe unittest
2874 int[10] sa = void;
2875 int[10] sb = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9];
2876 import std.algorithm.sorting : sort;
2877 foreach (RandomGen; PseudoRngTypes)
2879 sa[] = sb[];
2880 auto a = sa[];
2881 auto b = sb[];
2882 auto gen = RandomGen(123_456_789);
2883 randomShuffle(a, gen);
2884 sort(a);
2885 assert(a == b);
2886 randomShuffle(a);
2887 sort(a);
2888 assert(a == b);
2890 // For backwards compatibility verify randomShuffle(r, gen)
2891 // is equivalent to partialShuffle(r, 0, r.length, gen).
2892 auto gen1 = Xorshift(123_456_789);
2893 auto gen2 = gen1.save();
2894 sa[] = sb[];
2895 // @nogc std.random.randomShuffle.
2896 // https://issues.dlang.org/show_bug.cgi?id=19156
2897 () @nogc nothrow pure { randomShuffle(sa[], gen1); }();
2898 partialShuffle(sb[], sb.length, gen2);
2899 assert(sa[] == sb[]);
2902 // https://issues.dlang.org/show_bug.cgi?id=18501
2903 @safe unittest
2905 import std.algorithm.comparison : among;
2906 auto r = randomShuffle([0,1]);
2907 assert(r.among([0,1],[1,0]));
2911 Partially shuffles the elements of `r` such that upon returning $(D r[0 .. n])
2912 is a random subset of `r` and is randomly ordered. $(D r[n .. r.length])
2913 will contain the elements not in $(D r[0 .. n]). These will be in an undefined
2914 order, but will not be random in the sense that their order after
2915 `partialShuffle` returns will not be independent of their order before
2916 `partialShuffle` was called.
2918 `r` must be a random-access range with length. `n` must be less than
2919 or equal to `r.length`. If no RNG is specified, `rndGen` will be used.
2921 Params:
2922 r = random-access range whose elements are to be shuffled
2923 n = number of elements of `r` to shuffle (counting from the beginning);
2924 must be less than `r.length`
2925 gen = (optional) random number generator to use; if not
2926 specified, defaults to `rndGen`
2927 Returns:
2928 The shuffled random-access range.
2930 Range partialShuffle(Range, RandomGen)(Range r, in size_t n, ref RandomGen gen)
2931 if (isRandomAccessRange!Range && isUniformRNG!RandomGen)
2933 import std.algorithm.mutation : swapAt;
2934 import std.exception : enforce;
2935 enforce(n <= r.length, "n must be <= r.length for partialShuffle.");
2936 foreach (i; 0 .. n)
2938 r.swapAt(i, uniform(i, r.length, gen));
2940 return r;
2943 /// ditto
2944 Range partialShuffle(Range)(Range r, in size_t n)
2945 if (isRandomAccessRange!Range)
2947 return partialShuffle(r, n, rndGen);
2951 @safe unittest
2953 auto rnd = MinstdRand0(42);
2955 auto arr = [1, 2, 3, 4, 5, 6];
2956 arr = arr.dup.partialShuffle(1, rnd);
2958 version (X86_64) // https://issues.dlang.org/show_bug.cgi?id=15147
2959 assert(arr == [2, 1, 3, 4, 5, 6]); // 1<->2
2961 arr = arr.dup.partialShuffle(2, rnd);
2962 version (X86_64) // https://issues.dlang.org/show_bug.cgi?id=15147
2963 assert(arr == [1, 4, 3, 2, 5, 6]); // 1<->2, 2<->4
2965 arr = arr.dup.partialShuffle(3, rnd);
2966 version (X86_64) // https://issues.dlang.org/show_bug.cgi?id=15147
2967 assert(arr == [5, 4, 6, 2, 1, 3]); // 1<->5, 2<->4, 3<->6
2970 @safe unittest
2972 import std.algorithm;
2973 foreach (RandomGen; PseudoRngTypes)
2975 auto a = [0, 1, 1, 2, 3];
2976 auto b = a.dup;
2978 // Pick a fixed seed so that the outcome of the statistical
2979 // test below is deterministic.
2980 auto gen = RandomGen(12345);
2982 // NUM times, pick LEN elements from the array at random.
2983 immutable int LEN = 2;
2984 immutable int NUM = 750;
2985 int[][] chk;
2986 foreach (step; 0 .. NUM)
2988 partialShuffle(a, LEN, gen);
2989 chk ~= a[0 .. LEN].dup;
2992 // Check that each possible a[0 .. LEN] was produced at least once.
2993 // For a perfectly random RandomGen, the probability that each
2994 // particular combination failed to appear would be at most
2995 // 0.95 ^^ NUM which is approximately 1,962e-17.
2996 // As long as hardware failure (e.g. bit flip) probability
2997 // is higher, we are fine with this unittest.
2998 sort(chk);
2999 assert(equal(uniq(chk), [ [0,1], [0,2], [0,3],
3000 [1,0], [1,1], [1,2], [1,3],
3001 [2,0], [2,1], [2,3],
3002 [3,0], [3,1], [3,2], ]));
3004 // Check that all the elements are still there.
3005 sort(a);
3006 assert(equal(a, b));
3011 Rolls a dice with relative probabilities stored in $(D
3012 proportions). Returns the index in `proportions` that was chosen.
3014 Params:
3015 rnd = (optional) random number generator to use; if not
3016 specified, defaults to `rndGen`
3017 proportions = forward range or list of individual values
3018 whose elements correspond to the probabilities
3019 with which to choose the corresponding index
3020 value
3022 Returns:
3023 Random variate drawn from the index values
3024 [0, ... `proportions.length` - 1], with the probability
3025 of getting an individual index value `i` being proportional to
3026 `proportions[i]`.
3028 size_t dice(Rng, Num)(ref Rng rnd, Num[] proportions...)
3029 if (isNumeric!Num && isForwardRange!Rng)
3031 return diceImpl(rnd, proportions);
3034 /// Ditto
3035 size_t dice(R, Range)(ref R rnd, Range proportions)
3036 if (isForwardRange!Range && isNumeric!(ElementType!Range) && !isArray!Range)
3038 return diceImpl(rnd, proportions);
3041 /// Ditto
3042 size_t dice(Range)(Range proportions)
3043 if (isForwardRange!Range && isNumeric!(ElementType!Range) && !isArray!Range)
3045 return diceImpl(rndGen, proportions);
3048 /// Ditto
3049 size_t dice(Num)(Num[] proportions...)
3050 if (isNumeric!Num)
3052 return diceImpl(rndGen, proportions);
3056 @safe unittest
3058 auto x = dice(0.5, 0.5); // x is 0 or 1 in equal proportions
3059 auto y = dice(50, 50); // y is 0 or 1 in equal proportions
3060 auto z = dice(70, 20, 10); // z is 0 70% of the time, 1 20% of the time,
3061 // and 2 10% of the time
3065 @safe unittest
3067 auto rnd = MinstdRand0(42);
3068 auto z = rnd.dice(70, 20, 10);
3069 assert(z == 0);
3070 z = rnd.dice(30, 20, 40, 10);
3071 assert(z == 2);
3074 private size_t diceImpl(Rng, Range)(ref Rng rng, scope Range proportions)
3075 if (isForwardRange!Range && isNumeric!(ElementType!Range) && isForwardRange!Rng)
3078 import std.algorithm.searching : all;
3079 assert(proportions.save.all!"a >= 0");
3083 import std.algorithm.iteration : reduce;
3084 import std.exception : enforce;
3085 double sum = reduce!"a + b"(0.0, proportions.save);
3086 enforce(sum > 0, "Proportions in a dice cannot sum to zero");
3087 immutable point = uniform(0.0, sum, rng);
3088 assert(point < sum);
3089 auto mass = 0.0;
3091 size_t i = 0;
3092 foreach (e; proportions)
3094 mass += e;
3095 if (point < mass) return i;
3096 i++;
3098 // this point should not be reached
3099 assert(false);
3103 @safe unittest
3105 auto rnd = Xorshift(123_456_789);
3106 auto i = dice(rnd, 0.0, 100.0);
3107 assert(i == 1);
3108 i = dice(rnd, 100.0, 0.0);
3109 assert(i == 0);
3111 i = dice(100U, 0U);
3112 assert(i == 0);
3115 /+ @nogc bool array designed for RandomCover.
3116 - constructed with an invariable length
3117 - small length means 0 alloc and bit field (if up to 32(x86) or 64(x64) choices to cover)
3118 - bigger length means non-GC heap allocation(s) and dealloc. +/
3119 private struct RandomCoverChoices
3121 private size_t* buffer;
3122 private immutable size_t _length;
3123 private immutable bool hasPackedBits;
3124 private enum BITS_PER_WORD = typeof(buffer[0]).sizeof * 8;
3126 void opAssign(T)(T) @disable;
3128 this(this) pure nothrow @nogc @trusted
3130 import core.stdc.string : memcpy;
3131 import std.internal.memory : enforceMalloc;
3133 if (!hasPackedBits && buffer !is null)
3135 const nBytesToAlloc = size_t.sizeof * (_length / BITS_PER_WORD + int(_length % BITS_PER_WORD != 0));
3136 void* nbuffer = enforceMalloc(nBytesToAlloc);
3137 buffer = cast(size_t*) memcpy(nbuffer, buffer, nBytesToAlloc);
3141 this(size_t numChoices) pure nothrow @nogc @trusted
3143 import std.internal.memory : enforceCalloc;
3145 _length = numChoices;
3146 hasPackedBits = _length <= size_t.sizeof * 8;
3147 if (!hasPackedBits)
3149 const nWordsToAlloc = _length / BITS_PER_WORD + int(_length % BITS_PER_WORD != 0);
3150 buffer = cast(size_t*) enforceCalloc(nWordsToAlloc, BITS_PER_WORD / 8);
3154 size_t length() const pure nothrow @nogc @safe @property {return _length;}
3156 ~this() pure nothrow @nogc @trusted
3158 import core.memory : pureFree;
3160 if (!hasPackedBits && buffer !is null)
3161 pureFree(buffer);
3164 bool opIndex(size_t index) const pure nothrow @nogc @trusted
3166 assert(index < _length);
3167 import core.bitop : bt;
3168 if (!hasPackedBits)
3169 return cast(bool) bt(buffer, index);
3170 else
3171 return ((cast(size_t) buffer) >> index) & size_t(1);
3174 void opIndexAssign(bool value, size_t index) pure nothrow @nogc @trusted
3176 assert(index < _length);
3177 if (!hasPackedBits)
3179 import core.bitop : btr, bts;
3180 if (value)
3181 bts(buffer, index);
3182 else
3183 btr(buffer, index);
3185 else
3187 if (value)
3188 (*cast(size_t*) &buffer) |= size_t(1) << index;
3189 else
3190 (*cast(size_t*) &buffer) &= ~(size_t(1) << index);
3195 @safe @nogc nothrow unittest
3197 static immutable lengths = [3, 32, 65, 256];
3198 foreach (length; lengths)
3200 RandomCoverChoices c = RandomCoverChoices(length);
3201 assert(c.hasPackedBits == (length <= size_t.sizeof * 8));
3202 c[0] = true;
3203 c[2] = true;
3204 assert(c[0]);
3205 assert(!c[1]);
3206 assert(c[2]);
3207 c[0] = false;
3208 c[1] = true;
3209 c[2] = false;
3210 assert(!c[0]);
3211 assert(c[1]);
3212 assert(!c[2]);
3217 Covers a given range `r` in a random manner, i.e. goes through each
3218 element of `r` once and only once, just in a random order. `r`
3219 must be a random-access range with length.
3221 If no random number generator is passed to `randomCover`, the
3222 thread-global RNG rndGen will be used internally.
3224 Params:
3225 r = random-access range to cover
3226 rng = (optional) random number generator to use;
3227 if not specified, defaults to `rndGen`
3229 Returns:
3230 Range whose elements consist of the elements of `r`,
3231 in random order. Will be a forward range if both `r` and
3232 `rng` are forward ranges, an
3233 $(REF_ALTTEXT input range, isInputRange, std,range,primitives) otherwise.
3235 struct RandomCover(Range, UniformRNG = void)
3236 if (isRandomAccessRange!Range && (isUniformRNG!UniformRNG || is(UniformRNG == void)))
3238 private Range _input;
3239 private RandomCoverChoices _chosen;
3240 private size_t _current;
3241 private size_t _alreadyChosen = 0;
3242 private bool _isEmpty = false;
3244 static if (is(UniformRNG == void))
3246 this(Range input)
3248 _input = input;
3249 _chosen = RandomCoverChoices(_input.length);
3250 if (_input.empty)
3252 _isEmpty = true;
3254 else
3256 _current = _uniformIndex(_chosen.length, rndGen);
3260 else
3262 private UniformRNG _rng;
3264 this(Range input, ref UniformRNG rng)
3266 _input = input;
3267 _rng = rng;
3268 _chosen = RandomCoverChoices(_input.length);
3269 if (_input.empty)
3271 _isEmpty = true;
3273 else
3275 _current = _uniformIndex(_chosen.length, rng);
3279 this(Range input, UniformRNG rng)
3281 this(input, rng);
3285 static if (hasLength!Range)
3287 @property size_t length()
3289 return _input.length - _alreadyChosen;
3293 @property auto ref front()
3295 assert(!_isEmpty);
3296 return _input[_current];
3299 void popFront()
3301 assert(!_isEmpty);
3303 size_t k = _input.length - _alreadyChosen - 1;
3304 if (k == 0)
3306 _isEmpty = true;
3307 ++_alreadyChosen;
3308 return;
3311 size_t i;
3312 foreach (e; _input)
3314 if (_chosen[i] || i == _current) { ++i; continue; }
3315 // Roll a dice with k faces
3316 static if (is(UniformRNG == void))
3318 auto chooseMe = _uniformIndex(k, rndGen) == 0;
3320 else
3322 auto chooseMe = _uniformIndex(k, _rng) == 0;
3324 assert(k > 1 || chooseMe);
3325 if (chooseMe)
3327 _chosen[_current] = true;
3328 _current = i;
3329 ++_alreadyChosen;
3330 return;
3332 --k;
3333 ++i;
3337 static if (isForwardRange!UniformRNG)
3339 @property typeof(this) save()
3341 auto ret = this;
3342 ret._input = _input.save;
3343 ret._rng = _rng.save;
3344 return ret;
3348 @property bool empty() const { return _isEmpty; }
3351 /// Ditto
3352 auto randomCover(Range, UniformRNG)(Range r, auto ref UniformRNG rng)
3353 if (isRandomAccessRange!Range && isUniformRNG!UniformRNG)
3355 return RandomCover!(Range, UniformRNG)(r, rng);
3358 /// Ditto
3359 auto randomCover(Range)(Range r)
3360 if (isRandomAccessRange!Range)
3362 return RandomCover!(Range, void)(r);
3366 @safe unittest
3368 import std.algorithm.comparison : equal;
3369 import std.range : iota;
3370 auto rnd = MinstdRand0(42);
3372 version (X86_64) // https://issues.dlang.org/show_bug.cgi?id=15147
3373 assert(10.iota.randomCover(rnd).equal([7, 4, 2, 0, 1, 6, 8, 3, 9, 5]));
3376 @safe unittest // cover RandomCoverChoices postblit for heap storage
3378 import std.array : array;
3379 import std.range : iota;
3380 auto a = 1337.iota.randomCover().array;
3381 assert(a.length == 1337);
3384 @nogc nothrow pure @safe unittest
3386 // Optionally @nogc std.random.randomCover
3387 // https://issues.dlang.org/show_bug.cgi?id=14001
3388 auto rng = Xorshift(123_456_789);
3389 static immutable int[] sa = [1, 2, 3, 4, 5];
3390 auto r = randomCover(sa, rng);
3391 assert(!r.empty);
3392 const x = r.front;
3393 r.popFront();
3394 assert(!r.empty);
3395 const y = r.front;
3396 assert(x != y);
3399 @safe unittest
3401 import std.algorithm;
3402 import std.conv;
3403 int[] a = [ 0, 1, 2, 3, 4, 5, 6, 7, 8 ];
3404 int[] c;
3405 static foreach (UniformRNG; std.meta.AliasSeq!(void, PseudoRngTypes))
3407 static if (is(UniformRNG == void))
3409 auto rc = randomCover(a);
3410 static assert(isInputRange!(typeof(rc)));
3411 static assert(!isForwardRange!(typeof(rc)));
3413 else
3415 auto rng = UniformRNG(123_456_789);
3416 auto rc = randomCover(a, rng);
3417 static assert(isForwardRange!(typeof(rc)));
3418 // check for constructor passed a value-type RNG
3419 auto rc2 = RandomCover!(int[], UniformRNG)(a, UniformRNG(987_654_321));
3420 static assert(isForwardRange!(typeof(rc2)));
3421 auto rcEmpty = randomCover(c, rng);
3422 assert(rcEmpty.length == 0);
3425 int[] b = new int[9];
3426 uint i;
3427 foreach (e; rc)
3429 //writeln(e);
3430 b[i++] = e;
3432 sort(b);
3433 assert(a == b, text(b));
3437 @safe unittest
3439 // https://issues.dlang.org/show_bug.cgi?id=12589
3440 int[] r = [];
3441 auto rc = randomCover(r);
3442 assert(rc.length == 0);
3443 assert(rc.empty);
3445 // https://issues.dlang.org/show_bug.cgi?id=16724
3446 import std.range : iota;
3447 auto range = iota(10);
3448 auto randy = range.randomCover;
3450 for (int i=1; i <= range.length; i++)
3452 randy.popFront;
3453 assert(randy.length == range.length - i);
3457 // RandomSample
3459 Selects a random subsample out of `r`, containing exactly `n`
3460 elements. The order of elements is the same as in the original
3461 range. The total length of `r` must be known. If `total` is
3462 passed in, the total number of sample is considered to be $(D
3463 total). Otherwise, `RandomSample` uses `r.length`.
3465 Params:
3466 r = range to sample from
3467 n = number of elements to include in the sample;
3468 must be less than or equal to the total number
3469 of elements in `r` and/or the parameter
3470 `total` (if provided)
3471 total = (semi-optional) number of elements of `r`
3472 from which to select the sample (counting from
3473 the beginning); must be less than or equal to
3474 the total number of elements in `r` itself.
3475 May be omitted if `r` has the `.length`
3476 property and the sample is to be drawn from
3477 all elements of `r`.
3478 rng = (optional) random number generator to use;
3479 if not specified, defaults to `rndGen`
3481 Returns:
3482 Range whose elements consist of a randomly selected subset of
3483 the elements of `r`, in the same order as these elements
3484 appear in `r` itself. Will be a forward range if both `r`
3485 and `rng` are forward ranges, an input range otherwise.
3487 `RandomSample` implements Jeffrey Scott Vitter's Algorithm D
3488 (see Vitter $(HTTP dx.doi.org/10.1145/358105.893, 1984), $(HTTP
3489 dx.doi.org/10.1145/23002.23003, 1987)), which selects a sample
3490 of size `n` in O(n) steps and requiring O(n) random variates,
3491 regardless of the size of the data being sampled. The exception
3492 to this is if traversing k elements on the input range is itself
3493 an O(k) operation (e.g. when sampling lines from an input file),
3494 in which case the sampling calculation will inevitably be of
3495 O(total).
3497 RandomSample will throw an exception if `total` is verifiably
3498 less than the total number of elements available in the input,
3499 or if $(D n > total).
3501 If no random number generator is passed to `randomSample`, the
3502 thread-global RNG rndGen will be used internally.
3504 struct RandomSample(Range, UniformRNG = void)
3505 if (isInputRange!Range && (isUniformRNG!UniformRNG || is(UniformRNG == void)))
3507 private size_t _available, _toSelect;
3508 private enum ushort _alphaInverse = 13; // Vitter's recommended value.
3509 private double _Vprime;
3510 private Range _input;
3511 private size_t _index;
3512 private enum Skip { None, A, D }
3513 private Skip _skip = Skip.None;
3515 // If we're using the default thread-local random number generator then
3516 // we shouldn't store a copy of it here. UniformRNG == void is a sentinel
3517 // for this. If we're using a user-specified generator then we have no
3518 // choice but to store a copy.
3519 static if (is(UniformRNG == void))
3521 static if (hasLength!Range)
3523 this(Range input, size_t howMany)
3525 _input = input;
3526 initialize(howMany, input.length);
3530 this(Range input, size_t howMany, size_t total)
3532 _input = input;
3533 initialize(howMany, total);
3536 else
3538 UniformRNG _rng;
3540 static if (hasLength!Range)
3542 this(Range input, size_t howMany, ref scope UniformRNG rng)
3544 _rng = rng;
3545 _input = input;
3546 initialize(howMany, input.length);
3549 this(Range input, size_t howMany, UniformRNG rng)
3551 this(input, howMany, rng);
3555 this(Range input, size_t howMany, size_t total, ref scope UniformRNG rng)
3557 _rng = rng;
3558 _input = input;
3559 initialize(howMany, total);
3562 this(Range input, size_t howMany, size_t total, UniformRNG rng)
3564 this(input, howMany, total, rng);
3568 private void initialize(size_t howMany, size_t total)
3570 import std.conv : text;
3571 import std.exception : enforce;
3572 _available = total;
3573 _toSelect = howMany;
3574 enforce(_toSelect <= _available,
3575 text("RandomSample: cannot sample ", _toSelect,
3576 " items when only ", _available, " are available"));
3577 static if (hasLength!Range)
3579 enforce(_available <= _input.length,
3580 text("RandomSample: specified ", _available,
3581 " items as available when input contains only ",
3582 _input.length));
3586 private void initializeFront()
3588 assert(_skip == Skip.None);
3589 // We can save ourselves a random variate by checking right
3590 // at the beginning if we should use Algorithm A.
3591 if ((_alphaInverse * _toSelect) > _available)
3593 _skip = Skip.A;
3595 else
3597 _skip = Skip.D;
3598 _Vprime = newVprime(_toSelect);
3600 prime();
3604 Range primitives.
3606 @property bool empty() const
3608 return _toSelect == 0;
3611 /// Ditto
3612 @property auto ref front()
3614 assert(!empty);
3615 // The first sample point must be determined here to avoid
3616 // having it always correspond to the first element of the
3617 // input. The rest of the sample points are determined each
3618 // time we call popFront().
3619 if (_skip == Skip.None)
3621 initializeFront();
3623 return _input.front;
3626 /// Ditto
3627 void popFront()
3629 // First we need to check if the sample has
3630 // been initialized in the first place.
3631 if (_skip == Skip.None)
3633 initializeFront();
3636 _input.popFront();
3637 --_available;
3638 --_toSelect;
3639 ++_index;
3640 prime();
3643 /// Ditto
3644 static if (isForwardRange!Range && isForwardRange!UniformRNG)
3646 static if (is(typeof(((const UniformRNG* p) => (*p).save)(null)) : UniformRNG)
3647 && is(typeof(((const Range* p) => (*p).save)(null)) : Range))
3649 @property typeof(this) save() const
3651 auto ret = RandomSample.init;
3652 foreach (fieldIndex, ref val; this.tupleof)
3654 static if (is(typeof(val) == const(Range)) || is(typeof(val) == const(UniformRNG)))
3655 ret.tupleof[fieldIndex] = val.save;
3656 else
3657 ret.tupleof[fieldIndex] = val;
3659 return ret;
3662 else
3664 @property typeof(this) save()
3666 auto ret = this;
3667 ret._input = _input.save;
3668 ret._rng = _rng.save;
3669 return ret;
3674 /// Ditto
3675 @property size_t length() const
3677 return _toSelect;
3681 Returns the index of the visited record.
3683 @property size_t index()
3685 if (_skip == Skip.None)
3687 initializeFront();
3689 return _index;
3692 private size_t skip()
3694 assert(_skip != Skip.None);
3696 // Step D1: if the number of points still to select is greater
3697 // than a certain proportion of the remaining data points, i.e.
3698 // if n >= alpha * N where alpha = 1/13, we carry out the
3699 // sampling with Algorithm A.
3700 if (_skip == Skip.A)
3702 return skipA();
3704 else if ((_alphaInverse * _toSelect) > _available)
3706 // We shouldn't get here unless the current selected
3707 // algorithm is D.
3708 assert(_skip == Skip.D);
3709 _skip = Skip.A;
3710 return skipA();
3712 else
3714 assert(_skip == Skip.D);
3715 return skipD();
3720 Vitter's Algorithm A, used when the ratio of needed sample values
3721 to remaining data values is sufficiently large.
3723 private size_t skipA()
3725 size_t s;
3726 double v, quot, top;
3728 if (_toSelect == 1)
3730 static if (is(UniformRNG == void))
3732 s = uniform(0, _available);
3734 else
3736 s = uniform(0, _available, _rng);
3739 else
3741 v = 0;
3742 top = _available - _toSelect;
3743 quot = top / _available;
3745 static if (is(UniformRNG == void))
3747 v = uniform!"()"(0.0, 1.0);
3749 else
3751 v = uniform!"()"(0.0, 1.0, _rng);
3754 while (quot > v)
3756 ++s;
3757 quot *= (top - s) / (_available - s);
3761 return s;
3765 Randomly reset the value of _Vprime.
3767 private double newVprime(size_t remaining)
3769 static if (is(UniformRNG == void))
3771 double r = uniform!"()"(0.0, 1.0);
3773 else
3775 double r = uniform!"()"(0.0, 1.0, _rng);
3778 return r ^^ (1.0 / remaining);
3782 Vitter's Algorithm D. For an extensive description of the algorithm
3783 and its rationale, see:
3785 * Vitter, J.S. (1984), "Faster methods for random sampling",
3786 Commun. ACM 27(7): 703--718
3788 * Vitter, J.S. (1987) "An efficient algorithm for sequential random
3789 sampling", ACM Trans. Math. Softw. 13(1): 58-67.
3791 Variable names are chosen to match those in Vitter's paper.
3793 private size_t skipD()
3795 import std.math.traits : isNaN;
3796 import std.math.rounding : trunc;
3797 // Confirm that the check in Step D1 is valid and we
3798 // haven't been sent here by mistake
3799 assert((_alphaInverse * _toSelect) <= _available);
3801 // Now it's safe to use the standard Algorithm D mechanism.
3802 if (_toSelect > 1)
3804 size_t s;
3805 size_t qu1 = 1 + _available - _toSelect;
3806 double x, y1;
3808 assert(!_Vprime.isNaN());
3810 while (true)
3812 // Step D2: set values of x and u.
3813 while (1)
3815 x = _available * (1-_Vprime);
3816 s = cast(size_t) trunc(x);
3817 if (s < qu1)
3818 break;
3819 _Vprime = newVprime(_toSelect);
3822 static if (is(UniformRNG == void))
3824 double u = uniform!"()"(0.0, 1.0);
3826 else
3828 double u = uniform!"()"(0.0, 1.0, _rng);
3831 y1 = (u * (cast(double) _available) / qu1) ^^ (1.0/(_toSelect - 1));
3833 _Vprime = y1 * ((-x/_available)+1.0) * ( qu1/( (cast(double) qu1) - s ) );
3835 // Step D3: if _Vprime <= 1.0 our work is done and we return S.
3836 // Otherwise ...
3837 if (_Vprime > 1.0)
3839 size_t top = _available - 1, limit;
3840 double y2 = 1.0, bottom;
3842 if (_toSelect > (s+1))
3844 bottom = _available - _toSelect;
3845 limit = _available - s;
3847 else
3849 bottom = _available - (s+1);
3850 limit = qu1;
3853 foreach (size_t t; limit .. _available)
3855 y2 *= top/bottom;
3856 top--;
3857 bottom--;
3860 // Step D4: decide whether or not to accept the current value of S.
3861 if (_available/(_available-x) < y1 * (y2 ^^ (1.0/(_toSelect-1))))
3863 // If it's not acceptable, we generate a new value of _Vprime
3864 // and go back to the start of the for (;;) loop.
3865 _Vprime = newVprime(_toSelect);
3867 else
3869 // If it's acceptable we generate a new value of _Vprime
3870 // based on the remaining number of sample points needed,
3871 // and return S.
3872 _Vprime = newVprime(_toSelect-1);
3873 return s;
3876 else
3878 // Return if condition D3 satisfied.
3879 return s;
3883 else
3885 // If only one sample point remains to be taken ...
3886 return cast(size_t) trunc(_available * _Vprime);
3890 private void prime()
3892 if (empty)
3894 return;
3896 assert(_available && _available >= _toSelect);
3897 immutable size_t s = skip();
3898 assert(s + _toSelect <= _available);
3899 static if (hasLength!Range)
3901 assert(s + _toSelect <= _input.length);
3903 assert(!_input.empty);
3904 _input.popFrontExactly(s);
3905 _index += s;
3906 _available -= s;
3907 assert(_available > 0);
3911 /// Ditto
3912 auto randomSample(Range)(Range r, size_t n, size_t total)
3913 if (isInputRange!Range)
3915 return RandomSample!(Range, void)(r, n, total);
3918 /// Ditto
3919 auto randomSample(Range)(Range r, size_t n)
3920 if (isInputRange!Range && hasLength!Range)
3922 return RandomSample!(Range, void)(r, n, r.length);
3925 /// Ditto
3926 auto randomSample(Range, UniformRNG)(Range r, size_t n, size_t total, auto ref UniformRNG rng)
3927 if (isInputRange!Range && isUniformRNG!UniformRNG)
3929 return RandomSample!(Range, UniformRNG)(r, n, total, rng);
3932 /// Ditto
3933 auto randomSample(Range, UniformRNG)(Range r, size_t n, auto ref UniformRNG rng)
3934 if (isInputRange!Range && hasLength!Range && isUniformRNG!UniformRNG)
3936 return RandomSample!(Range, UniformRNG)(r, n, r.length, rng);
3940 @safe unittest
3942 import std.algorithm.comparison : equal;
3943 import std.range : iota;
3944 auto rnd = MinstdRand0(42);
3945 assert(10.iota.randomSample(3, rnd).equal([7, 8, 9]));
3948 @system unittest
3950 // @system because it takes the address of a local
3951 import std.conv : text;
3952 import std.exception;
3953 import std.range;
3954 // For test purposes, an infinite input range
3955 struct TestInputRange
3957 private auto r = recurrence!"a[n-1] + 1"(0);
3958 bool empty() @property const pure nothrow { return r.empty; }
3959 auto front() @property pure nothrow { return r.front; }
3960 void popFront() pure nothrow { r.popFront(); }
3962 static assert(isInputRange!TestInputRange);
3963 static assert(!isForwardRange!TestInputRange);
3965 const(int)[] a = [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 ];
3967 foreach (UniformRNG; PseudoRngTypes)
3968 (){ // avoid workaround optimizations for large functions
3969 // https://issues.dlang.org/show_bug.cgi?id=2396
3970 auto rng = UniformRNG(1234);
3971 /* First test the most general case: randomSample of input range, with and
3972 * without a specified random number generator.
3974 static assert(isInputRange!(typeof(randomSample(TestInputRange(), 5, 10))));
3975 static assert(isInputRange!(typeof(randomSample(TestInputRange(), 5, 10, rng))));
3976 static assert(!isForwardRange!(typeof(randomSample(TestInputRange(), 5, 10))));
3977 static assert(!isForwardRange!(typeof(randomSample(TestInputRange(), 5, 10, rng))));
3978 // test case with range initialized by direct call to struct
3980 auto sample =
3981 RandomSample!(TestInputRange, UniformRNG)
3982 (TestInputRange(), 5, 10, UniformRNG(987_654_321));
3983 static assert(isInputRange!(typeof(sample)));
3984 static assert(!isForwardRange!(typeof(sample)));
3987 /* Now test the case of an input range with length. We ignore the cases
3988 * already covered by the previous tests.
3990 static assert(isInputRange!(typeof(randomSample(TestInputRange().takeExactly(10), 5))));
3991 static assert(isInputRange!(typeof(randomSample(TestInputRange().takeExactly(10), 5, rng))));
3992 static assert(!isForwardRange!(typeof(randomSample(TestInputRange().takeExactly(10), 5))));
3993 static assert(!isForwardRange!(typeof(randomSample(TestInputRange().takeExactly(10), 5, rng))));
3994 // test case with range initialized by direct call to struct
3996 auto sample =
3997 RandomSample!(typeof(TestInputRange().takeExactly(10)), UniformRNG)
3998 (TestInputRange().takeExactly(10), 5, 10, UniformRNG(654_321_987));
3999 static assert(isInputRange!(typeof(sample)));
4000 static assert(!isForwardRange!(typeof(sample)));
4003 // Now test the case of providing a forward range as input.
4004 static assert(!isForwardRange!(typeof(randomSample(a, 5))));
4005 static if (isForwardRange!UniformRNG)
4007 static assert(isForwardRange!(typeof(randomSample(a, 5, rng))));
4008 // ... and test with range initialized directly
4010 auto sample =
4011 RandomSample!(const(int)[], UniformRNG)
4012 (a, 5, UniformRNG(321_987_654));
4013 static assert(isForwardRange!(typeof(sample)));
4016 else
4018 static assert(isInputRange!(typeof(randomSample(a, 5, rng))));
4019 static assert(!isForwardRange!(typeof(randomSample(a, 5, rng))));
4020 // ... and test with range initialized directly
4022 auto sample =
4023 RandomSample!(const(int)[], UniformRNG)
4024 (a, 5, UniformRNG(789_123_456));
4025 static assert(isInputRange!(typeof(sample)));
4026 static assert(!isForwardRange!(typeof(sample)));
4030 /* Check that randomSample will throw an error if we claim more
4031 * items are available than there actually are, or if we try to
4032 * sample more items than are available. */
4033 assert(collectExceptionMsg(
4034 randomSample(a, 5, 15)
4035 ) == "RandomSample: specified 15 items as available when input contains only 10");
4036 assert(collectExceptionMsg(
4037 randomSample(a, 15)
4038 ) == "RandomSample: cannot sample 15 items when only 10 are available");
4039 assert(collectExceptionMsg(
4040 randomSample(a, 9, 8)
4041 ) == "RandomSample: cannot sample 9 items when only 8 are available");
4042 assert(collectExceptionMsg(
4043 randomSample(TestInputRange(), 12, 11)
4044 ) == "RandomSample: cannot sample 12 items when only 11 are available");
4046 /* Check that sampling algorithm never accidentally overruns the end of
4047 * the input range. If input is an InputRange without .length, this
4048 * relies on the user specifying the total number of available items
4049 * correctly.
4052 uint i = 0;
4053 foreach (e; randomSample(a, a.length))
4055 assert(e == i);
4056 ++i;
4058 assert(i == a.length);
4060 i = 0;
4061 foreach (e; randomSample(TestInputRange(), 17, 17))
4063 assert(e == i);
4064 ++i;
4066 assert(i == 17);
4070 // Check length properties of random samples.
4071 assert(randomSample(a, 5).length == 5);
4072 assert(randomSample(a, 5, 10).length == 5);
4073 assert(randomSample(a, 5, rng).length == 5);
4074 assert(randomSample(a, 5, 10, rng).length == 5);
4075 assert(randomSample(TestInputRange(), 5, 10).length == 5);
4076 assert(randomSample(TestInputRange(), 5, 10, rng).length == 5);
4078 // ... and emptiness!
4079 assert(randomSample(a, 0).empty);
4080 assert(randomSample(a, 0, 5).empty);
4081 assert(randomSample(a, 0, rng).empty);
4082 assert(randomSample(a, 0, 5, rng).empty);
4083 assert(randomSample(TestInputRange(), 0, 10).empty);
4084 assert(randomSample(TestInputRange(), 0, 10, rng).empty);
4086 /* Test that the (lazy) evaluation of random samples works correctly.
4088 * We cover 2 different cases: a sample where the ratio of sample points
4089 * to total points is greater than the threshold for using Algorithm, and
4090 * one where the ratio is small enough (< 1/13) for Algorithm D to be used.
4092 * For each, we also cover the case with and without a specified RNG.
4095 // Small sample/source ratio, no specified RNG.
4096 uint i = 0;
4097 foreach (e; randomSample(randomCover(a), 5))
4099 ++i;
4101 assert(i == 5);
4103 // Small sample/source ratio, specified RNG.
4104 i = 0;
4105 foreach (e; randomSample(randomCover(a), 5, rng))
4107 ++i;
4109 assert(i == 5);
4111 // Large sample/source ratio, no specified RNG.
4112 i = 0;
4113 foreach (e; randomSample(TestInputRange(), 123, 123_456))
4115 ++i;
4117 assert(i == 123);
4119 // Large sample/source ratio, specified RNG.
4120 i = 0;
4121 foreach (e; randomSample(TestInputRange(), 123, 123_456, rng))
4123 ++i;
4125 assert(i == 123);
4127 /* Sample/source ratio large enough to start with Algorithm D,
4128 * small enough to switch to Algorithm A.
4130 i = 0;
4131 foreach (e; randomSample(TestInputRange(), 10, 131))
4133 ++i;
4135 assert(i == 10);
4138 // Test that the .index property works correctly
4140 auto sample1 = randomSample(TestInputRange(), 654, 654_321);
4141 for (; !sample1.empty; sample1.popFront())
4143 assert(sample1.front == sample1.index);
4146 auto sample2 = randomSample(TestInputRange(), 654, 654_321, rng);
4147 for (; !sample2.empty; sample2.popFront())
4149 assert(sample2.front == sample2.index);
4152 /* Check that it also works if .index is called before .front.
4153 * See: https://issues.dlang.org/show_bug.cgi?id=10322
4155 auto sample3 = randomSample(TestInputRange(), 654, 654_321);
4156 for (; !sample3.empty; sample3.popFront())
4158 assert(sample3.index == sample3.front);
4161 auto sample4 = randomSample(TestInputRange(), 654, 654_321, rng);
4162 for (; !sample4.empty; sample4.popFront())
4164 assert(sample4.index == sample4.front);
4168 /* Test behaviour if .popFront() is called before sample is read.
4169 * This is a rough-and-ready check that the statistical properties
4170 * are in the ballpark -- not a proper validation of statistical
4171 * quality! This incidentally also checks for reference-type
4172 * initialization bugs, as the foreach () loop will operate on a
4173 * copy of the popFronted (and hence initialized) sample.
4176 size_t count0, count1, count99;
4177 foreach (_; 0 .. 50_000)
4179 auto sample = randomSample(iota(100), 5, &rng);
4180 sample.popFront();
4181 foreach (s; sample)
4183 if (s == 0)
4185 ++count0;
4187 else if (s == 1)
4189 ++count1;
4191 else if (s == 99)
4193 ++count99;
4197 /* Statistical assumptions here: this is a sequential sampling process
4198 * so (i) 0 can only be the first sample point, so _can't_ be in the
4199 * remainder of the sample after .popFront() is called. (ii) By similar
4200 * token, 1 can only be in the remainder if it's the 2nd point of the
4201 * whole sample, and hence if 0 was the first; probability of 0 being
4202 * first and 1 second is 5/100 * 4/99 (thank you, Algorithm S:-) and
4203 * so the mean count of 1 should be about 202. Finally, 99 can only
4204 * be the _last_ sample point to be picked, so its probability of
4205 * inclusion should be independent of the .popFront() and it should
4206 * occur with frequency 5/100, hence its count should be about 5000.
4207 * Unfortunately we have to set quite a high tolerance because with
4208 * sample size small enough for unittests to run in reasonable time,
4209 * the variance can be quite high.
4211 assert(count0 == 0);
4212 assert(count1 < 150, text("1: ", count1, " > 150."));
4213 assert(2_200 < count99, text("99: ", count99, " < 2200."));
4214 assert(count99 < 2_800, text("99: ", count99, " > 2800."));
4217 /* Odd corner-cases: RandomSample has 2 constructors that are not called
4218 * by the randomSample() helper functions, but that can be used if the
4219 * constructor is called directly. These cover the case of the user
4220 * specifying input but not input length.
4223 auto input1 = TestInputRange().takeExactly(456_789);
4224 static assert(hasLength!(typeof(input1)));
4225 auto sample1 = RandomSample!(typeof(input1), void)(input1, 789);
4226 static assert(isInputRange!(typeof(sample1)));
4227 static assert(!isForwardRange!(typeof(sample1)));
4228 assert(sample1.length == 789);
4229 assert(sample1._available == 456_789);
4230 uint i = 0;
4231 for (; !sample1.empty; sample1.popFront())
4233 assert(sample1.front == sample1.index);
4234 ++i;
4236 assert(i == 789);
4238 auto input2 = TestInputRange().takeExactly(456_789);
4239 static assert(hasLength!(typeof(input2)));
4240 auto sample2 = RandomSample!(typeof(input2), typeof(rng))(input2, 789, rng);
4241 static assert(isInputRange!(typeof(sample2)));
4242 static assert(!isForwardRange!(typeof(sample2)));
4243 assert(sample2.length == 789);
4244 assert(sample2._available == 456_789);
4245 i = 0;
4246 for (; !sample2.empty; sample2.popFront())
4248 assert(sample2.front == sample2.index);
4249 ++i;
4251 assert(i == 789);
4254 /* Test that the save property works where input is a forward range,
4255 * and RandomSample is using a (forward range) random number generator
4256 * that is not rndGen.
4258 static if (isForwardRange!UniformRNG)
4260 auto sample1 = randomSample(a, 5, rng);
4261 // https://issues.dlang.org/show_bug.cgi?id=15853
4262 auto sample2 = ((const ref typeof(sample1) a) => a.save)(sample1);
4263 assert(sample1.array() == sample2.array());
4266 // https://issues.dlang.org/show_bug.cgi?id=8314
4268 auto sample(RandomGen)(uint seed) { return randomSample(a, 1, RandomGen(seed)).front; }
4270 // Start from 1 because not all RNGs accept 0 as seed.
4271 immutable fst = sample!UniformRNG(1);
4272 uint n = 1;
4273 while (sample!UniformRNG(++n) == fst && n < n.max) {}
4274 assert(n < n.max);
4276 }();