The OpenD Programming Language

1 /++
2  + Permuted Congruential Generator (PCG)
3  +
4  + Implemented as per the C++ version of PCG, $(HTTP _pcg-random.org).
5  +
6  + Paper available $(HTTP _pcg-random.org/paper.html)
7  +
8  + Author:  Melissa O'Neill (C++). D translation Nicholas Wilson.
9  +
10  + PCG Random Number Generation for C++
11  +
12  + Copyright 2014 Melissa O'Neill <oneill@pcg-random.org>
13  +
14  + Licensed under the Apache License, Version 2.0 (the "License");
15  + you may not use this file except in compliance with the License.
16  + You may obtain a copy of the License at
17  +
18  +     $(HTTP www.apache.org/licenses/LICENSE-2.0)
19  +
20  + Unless required by applicable law or agreed to in writing, software
21  + distributed under the License is distributed on an "AS IS" BASIS,
22  + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
23  + See the License for the specific language governing permissions and
24  + limitations under the License.
25  +
26  + For additional information about the PCG random number generation scheme,
27  + including its license and other licensing options, visit
28  +
29  +     $(HTTP _pcg-random.org)
30  +/
31 module mir.random.engine.pcg;
32 
33 import mir.random.engine;
34 import std.traits : ReturnType, TemplateArgsOf;
35 
36 @safe:
37 nothrow:
38 @nogc:
39 
40  ///
41 @nogc nothrow pure @safe version(mir_random_test) unittest
42 {
43     import mir.random /+: isSaturatedRandomEngine, rand+/;
44     import mir.random.engine.pcg : pcg32;
45     import mir.math.common: fabs;
46 
47     static assert(isRandomEngine!pcg32);
48     static assert(isSaturatedRandomEngine!pcg32);
49     auto gen = pcg32(1234u);//Seed with constant.
50     assert(gen.rand!double.fabs == 0x1.e5fe29e2942a8p-1);//Generate number from 0 inclusive to 1 exclusive.
51     assert(gen.rand!ulong == 13957536084079755083UL);
52 }
53 
54 /// Select the above mixin templates.
55 enum stream_t
56 {
57     /// Increment is cast(size_t) &RNG.
58     unique,
59     /// Increment is 0.
60     none,
61     /// Increment is the default increment.
62     oneseq,
63     /// Increment is runtime setable and defaults to the default (same as oneseq)
64     specific,
65 }
66 
67 /// 32-bit output PCGs with 64 bits of state.
68 alias pcg32        = PermutedCongruentialEngine!(xsh_rr!(uint,ulong),stream_t.specific,true);
69 /// ditto
70 alias pcg32_unique = PermutedCongruentialEngine!(xsh_rr!(uint,ulong),stream_t.unique,true);
71 /// ditto
72 alias pcg32_oneseq = PermutedCongruentialEngine!(xsh_rr!(uint,ulong),stream_t.oneseq,true);
73 /// ditto
74 alias pcg32_fast   = PermutedCongruentialEngine!(xsh_rr!(uint,ulong),stream_t.none,true);
75 
76 static if (__traits(compiles, ucent.max))
77 {
78     /// 64-bit output PCGs with 128 bits of state. Requires `ucent` type.
79     alias pcg64        = PermutedCongruentialEngine!(xsh_rr!(ulong,ucent),stream_t.specific,true);
80     ///
81     alias pcg64_unique = PermutedCongruentialEngine!(xsh_rr!(ulong,ucent),stream_t.unique,true);
82     ///
83     alias pcg64_oneseq = PermutedCongruentialEngine!(xsh_rr!(ulong,ucent),stream_t.oneseq,true);
84     ///
85     alias pcg64_fast   = PermutedCongruentialEngine!(xsh_rr!(ulong,ucent),stream_t.none,true);
86 }
87 
88 /// PCGs with n bits output and n bits of state.
89 alias pcg8_once_insecure  = PermutedCongruentialEngine!(rxs_m_xs_forward!(ubyte ,ubyte ),stream_t.specific,true);
90 /// ditto
91 alias pcg16_once_insecure = PermutedCongruentialEngine!(rxs_m_xs_forward!(ushort,ushort),stream_t.specific,true);
92 /// ditto
93 alias pcg32_once_insecure = PermutedCongruentialEngine!(rxs_m_xs_forward!(uint  ,uint  ),stream_t.specific,true);
94 /// ditto
95 alias pcg64_once_insecure = PermutedCongruentialEngine!(rxs_m_xs_forward!(ulong,ulong  ),stream_t.specific,true);
96 //alias pcg128_once_insecure = PermutedCongruentialEngine!(rxs_m_xs_forward!(ucent,ucent,stream_t.specific,true);
97 
98 /// As above but the increment is not dynamically setable.
99 alias pcg8_oneseq_once_insecure = PermutedCongruentialEngine!(rxs_m_xs_forward!(ubyte ,ubyte ),stream_t.oneseq,true);
100 /// ditto
101 alias pcg16_oneseq_once_insecure = PermutedCongruentialEngine!(rxs_m_xs_forward!(ushort,ushort),stream_t.oneseq,true);
102 /// ditto
103 alias pcg32_oneseq_once_insecure = PermutedCongruentialEngine!(rxs_m_xs_forward!(uint  ,uint  ),stream_t.oneseq,true);
104 /// ditto
105 alias pcg64_oneseq_once_insecure = PermutedCongruentialEngine!(rxs_m_xs_forward!(ulong,ulong  ),stream_t.oneseq,true);
106 /// ditto
107 /// Requires `ucent` type.
108 static if (__traits(compiles, ucent.max))
109 alias pcg128_oneseq_once_insecure = PermutedCongruentialEngine!(rxs_m_xs_forward!(ucent,ucent  ),stream_t.specific,true);
110 
111 /++
112  + The PermutedCongruentialEngine:
113  + Params:
114  +  output = should be one of the above functions.
115  +      Controls the output permutation of the state.
116  +  streamType = one of unique, none, oneseq, specific.
117  +      Controls the Increment of the LCG portion of the PCG.
118  +  output_previous = if true then the pre-advance version (increasing instruction-level parallelism)
119  +      if false then use the post-advance version (reducing register pressure)
120  +  mult_ = optionally set the multiplier for the LCG.
121  +/
122 struct PermutedCongruentialEngine(alias output,        // Output function
123                                   stream_t streamType, // The stream type
124                                   bool output_previous,
125                                   mult_...) if (mult_.length <= 1)
126 {
127     ///
128     enum isRandomEngine = true;
129 
130     ///
131     alias Uint  = TemplateArgsOf!output[1];
132 
133     static if (mult_.length == 0)
134         enum mult = default_multiplier!Uint;
135     else
136     {
137         static assert(is(typeof(mult_[0]) == Uint),
138             "The specified multiplier must be the state type of the output function");
139         enum mult = mult_[0];
140     }
141         
142     @disable this(this);
143     @disable this();
144     static if (streamType == stream_t.none)
145         mixin no_stream!Uint;
146     else static if (streamType == stream_t.unique)
147         mixin unique_stream!Uint;
148     else static if (streamType == stream_t.specific)
149         mixin specific_stream!Uint;
150     else static if (streamType == stream_t.oneseq)
151         mixin oneseq_stream!Uint;
152     else
153         static assert(0);
154 
155     ///
156     Uint state;
157     
158     ///
159     enum period_pow2 = Uint.sizeof*8 - 2*is_mcg;
160 
161     ///
162     enum max = (ReturnType!output).max;
163 
164 private:
165 
166     static if (__traits(compiles, { enum e = mult + increment; }))
167     {
168         static Uint bump()(Uint state_)
169         {
170             return cast(Uint)(state_ * mult + increment);
171         }
172     }
173     else
174     {
175         Uint bump()(Uint state_)
176         {
177             return cast(Uint)(state_ * mult + increment);
178         }
179     }
180 
181     Uint base_generate()()
182     {
183         return state = bump(state);
184     }
185 
186     Uint base_generate0()()
187     {
188         Uint old_state = state;
189         state = bump(state);
190         return old_state;
191     }
192 
193 public:
194     static if (can_specify_stream)
195     ///
196     this()(Uint seed, Uint stream_ = default_increment_unset_stream!Uint)
197     {
198         state = bump(cast(Uint)(seed + increment));
199         set_stream(stream_);
200     }
201     else
202     ///
203     this()(Uint seed)
204     {
205         static if (is_mcg)
206             state = seed | 3u;
207         else
208             state = bump(cast(Uint)(seed + increment));
209     }
210 
211     ///
212     ReturnType!output opCall()()
213     {
214         static if(output_previous)
215             return output(base_generate0());
216         else
217             return output(base_generate());
218     }
219 
220     /++
221     Skip forward in the random sequence in $(BIGOH log n) time.
222     Even though delta is an unsigned integer, we can pass a
223     signed integer to go backwards, it just goes "the long way round".
224     +/
225     void skip()(Uint delta)
226     {
227         // The method used here is based on Brown, "Random Number Generation
228         // with Arbitrary Stride,", Transactions of the American Nuclear
229         // Society (Nov. 1994).  The algorithm is very similar to fast
230         // exponentiation.
231         //
232         // Even though delta is an unsigned integer, we can pass a
233         // signed integer to go backwards, it just goes "the long way round".
234         
235         Uint acc_mult = 1, acc_plus = 0;
236         Uint cur_plus = increment, cur_mult = mult;
237         while (delta > 0)
238         {
239             if (delta & 1)
240             {
241                 acc_mult *= cur_mult;
242                 acc_plus = cast(Uint)(acc_plus * cur_mult + cur_plus);
243             }
244             cur_plus  *= cur_mult + 1;
245             cur_mult *= cur_mult;
246             delta >>= 1;
247         }
248         state = cast(Uint)(acc_mult*state + acc_plus);
249     }
250 
251     static if (output_previous)
252     {
253         /++
254         Compatibility with $(LINK2 https://dlang.org/phobos/std_random.html#.isUniformRNG,
255         Phobos library methods). Presents this RNG as an InputRange.
256         Only available if `output_previous == true`.
257 
258         The reason that this is enabled when `output_previous == true` is because
259         `front` can be implemented without additional cost.
260 
261         `save` is implemented if `streamType` is not [stream_t.unique].
262 
263         This struct disables its default copy constructor and so will only work with
264         Phobos functions that "do the right thing" and take RNGs by reference and
265         do not accidentally make implicit copies.
266         +/
267         enum bool isUniformRandom = true;
268         /// ditto
269         enum ReturnType!output min = (ReturnType!output).min;
270         /// ditto
271         enum bool empty = false;
272         /// ditto
273         @property ReturnType!output front()() const { return output(state); }
274         /// ditto
275         void popFront()() { state = bump(state); }
276         /// ditto
277         void seed()(Uint seed) { this.__ctor(seed); }
278         /// ditto
279         static if (streamType != stream_t.unique)
280         @property typeof(this) save()() const
281         {
282             typeof(return) result = void;
283             foreach (i, e; this.tupleof)
284                 result.tupleof[i] = e;
285             return result;
286         }
287     }
288 }
289 
290 @nogc nothrow pure @safe version(mir_random_test) unittest
291 {
292     //Test that the default generators (all having output_previous==true)
293     //can be used as Phobos-style randoms.
294     import std.meta: AliasSeq;
295     import std.random: isSeedable, isPhobosUniformRNG = isUniformRNG;
296     import std.range.primitives: isForwardRange;
297     foreach(RNG; AliasSeq!(pcg32, pcg32_oneseq, pcg32_fast,
298                            pcg32_once_insecure, pcg32_oneseq_once_insecure,
299                            pcg64_once_insecure, pcg64_oneseq_once_insecure))
300     {
301         static assert(isPhobosUniformRNG!(RNG, typeof(RNG.max)));
302         static assert(isSeedable!(RNG, RNG.Uint));
303         static assert(isForwardRange!RNG);
304         auto gen1 = RNG(1);
305         auto gen2 = RNG(2);
306         auto gen3 = gen1.save;
307         gen2.seed(1);
308         assert(gen1 == gen2);
309         immutable a = gen1.front;
310         gen1.popFront();
311         assert(a == gen2());
312         assert(gen1.front == gen2());
313         assert(a == gen3());
314         assert(gen1.front == gen3());
315     }
316 
317     foreach(RNG; AliasSeq!(pcg32_unique))
318     {
319         static assert(isPhobosUniformRNG!(RNG, typeof(RNG.max)));
320         static assert(isSeedable!(RNG, RNG.Uint));
321     }
322 }
323 
324 // Default multiplier to use for the LCG portion of the PCG
325 private template default_multiplier(Uint)
326 {
327     static if (is(Uint == ubyte))
328         enum ubyte default_multiplier = 141u;
329     else static if (is(Uint == ushort))
330         enum ushort default_multiplier = 12829u;
331     else static if (is(Uint == uint))
332         enum uint default_multiplier = 747796405u;
333     else static if (is(Uint == ulong))
334         enum ulong default_multiplier = 6364136223846793005u;
335     else static if (is(ucent) && is(Uint == ucent))
336         mixin("enum ucent default_multiplier = 0x2360ED051FC65DA44385DF649FCCF645;");
337     else
338         static assert(0);
339 }
340 
341 // Default increment to use for the LCG portion of the PCG
342 private template default_increment(Uint)
343 {
344     static if (is(Uint == ubyte))
345         enum ubyte default_increment = 77u;
346     else static if (is(Uint == ushort))
347         enum ushort default_increment = 47989u;
348     else static if (is(Uint == uint))
349         enum uint default_increment = 2891336453u;
350     else static if (is(Uint == ulong))
351         enum ulong default_increment = 1442695040888963407u;
352     else static if (is(ucent) && is(Uint == ucent))
353         mixin("enum ucent default_increment = 0x5851F42D4C957F2D14057B7EF767814F;");
354     else
355         static assert(0);
356 }
357 
358 private template mcg_multiplier(Uint)
359 {
360     static if (is(Uint == ubyte))
361         enum ubyte mcg_multiplier = 217u;
362     else static if (is(Uint == ushort))
363         enum ushort mcg_multiplier = 62169u;
364     else static if (is(Uint == uint))
365         enum uint mcg_multiplier = 277803737u;
366     else static if (is(Uint == ulong))
367         enum ulong mcg_multiplier = 12605985483714917081u;
368     else static if (is(ucent) && is(Uint == ucent))
369         mixin("enum ucent mcg_multiplier = 0x6BC8F622C397699CAEF17502108EF2D9;");
370     else
371         static assert(0);
372 }
373 
374 private template mcg_unmultiplier(Uint)
375 {
376     static if (is(Uint == ubyte))
377         enum ubyte mcg_unmultiplier = 105u;
378     else static if (is(Uint == ushort))
379         enum ushort mcg_unmultiplier = 28009u;
380     else static if (is(Uint == uint))
381         enum uint mcg_unmultiplier = 2897767785u;
382     else static if (is(Uint == ulong))
383         enum ulong mcg_unmultiplier = 15009553638781119849u;
384     else static if (is(ucent) && is(Uint == ucent))
385         mixin("enum ucent mcg_unmultiplier = 0xC827645E182BC965D04CA582ACB86D69;");
386     else
387         static assert(0);
388 }
389 
390 private template default_increment_unset_stream(Uint)
391 {
392     enum default_increment_unset_stream = (default_increment!Uint & ~1) >> 1;
393 }
394 
395 /// Increment for LCG portion of the PCG is the address of the RNG
396 mixin template unique_stream(Uint)
397 {
398     ///
399     enum is_mcg = false;
400     ///
401     @property Uint increment()() const @trusted
402     {
403         return cast(Uint) (&this) | 1;
404     }
405     ///
406     Uint stream()()
407     {
408         return increment >> 1;
409     }
410     ///
411     enum can_specify_stream = false;
412     ///
413     enum size_t streams_pow2 = Uint.sizeof < size_t.sizeof ? Uint.sizeof : size_t.sizeof - 1u;
414 }
415 
416 @nogc nothrow pure @system unittest
417 {
418     pcg32_unique gen = pcg32_unique(1);
419     void* address = &gen;
420     assert(gen.increment == (1 | cast(ulong) address));
421 }
422 
423 
424 /// Increment is 0. The LCG portion of the PCG is an MCG.
425 mixin template no_stream(Uint)
426 {
427     ///
428     enum is_mcg = true;
429     ///
430     enum Uint increment = 0;
431     
432     ///
433     enum can_specify_stream = false;
434     ///
435     enum size_t streams_pow2 = 0;
436 }
437 
438 /// Increment of the LCG portion of the PCG is default_increment.
439 mixin template oneseq_stream(Uint)
440 {
441     ///
442     enum is_mcg = false;
443     ///
444     enum Uint increment = default_increment!Uint;
445     ///
446     enum can_specify_stream = false;
447     ///
448     enum size_t streams_pow2 = 0;
449 }
450 
451 /// The increment is dynamically settable and defaults to default_increment!T.
452 mixin template specific_stream(Uint)
453 {
454     ///
455     enum is_mcg = false;
456     ///
457     Uint inc_ = default_increment!Uint;
458     ///
459     alias increment = inc_;
460     ///
461     enum can_specify_stream = true;
462     ///
463     void set_stream()(Uint u)
464     {
465         inc_ = cast(Uint)((u << 1) | 1);
466     }
467     ///
468     enum size_t streams_pow2 = size_t.sizeof*8 -1u;
469 }
470 
471 /++
472  + XorShifts are invertible, but they are someting of a pain to invert.
473  + This function backs them out.  It's used by the whacky "inside out"
474  + generator defined later.
475  +/
476 Uint unxorshift(Uint)(Uint x, size_t bits, size_t shift)
477 {
478     if (2*shift >= bits) {
479         return x ^ (x >> shift);
480     }
481     Uint lowmask1 = (itype(1U) << (bits - shift*2)) - 1;
482     Uint highmask1 = ~lowmask1;
483     Uint top1 = x;
484     Uint bottom1 = x & lowmask1;
485     top1 ^= top1 >> shift;
486     top1 &= highmask1;
487     x = top1 | bottom1;
488     Uint lowmask2 = (itype(1U) << (bits - shift)) - 1;
489     Uint bottom2 = x & lowmask2;
490     bottom2 = unxorshift(bottom2, bits - shift, shift);
491     bottom2 &= lowmask1;
492     return top1 | bottom2;
493 }
494 
495 
496 /++
497  + OUTPUT FUNCTIONS.
498  +
499  + These are the core of the PCG generation scheme.  They specify how to
500  + turn the base LCG's internal state into the output value of the final
501  + generator.
502  +
503  + All of the classes have code that is written to allow it to be applied
504  + at *arbitrary* bit sizes.
505  +/
506 
507 /++
508  + XSH RS -- high xorshift, followed by a random shift
509  +
510  + Fast.  A good performer.
511  +/
512 O xsh_rs(O, Uint)(Uint state)
513 {
514     enum bits        = Uint.sizeof * 8;
515     enum xtypebits   = O.sizeof * 8;
516     enum sparebits   = bits - xtypebits;
517     enum opbits = sparebits-5 >= 64 ? 5
518                 : sparebits-4 >= 32 ? 4
519                 : sparebits-3 >= 16 ? 3
520                 : sparebits-2 >= 4  ? 2
521                 : sparebits-1 >= 1  ? 1
522                 :                     0;
523     enum mask           = (1 << opbits) - 1;
524     enum maxrandshift   = mask;
525     enum topspare       = opbits;
526     enum bottomspare    = sparebits - topspare;
527     enum xshift         = topspare + (xtypebits+maxrandshift)/2;
528     
529     auto rshift = opbits ? size_t(state >> (bits - opbits)) & mask : 0;
530     state ^= state >> xshift;
531     O result = O(s >> (bottomspare - maxrandshift + rshift));
532     return result;
533 }
534 /++
535  + XSH RR -- high xorshift, followed by a random rotate
536  +
537  + Fast.  A good performer.  Slightly better statistically than XSH RS.
538  +/
539 O xsh_rr(O, Uint)(Uint state)
540 {
541     enum bits        = Uint.sizeof * 8;
542     enum xtypebits   = O.sizeof * 8;
543     enum sparebits   = bits - xtypebits;
544     enum wantedopbits =   xtypebits >= 128 ? 7
545                         : xtypebits >=  64 ? 6
546                         : xtypebits >=  32 ? 5
547                         : xtypebits >=  16 ? 4
548                         :                    3;
549     enum opbits = sparebits >= wantedopbits ? wantedopbits : sparebits;
550     enum amplifier = wantedopbits - opbits;
551     enum mask = (1 << opbits) - 1;
552     enum topspare    = opbits;
553     enum bottomspare = sparebits - topspare;
554     enum xshift      = (topspare + xtypebits)/2;
555     
556     auto rot = opbits ? size_t(state >> (bits - opbits)) & mask : 0;
557     auto amprot = (rot << amplifier) & mask;
558     state ^= state >> xshift;
559     O result = cast(O)(state >> bottomspare);
560     import core.bitop: ror;
561     result = ror(result, cast(uint)amprot);
562     return result;
563 }
564 /++
565  + RXS -- random xorshift
566  +/
567 O rxs(O, Uint)(Uint state)
568 {
569     enum bits        = Uint.sizeof * 8;
570     enum xtypebits   = O.sizeof * 8;
571     enum shift       = bits - xtypebits;
572     enum extrashift  = (xtypebits - shift)/2;
573     enum rshift = shift > 64+8 ? (s >> (bits - 6)) & 63
574                   : shift > 32+4 ? (s >> (bits - 5)) & 31
575                   : shift > 16+2 ? (s >> (bits - 4)) & 15
576                   : shift >  8+1 ? (s >> (bits - 3)) & 7
577                   : shift >  4+1 ? (s >> (bits - 2)) & 3
578                   : shift >  2+1 ? (s >> (bits - 1)) & 1
579                   : 0;
580     state ^= state >> (shift + extrashift - rshift);
581     O result = state >> rshift;
582     return result;
583 }
584 /++
585  + RXS M XS -- random xorshift, mcg multiply, fixed xorshift
586  +
587  + The most statistically powerful generator, but all those steps
588  + make it slower than some of the others.  We give it the rottenest jobs.
589  +
590  + Because it's usually used in contexts where the state type and the
591  + result type are the same, it is a permutation and is thus invertible.
592  + We thus provide a function to invert it.  This function is used to
593  + for the "inside out" generator used by the extended generator.
594  +/
595 O rxs_m_xs_forward(O, Uint)(Uint state)
596     if(is(O == Uint))
597 {
598     enum bits        = Uint.sizeof * 8;
599     enum xtypebits   = O.sizeof * 8;
600     enum opbits = xtypebits >= 128 ? 6
601                 : xtypebits >=  64 ? 5
602                 : xtypebits >=  32 ? 4
603                 : xtypebits >=  16 ? 3
604                 :                    2;
605     enum shift = bits - xtypebits;
606     enum mask = (1 << opbits) - 1;
607     size_t rshift = opbits ? (state >> (bits - opbits)) & mask : 0;
608     state ^= state >> (opbits + rshift);
609     state *= mcg_multiplier!Uint;
610     O result = state >> shift;
611     result ^= result >> ((2U*xtypebits+2U)/3U);
612     return result;
613 }
614 /// ditto
615 O rxs_m_xs_reverse(O, Uint)(Uint state)
616     if(is(O == Uint))
617 {
618     enum bits        = Uint.sizeof * 8;
619     enum opbits = bits >= 128 ? 6
620                 : bits >=  64 ? 5
621                 : bits >=  32 ? 4
622                 : bits >=  16 ? 3
623                 :               2;
624     enum mask = (1 << opbits) - 1;
625     
626     state = unxorshift(state, bits, (2U*bits+2U)/3U);
627     
628     state *= mcg_unmultiplier!Uint;
629     
630     auto rshift = opbits ? (state >> (bits - opbits)) & mask : 0;
631     state = unxorshift(s, bits, opbits + rshift);
632     
633     return s;
634 }
635 /++
636  + XSL RR -- fixed xorshift (to low bits), random rotate
637  +
638  + Useful for 128-bit types that are split across two CPU registers.
639  +/
640 O xsl_rr(O, Uint)(Uint state)
641 {
642     enum bits        = Uint.sizeof * 8;
643     enum xtypebits   = O.sizeof * 8;
644     enum sparebits = bits - xtypebits;
645     enum wantedopbits =   xtypebits >= 128 ? 7
646                         : xtypebits >=  64 ? 6
647                         : xtypebits >=  32 ? 5
648                         : xtypebits >=  16 ? 4
649                         :                    3;
650     enum opbits = sparebits >= wantedopbits ? wantedopbits : sparebits;
651     enum amplifier = wantedopbits - opbits;
652     enum mask = (1 << opbits) - 1;
653     enum topspare = sparebits;
654     enum bottomspare = sparebits - topspare;
655     enum xshift = (topspare + xtypebits) / 2;
656     
657     auto rot = opbits ? size_t(state >> (bits - opbits)) & mask : 0;
658     auto amprot = (rot << amplifier) & mask;
659     state ^= state >> xshift;
660     O result = state >> bottomspare;
661     result = rotr(result, amprot);
662     return result;
663 }
664 
665 private template half_size(Uint)
666 {
667     static if (is(Uint == ucent))
668         alias half_size = ulong;
669     else static if (is(Uint == ulong))
670         alias half_size = uint;
671     else static if (is(Uint == uint))
672         alias half_size = ushort;
673     else static if (is(Uint == ushort))
674         alias half_size = ubyte;
675     else
676         static assert(0);
677 }
678 
679 /++
680  + XSL RR RR -- fixed xorshift (to low bits), random rotate (both parts)
681  +
682  + Useful for 128-bit types that are split across two CPU registers.
683  + If you really want an invertible 128-bit RNG, I guess this is the one.
684  +/
685 
686 O xsl_rr_rr(O, Uint)(Uint state)
687     if(is(O == Uint))
688 {
689     alias H = half_size!Uint;
690     enum htypebits = H.sizeof * 8;
691     enum bits      = Uint.sizeof * 8;
692     enum sparebits = bits - htypebits;
693     enum wantedopbits =   htypebits >= 128 ? 7
694                         : htypebits >=  64 ? 6
695                         : htypebits >=  32 ? 5
696                         : htypebits >=  16 ? 4
697                         :                    3;
698     enum opbits = sparebits >= wantedopbits ? wantedopbits : sparebits;
699     enum amplifier = wantedopbits - opbits;
700     enum mask = (1 << opbits) - 1;
701     enum topspare = sparebits;
702     enum xshift = (topspare + htypebits) / 2;
703     
704     auto rot = opbits ? size_t(s >> (bits - opbits)) & mask : 0;
705     auto amprot = (rot << amplifier) & mask;
706     state ^= state >> xshift;
707     H lowbits = cast(H)s;
708     lowbits = rotr(lowbits, amprot);
709     H highbits = cast(H)(s >> topspare);
710     auto rot2 = lowbits & mask;
711     auto amprot2 = (rot2 << amplifier) & mask;
712     highbits = rotr(highbits, amprot2);
713     return (O(highbits) << topspare) ^ O(lowbits);
714     
715 }
716 
717 /++
718  + XSH -- fixed xorshift (to high bits)
719  +
720  + Not available at 64-bits or less.
721  +/
722 
723 O xsh(O, Uint)(Uint state) if(Uint.sizeof > 8)
724 {
725     enum bits        = Uint.sizeof * 8;
726     enum xtypebits   = O.sizeof * 8;
727     enum sparebits = bits - xtypebits;
728     enum topspare = 0;
729     enum bottomspare = sparebits - topspare;
730     enum xshift = (topspare + xtypebits) / 2;
731     
732     state ^= state >> xshift;
733     O result = state >> bottomspare;
734     return result;
735 }
736 
737 /++
738  + XSL -- fixed xorshift (to low bits)
739  +
740  + Not available at 64-bits or less.
741  +/
742 
743 O xsl(O, Uint)(Uint state) if(Uint.sizeof > 8)
744 {
745     enum bits        = Uint.sizeof * 8;
746     enum xtypebits   = O.sizeof * 8;
747     enum sparebits = bits - xtypebits;
748     enum topspare = sparebits;
749     enum bottomspare = sparebits - topspare;
750     enum xshift = (topspare + xtypebits) / 2;
751     
752     state ^= state >> xshift;
753     O result = state >> bottomspare;
754     return result;
755 }
756 
757 private alias AliasSeq(T...) = T;
758 
759 @safe version(mir_random_test) unittest
760 {
761     
762     foreach(RNG; AliasSeq!(pcg32,pcg32_unique,pcg32_oneseq,pcg32_fast,
763                            pcg8_once_insecure,pcg16_once_insecure,pcg32_once_insecure,pcg64_once_insecure,
764                            pcg8_oneseq_once_insecure,pcg16_oneseq_once_insecure,pcg32_oneseq_once_insecure,
765                            pcg64_oneseq_once_insecure))
766     {
767         static assert(isSaturatedRandomEngine!RNG);
768         auto gen = RNG(cast(RNG.Uint)unpredictableSeed);
769         gen();
770     }
771 }
772 @safe version(mir_random_test) unittest
773 {
774     auto gen = pcg32(0x198c8585);
775     gen.skip(1000);
776     assert(gen()== 0xd187a760);
777     auto gen2 = pcg32(0x198c8585);
778     
779     foreach(_; 0 .. 1000)
780         gen2();
781     assert(gen2() == 0xd187a760);
782     assert(gen() == gen2());
783 }
784 
785 @nogc nothrow pure @safe unittest
786 {
787     foreach (ShouldHaveStaticBump; AliasSeq!(pcg32_oneseq, pcg32_fast, pcg32_oneseq_once_insecure))
788         static assert (__traits(compiles, { enum e = ShouldHaveStaticBump.bump(1u); }));
789 
790     foreach (ShouldLackStaticBump; AliasSeq!(pcg32, pcg32_unique, pcg32_once_insecure))
791         static assert (!__traits(compiles, { enum e = ShouldLackStaticBump.bump(1u); }));
792 }