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