The OpenD Programming Language

1 /++
2 	A random number generator that can work with [std.random] but does not have to.
3 
4 	It is designed to be reasonably good and fast, more for fun than for security.
5 
6 
7 	Authors:
8 		Forked from Herringway's pcg.d file:
9 		https://github.com/Herringway/unexpected/blob/main/pcg/source/unexpected/pcg.d
10 
11 		Modified by Adam D. Ruppe
12 	Copyright:
13 		Original version copyright Herringway, 2023
14 
15 	License: BSL-1.0
16 
17 		Boost Software License - Version 1.0 - August 17th, 2003
18 
19 		Permission is hereby granted, free of charge, to any person or organization
20 		obtaining a copy of the software and accompanying documentation covered by
21 		this license (the "Software") to use, reproduce, display, distribute,
22 		execute, and transmit the Software, and to prepare derivative works of the
23 		Software, and to permit third-parties to whom the Software is furnished to
24 		do so, all subject to the following:
25 
26 		The copyright notices in the Software and this entire statement, including
27 		the above license grant, this restriction and the following disclaimer,
28 		must be included in all copies of the Software, in whole or in part, and
29 		all derivative works of the Software, unless such copies or derivative
30 		works are solely in the form of machine-executable object code generated by
31 		a source language processor.
32 
33 		THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
34 		IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
35 		FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
36 		SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
37 		FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
38 		ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
39 		DEALINGS IN THE SOFTWARE.
40 +/
41 module arsd.random;
42 
43 /+
44 Herringway: adam_d_ruppe: when you get back, there're a few other things you might wanna consider for your std.random
45 Herringway: like seeding with ranges instead of single values (mersenne twister has a looooot of state that needs seeding, and a single value isn't doing to do a very good job)
46 Herringway: as well as providing more sources of data to seed with, ike OS APIs n such
47 +/
48 
49 // desired functions:
50 // https://phobos.dpldocs.info/source/std.random.d.html#L2119
51 
52 /++
53 	Gets a random number from a uniform distribution including min and up to (but not including) max from the reasonable default generator.
54 
55 	History:
56 		Added April 17, 2025
57 +/
58 int uniform(int min, int max) {
59 	return uniform(getReasonableDefaultGenerator(), min, max);
60 }
61 
62 /// ditto
63 int uniform(Rng gen, int min, int max) {
64 	auto f = cast(uint) gen.next;
65 
66 	// FIXME i think this is biased but also meh
67 	return f % (max - min) + min;
68 }
69 
70 /// ditto
71 alias randomInteger = uniform;
72 
73 /+
74 unittest {
75 	import arsd.core;
76 	writeln(uniform(-10, 0));
77 }
78 +/
79 
80 /++
81 	Gets a random number between 0.0 and 1.0, including 0.0, but not including 1.0.
82 
83 	History:
84 		Added April 18, 2025
85 +/
86 float randomFloat() {
87 	return randomFloat(getReasonableDefaultGenerator());
88 }
89 
90 /// ditto
91 float randomFloat(Rng gen) {
92 	auto max = (1 << float.mant_dig) - 1;
93 	float n = uniform(gen, 0, max);
94 	return n / max;
95 }
96 
97 // might do a long uniform and maybe double too? idk
98 
99 /++
100 	Shuffles the contents of the array, in place. Assumes elements can be easily swapped.
101 	(the current implementation is an in-place Fisher-Yates algorithm)
102 
103 	History:
104 		Added April 19, 2025
105 +/
106 void shuffle(T)(T[] array) {
107 	shuffle(getReasonableDefaultGenerator(), array);
108 }
109 
110 /// ditto
111 void shuffle(T)(Rng gen, T[] array) {
112 	assert(array.length < int.max);
113 
114 	foreach(index, item; array) {
115 		auto ridx = uniform(gen, cast(int) index, cast(int) array.length);
116 		if(ridx == index)
117 			continue;
118 		array[index] = array[ridx];
119 		array[ridx] = item;
120 	}
121 }
122 
123 version(arsd_random_unittest)
124 unittest {
125 	auto array = [1,2,3,4,5,6,7,8,9,0];
126 	auto results = new int[](array.length);
127 
128 	foreach(i; 0 .. 1_000_000) {
129 	shuffle(array);
130 	auto searchingFor = 9;
131 	foreach(where, item; array)
132 		if(searchingFor == item)
133 			results[where]++;
134 	}
135 
136 	import arsd.core; writeln(results);
137 }
138 
139 /++
140 	Returns an index of the weights, with the proportional odds given by the weights.
141 
142 	So weightedChoice([1, 2, 1]) is twice as likely to return 1 as it is 0 or 2.
143 
144 	History:
145 		Added April 19, 2025
146 +/
147 int weightedChoice(scope const int[] weights...) {
148 	return weightedChoice(getReasonableDefaultGenerator(), weights);
149 }
150 
151 /// ditto
152 int weightedChoice(Rng gen, scope const int[] weights...) {
153 	int sum = 0;
154 	foreach(weight; weights)
155 		sum += weight;
156 
157 	int val = uniform(gen, 0, sum);
158 
159 	sum = 0;
160 	foreach(idx, weight; weights) {
161 		sum += weight;
162 
163 		if(val < sum)
164 			return cast(int) idx;
165 	}
166 
167 	assert(0);
168 }
169 
170 /++
171 	Pick a random number off the normal (aka gaussian) distribution bell curve.
172 
173 	Parameters:
174 		median = median
175 		standardDeviation = standard deviation
176 		min = minimum value to ever return
177 		max = one above the highest value to ever return; an exclusive endpoint
178 
179 	History:
180 		Added April 18, 2025
181 +/
182 int bellCurve(int median, int standardDeviation, int min = int.min, int max = int.max) {
183 	return bellCurve(getReasonableDefaultGenerator(), median, standardDeviation, min, max);
184 }
185 
186 /// ditto
187 int bellCurve(Rng gen, int median, int standardDeviation, int min = int.min, int max = int.max) {
188 	import core.stdc.math;
189 
190 	auto mag =  standardDeviation * sqrt(-2.0 * log(randomFloat(gen)));
191 	int value = cast(int) (mag * cos(2 * 3.14159268f * randomFloat(gen)) + median);
192 	if(value < min)
193 		value = min;
194 	if(value >= max)
195 		value = max - 1;
196 	return value;
197 }
198 // bimodal distribution?
199 // maybe a pareto distribution too idk tho
200 
201 
202 version(arsd_random_unittest)
203 unittest {
204 	int[21] results;
205 	foreach(i; 0 .. 1_000_00) {
206 		//results[uniform(0, cast(int) results.length)] += 1;
207 		//results[bellCurve(10, 3, 0, cast(int) results.length)] += 1;
208 
209 		results[weightedChoice([0, 2, 1, 0, 2, 6, 0, 6, 6])] += 1;
210 	}
211 	import std.stdio; writeln(results);
212 
213 	// foreach(i; 0 .. 10) writeln(bellCurve(100, 10));
214 }
215 
216 /++
217 	A simple generic interface to a random number generator.
218 +/
219 interface Rng {
220 	/++
221 		Seeds the generator, calling the delegate zero (if it is a true rng), one, or more times to get all the state it needs.
222 	+/
223 	void seed(scope ulong delegate() getEntropy);
224 	/++
225 		Get the next number in the sequence. Some may not actually use all 64 bits of the return type.
226 	+/
227 	ulong next();
228 
229 	/+
230 	/++
231 		Saves a copy of the current generator state to a fresh object.
232 
233 		See_Also:
234 			[saveState], which returns an array of bytes you can save to a file (or whatever)
235 	+/
236 	Rng save() const;
237 	+/
238 }
239 
240 /+
241 interface RestorableRng {
242 	/++
243 		Saves the rng state to an array.
244 
245 		To restore state, you must first construct an object of the same type, then call `restoreState`
246 		on that fresh object. If you get the wrong type, it won't work right (and may or may not throw an exception).
247 	+/
248 	ubyte[] saveState() const;
249 
250 	/// ditto
251 	void restoreState(in ubyte[]);
252 }
253 +/
254 
255 class RngFromRange(R) : Rng {
256 	private R r;
257 
258 	void seed(scope ulong delegate() getEntropy) {
259 		r = R(getEntropy());
260 	}
261 	ulong next() {
262 		auto f = r.front;
263 		r.popFront;
264 		return f;
265 	}
266 	Rng save() const {
267 		auto t = new RngFromRange();
268 		t.r = this.r.save;
269 		return t;
270 	}
271 }
272 
273 alias reasonableDefault = PCG!(uint, ulong, xslrr);
274 
275 /++
276 	Gets a "reasonable default" random number generator, one good enough
277 	for my casual use. This is the object used by the other functions when
278 	you don't explicitly use your own generator.
279 
280 	It will be automatically seeded from the operating system random number
281 	pool if you don't pass one of your own.
282 
283 	History:
284 		Added April 17, 2025
285 +/
286 Rng getReasonableDefaultGenerator(lazy ulong seed) @trusted {
287 	static Rng generator;
288 	if(generator is null) {
289 		generator = new RngFromRange!reasonableDefault();
290 		generator.seed(&seed);
291 	}
292 	return generator;
293 }
294 
295 /// ditto
296 Rng getReasonableDefaultGenerator() {
297 	return getReasonableDefaultGenerator(unpredictableSeed());
298 }
299 
300 private ulong unpredictableSeed() {
301 	ulong r;
302 	osRandom(r);
303 	return r;
304 }
305 
306 version (none) {
307 } else version (linux) {
308 	private bool osRandom(out ulong result) @trusted {
309 		import core.sys.posix.unistd;
310 		import core.sys.posix.fcntl;
311 		int fd = open("/dev/urandom", O_RDONLY);
312 		if(fd == -1)
313 			return false;
314 		auto ret = read(fd, &result, typeof(result).sizeof);
315 		if(ret != typeof(result).sizeof) {
316 			close(fd);
317 			return false;
318 		}
319 		close(fd);
320 		return true;
321 	}
322 
323 } else version (Windows) {
324 	pragma(lib, "Bcrypt.lib");
325 
326 	private bool osRandom(out ulong result) @trusted {
327 		import core.sys.windows.windef : PUCHAR, ULONG;
328 		import core.sys.windows.ntdef : NT_SUCCESS;
329 		import core.sys.windows.bcrypt : BCryptGenRandom, BCRYPT_USE_SYSTEM_PREFERRED_RNG;
330 
331 		const gotRandom = BCryptGenRandom(
332 			null,
333 			cast(PUCHAR) &result,
334 			ULONG(typeof(result).sizeof),
335 			BCRYPT_USE_SYSTEM_PREFERRED_RNG,
336 		);
337 
338 		return NT_SUCCESS(gotRandom);
339 	}
340 } else version (all) {
341 	private bool osRandom(out ulong result) @trusted {
342 		import std.random;
343 		result = std.random.unpredictableSeed;
344 		return false;
345 	}
346 }
347 
348 private V rotr(V)(V value, uint r) {
349 	return cast(V)(value >> r | value << (-r & (V.sizeof * 8 - 1)));
350 }
351 
352 private int log2(int d) {
353 	assert(__ctfe);
354 	if(d == 8) return 3;
355 	if(d == 16) return 4;
356 	if(d == 32) return 5;
357 	if(d == 64) return 6;
358 	if(d == 128) return 7;
359 	assert(0);
360 }
361 
362 struct PCGConsts(X, I) {
363 	enum spareBits = (I.sizeof - X.sizeof) * 8;
364 	enum wantedOpBits = log2(X.sizeof * 8);
365 	struct xshrr {
366 		enum opBits = spareBits >= wantedOpBits ? wantedOpBits : spareBits;
367 		enum amplifier = wantedOpBits - opBits;
368 		enum xShift = (opBits + X.sizeof * 8) / 2;
369 		enum mask = (1 << opBits) - 1;
370 		enum bottomSpare = spareBits - opBits;
371 	}
372 	struct xshrs {
373 		// there must be a simpler way to express this
374 		static if (spareBits - 5 >= 64) {
375 			enum opBits = 5;
376 		} else static if (spareBits - 4 >= 32) {
377 			enum opBits = 4;
378 		} else static if (spareBits - 3 >= 16) {
379 			enum opBits = 3;
380 		} else static if (spareBits - 2 >= 4) {
381 			enum opBits = 2;
382 		} else static if (spareBits - 1 >= 1) {
383 			enum opBits = 1;
384 		} else {
385 			enum opBits = 0;
386 		}
387 		enum xShift = opBits + ((X.sizeof * 8) + mask) / 2;
388 		enum mask = (1 << opBits) - 1;
389 		enum bottomSpare = spareBits - opBits;
390 	}
391 	struct xsh {
392 		enum topSpare = 0;
393 		enum bottomSpare = spareBits - topSpare;
394 		enum xShift = (topSpare + X.sizeof * 8) / 2;
395 	}
396 	struct xsl {
397 		enum topSpare = spareBits;
398 		enum bottomSpare = spareBits - topSpare;
399 		enum xShift = (topSpare + X.sizeof * 8) / 2;
400 	}
401 	struct rxs {
402 		enum shift = (I.sizeof - X.sizeof) * 8;
403 		// there must be a simpler way to express this
404 		static if (shift > 64 + 8) {
405 			enum rShiftAmount = I.sizeof - 6;
406 			enum rShiftMask = 63;
407 		} else static if (shift > 32 + 4) {
408 			enum rShiftAmount = I.sizeof - 5;
409 			enum rShiftMask = 31;
410 		} else static if (shift > 16 + 2) {
411 			enum rShiftAmount = I.sizeof - 4;
412 			enum rShiftMask = 15;
413 		} else static if (shift > 8 + 1) {
414 			enum rShiftAmount = I.sizeof - 3;
415 			enum rShiftMask = 7;
416 		} else static if (shift > 4 + 1) {
417 			enum rShiftAmount = I.sizeof - 2;
418 			enum rShiftMask = 3;
419 		} else static if (shift > 2 + 1) {
420 			enum rShiftAmount = I.sizeof - 1;
421 			enum rShiftMask = 1;
422 		} else {
423 			enum rShiftAmount = 0;
424 			enum rShiftMask = 0;
425 		}
426 		enum extraShift = (X.sizeof - shift)/2;
427 	}
428 	struct rxsm {
429 		enum opBits = log2(X.sizeof * 8) - 1;
430 		enum shift = (I.sizeof - X.sizeof) * 8;
431 		enum mask = (1 << opBits) - 1;
432 	}
433 	struct xslrr {
434 		enum opBits = spareBits >= wantedOpBits ? wantedOpBits : spareBits;
435 		enum amplifier = wantedOpBits - opBits;
436 		enum mask = (1 << opBits) - 1;
437 		enum topSpare = spareBits;
438 		enum bottomSpare = spareBits - topSpare;
439 		enum xShift = (topSpare + X.sizeof * 8) / 2;
440 	}
441 }
442 
443 private X xorshift(X, I)(I tmp, uint amt1, uint amt2) {
444 	tmp ^= tmp >> amt1;
445 	return cast(X)(tmp >> amt2);
446 }
447 
448 /// XSH RR (xorshift high, random rotate) - decent performance, slightly better results
449 private X xshrr(X, I)(const I state) {
450 	alias constants = PCGConsts!(X, I).xshrr;
451 	static if (constants.opBits > 0) {
452 		auto rot = (state >> (I.sizeof * 8 - constants.opBits)) & constants.mask;
453 	} else {
454 		enum rot = 0;
455 	}
456 	uint amprot = cast(uint)((rot << constants.amplifier) & constants.mask);
457 	I tmp = state;
458 	return rotr(xorshift!X(tmp, constants.xShift, constants.bottomSpare), amprot);
459 }
460 
461 /// XSH RS (xorshift high, random shift) - decent performance
462 private X xshrs(X, I)(const I state) {
463 	alias constants = PCGConsts!(X, I).xshrs;
464 	static if (constants.opBits > 0) {
465 		uint rshift = (state >> (I.sizeof * 8 - constants.opBits)) & constants.mask;
466 	} else {
467 		uint rshift = 0;
468 	}
469 	I tmp = state;
470 	return xorshift!X(tmp, constants.xShift, cast(uint)(constants.bottomSpare - constants.mask + rshift));
471 }
472 
473 /// XSH (fixed xorshift, high) - don't use this for anything smaller than ulong
474 private X xsh(X, I)(const I state) {
475 	alias constants = PCGConsts!(X, I).xsh;
476 
477 	I tmp = state;
478 	return xorshift!X(tmp, constants.xShift, constants.bottomSpare);
479 }
480 
481 /// XSL (fixed xorshift, low) - don't use this for anything smaller than ulong
482 private X xsl(X, I)(const I state) {
483 	alias constants = PCGConsts!(X, I).xsl;
484 
485 	I tmp = state;
486 	return xorshift!X(tmp, constants.xShift, constants.bottomSpare);
487 }
488 
489 /// RXS (random xorshift)
490 private X rxs(X, I)(const I state) {
491 	alias constants = PCGConsts!(X, I).rxs;
492 	uint rshift = (state >> constants.rShiftAmount) & constants.rShiftMask;
493 	I tmp = state;
494 	return xorshift!X(tmp, cast(uint)(constants.shift + constants.extraShift - rshift), rshift);
495 }
496 
497 /++
498 	RXS M XS (random xorshift, multiply, fixed xorshift)
499 	This has better statistical properties, but supposedly performs worse. This
500 	was not reproducible, however.
501 +/
502 private X rxsmxs(X, I)(const I state) {
503 	X result = rxsm!X(state);
504 	result ^= result >> ((2 * X.sizeof * 8 + 2) / 3);
505 	return result;
506 }
507 
508 /// RXS M (random xorshift, multiply)
509 private X rxsm(X, I)(const I state) {
510 	alias constants = PCGConsts!(X, I).rxsm;
511 	I tmp = state;
512 	static if (constants.opBits > 0) {
513 		uint rshift = (tmp >> (I.sizeof * 8 - constants.opBits)) & constants.mask;
514 	} else {
515 		uint rshift = 0;
516 	}
517 	tmp ^= tmp >> (constants.opBits + rshift);
518 	tmp *= PCGMMultiplier!I;
519 	return cast(X)(tmp >> constants.shift);
520 }
521 
522 /// DXSM (double xorshift, multiply) - newer, better performance for types 2x the size of the largest type the cpu can handle
523 private X dxsm(X, I)(const I state) {
524 	static assert(X.sizeof <= I.sizeof/2, "Output type must be half the size of the state type.");
525 	X hi = cast(X)(state >> ((I.sizeof - X.sizeof) * 8));
526 	X lo = cast(X)state;
527 
528 	lo |= 1;
529 	hi ^= hi >> (X.sizeof * 8 / 2);
530 	hi *= PCGMMultiplier!I;
531 	hi ^= hi >> (3*(X.sizeof * 8 / 4));
532 	hi *= lo;
533 	return hi;
534 }
535 /// XSL RR (fixed xorshift, random rotate) - better performance for types 2x the size of the largest type the cpu can handle
536 private X xslrr(X, I)(const I state) {
537 	alias constants = PCGConsts!(X, I).xslrr;
538 
539 	I tmp = state;
540 	static if (constants.opBits > 0) {
541 		uint rot = (tmp >> (I.sizeof * 8 - constants.opBits)) & constants.mask;
542 	} else {
543 		uint rot = 0;
544 	}
545 	uint amprot = (rot << constants.amplifier) & constants.mask;
546 	return rotr(xorshift!X(tmp, constants.xShift, constants.bottomSpare), amprot);
547 }
548 
549 struct PCG(T, S, alias func, S multiplier = DefaultPCGMultiplier!S, S increment = DefaultPCGIncrement!S) {
550 	private S state;
551 
552 	this(S val) @safe pure nothrow @nogc {
553 		seed(val);
554 	}
555 	void seed(S val) @safe pure nothrow @nogc {
556 		state = cast(S)(val + increment);
557 		popFront();
558 	}
559 	void popFront() @safe pure nothrow @nogc {
560 		state = cast(S)(state * multiplier + increment);
561 	}
562 	T front() const @safe pure nothrow @nogc @property {
563 		return func!T(state);
564 	}
565 	typeof(this) save() @safe pure nothrow @nogc const {
566 		return this;
567 	}
568 	enum bool empty = false;
569 	enum bool isUniformRandom = true;
570 	enum T min = T.min;
571 	enum T max = T.max;
572 	const(S) toSiryulType()() const @safe {
573 		return state;
574 	}
575 	static PCG fromSiryulType()(S val) @safe {
576 		PCG result;
577 		result.state = val;
578 		return result;
579 	}
580 }
581 
582 template DefaultPCGMultiplier(T) {
583 	static if (is(T == ubyte)) {
584 		enum DefaultPCGMultiplier = 141;
585 	} else static if (is(T == ushort)) {
586 		enum DefaultPCGMultiplier = 12829;
587 	} else static if (is(T == uint)) {
588 		enum DefaultPCGMultiplier = 747796405;
589 	} else static if (is(T == ulong)) {
590 		enum DefaultPCGMultiplier = 6364136223846793005;
591 	} else static if (is(T == ucent)) {
592 		//enum DefaultPCGMultiplier = 47026247687942121848144207491837523525;
593 	}
594 }
595 
596 template DefaultPCGIncrement(T) {
597 	static if (is(T == ubyte)) {
598 		enum DefaultPCGIncrement = 77;
599 	} else static if (is(T == ushort)) {
600 		enum DefaultPCGIncrement = 47989;
601 	} else static if (is(T == uint)) {
602 		enum DefaultPCGIncrement = 2891336453;
603 	} else static if (is(T == ulong)) {
604 		enum DefaultPCGIncrement = 1442695040888963407;
605 	} else static if (is(T == ucent)) {
606 		//enum DefaultPCGIncrement = 117397592171526113268558934119004209487;
607 	}
608 }
609 
610 private template PCGMMultiplier(T) {
611 	static if (is(T : ubyte)) {
612 		enum PCGMMultiplier = 217;
613 	} else static if (is(T : ushort)) {
614 		enum PCGMMultiplier = 62169;
615 	} else static if (is(T : uint)) {
616 		enum PCGMMultiplier = 277803737;
617 	} else static if (is(T : ulong)) {
618 		enum PCGMMultiplier = 12605985483714917081;
619 	//} else static if (is(T == ucent)) {
620 		//enum PCGMMultiplier = 327738287884841127335028083622016905945;
621 	}
622 }
623 
624 version(arsd_random_unittest) {
625 	private alias AliasSeq(T...) = T;
626 
627 	alias SupportedTypes = AliasSeq!(ubyte, ushort, uint, ulong);
628 	alias SupportedFunctions = AliasSeq!(xshrr, xshrs, xsh, xsl, rxs, rxsmxs, rxsm, xslrr);
629 
630 	static foreach (ResultType; SupportedTypes) {
631 		static foreach (StateType; SupportedTypes) {
632 			static if (StateType.sizeof >= ResultType.sizeof) {
633 				static foreach (Function; SupportedFunctions) {
634 					mixin("alias PCG", int(StateType.sizeof * 8), int(ResultType.sizeof * 8), __traits(identifier, Function), " = PCG!(ResultType, StateType, Function);");
635 				}
636 			}
637 		}
638 	}
639 	alias PCG6432dxsm = PCG!(uint, ulong, dxsm);
640 
641 	@safe unittest {
642 		import std.algorithm : reduce;
643 		import std.datetime.stopwatch : benchmark;
644 		import std.math : pow, sqrt;
645 		import std.random : isSeedable, Mt19937, uniform, uniform01, unpredictableSeed;
646 		import std.stdio : writefln, writeln;
647 		auto seed = unpredictableSeed;
648 
649 		void testRNG(RNG, string name)(uint seed) {
650 			static if (isSeedable!(RNG, uint)) {
651 				auto rng = RNG(seed);
652 			} else static if (isSeedable!(RNG, ushort)) {
653 				auto rng = RNG(cast(ushort)seed);
654 			} else static if (isSeedable!(RNG, ubyte)) {
655 				auto rng = RNG(cast(ubyte)seed);
656 			}
657 			writefln!"--%s--"(name);
658 			double total = 0;
659 			ulong[ubyte] distribution;
660 			void test() {
661 				total += uniform01(rng);
662 				distribution.require(uniform!ubyte(rng), 0)++;
663 			}
664 			auto result = benchmark!(test)(1000000)[0];
665 			writefln!"Benchmark completed in %s"(result);
666 			writeln(total);
667 			double avg = reduce!((a, b) => a + b / distribution.length)(0.0f, distribution);
668 			auto var = reduce!((a, b) => a + pow(b - avg, 2) / distribution.length)(0.0f, distribution);
669 			auto sd = sqrt(var);
670 			writefln!"Average: %s, Standard deviation: %s"(avg, sd);
671 		}
672 
673 		testRNG!(PCG168xshrr, "PCG168xshrr")(seed);
674 		testRNG!(PCG3216xshrr, "PCG3216xshrr")(seed);
675 		testRNG!(PCG6432xslrr, "PCG6432xslrr")(seed);
676 		testRNG!(PCG648rxsmxs, "PCG648rxsmxs")(seed);
677 		testRNG!(PCG6432dxsm, "PCG6432dxsm")(seed);
678 		testRNG!(Mt19937, "Mt19937")(seed);
679 		testRNG!(reasonableDefault, "reasonableDefault")(seed);
680 	}
681 }