arsd/random.d

/++
	A random number generator that can work with [std.random] but does not have to.


	Authors:
		Forked from Herringway's pcg.d file:
		https://github.com/Herringway/unexpected/blob/main/pcg/source/unexpected/pcg.d

		Modified by Adam D. Ruppe
	Copyright:
		Original version copyright Herringway, 2023

	License: BSL-1.0

		Boost Software License - Version 1.0 - August 17th, 2003

		Permission is hereby granted, free of charge, to any person or organization
		obtaining a copy of the software and accompanying documentation covered by
		this license (the "Software") to use, reproduce, display, distribute,
		execute, and transmit the Software, and to prepare derivative works of the
		Software, and to permit third-parties to whom the Software is furnished to
		do so, all subject to the following:

		The copyright notices in the Software and this entire statement, including
		the above license grant, this restriction and the following disclaimer,
		must be included in all copies of the Software, in whole or in part, and
		all derivative works of the Software, unless such copies or derivative
		works are solely in the form of machine-executable object code generated by
		a source language processor.

		THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
		IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
		FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
		SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
		FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
		ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
		DEALINGS IN THE SOFTWARE.
+/
module arsd.random;

/+
Herringway: adam_d_ruppe: when you get back, there're a few other things you might wanna consider for your std.random
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)
Herringway: as well as providing more sources of data to seed with, ike OS APIs n such
+/


alias reasonableDefault = PCG!(uint, ulong, xslrr);

// desired functions:
// https://phobos.dpldocs.info/source/std.random.d.html#L2119
int uniform(int min, int max) { return 0; }
// might do a long uniform and maybe float? but idk divide that mebbe
void shuffle(T)(T[] array) {} // fisher-yates algorithm
int weightedChoice(scope const int[] weights...) { return 0; } // std.random.dice
// the normal / gaussian distribution
int bellCurve(int median, int standardDeviation) { return 0; }
// bimodal distribution?
// maybe a pareto distribution too idk tho


private V rotr(V)(V value, uint r) {
	return cast(V)(value >> r | value << (-r & (V.sizeof * 8 - 1)));
}

private int log2(int d) {
	assert(__ctfe);
	if(d == 8) return 3;
	if(d == 16) return 4;
	if(d == 32) return 5;
	if(d == 64) return 6;
	if(d == 128) return 7;
	assert(0);
}

struct PCGConsts(X, I) {
	enum spareBits = (I.sizeof - X.sizeof) * 8;
	enum wantedOpBits = log2(X.sizeof * 8);
	struct xshrr {
		enum opBits = spareBits >= wantedOpBits ? wantedOpBits : spareBits;
		enum amplifier = wantedOpBits - opBits;
		enum xShift = (opBits + X.sizeof * 8) / 2;
		enum mask = (1 << opBits) - 1;
		enum bottomSpare = spareBits - opBits;
	}
	struct xshrs {
		// there must be a simpler way to express this
		static if (spareBits - 5 >= 64) {
			enum opBits = 5;
		} else static if (spareBits - 4 >= 32) {
			enum opBits = 4;
		} else static if (spareBits - 3 >= 16) {
			enum opBits = 3;
		} else static if (spareBits - 2 >= 4) {
			enum opBits = 2;
		} else static if (spareBits - 1 >= 1) {
			enum opBits = 1;
		} else {
			enum opBits = 0;
		}
		enum xShift = opBits + ((X.sizeof * 8) + mask) / 2;
		enum mask = (1 << opBits) - 1;
		enum bottomSpare = spareBits - opBits;
	}
	struct xsh {
		enum topSpare = 0;
		enum bottomSpare = spareBits - topSpare;
		enum xShift = (topSpare + X.sizeof * 8) / 2;
	}
	struct xsl {
		enum topSpare = spareBits;
		enum bottomSpare = spareBits - topSpare;
		enum xShift = (topSpare + X.sizeof * 8) / 2;
	}
	struct rxs {
		enum shift = (I.sizeof - X.sizeof) * 8;
		// there must be a simpler way to express this
		static if (shift > 64 + 8) {
			enum rShiftAmount = I.sizeof - 6;
			enum rShiftMask = 63;
		} else static if (shift > 32 + 4) {
			enum rShiftAmount = I.sizeof - 5;
			enum rShiftMask = 31;
		} else static if (shift > 16 + 2) {
			enum rShiftAmount = I.sizeof - 4;
			enum rShiftMask = 15;
		} else static if (shift > 8 + 1) {
			enum rShiftAmount = I.sizeof - 3;
			enum rShiftMask = 7;
		} else static if (shift > 4 + 1) {
			enum rShiftAmount = I.sizeof - 2;
			enum rShiftMask = 3;
		} else static if (shift > 2 + 1) {
			enum rShiftAmount = I.sizeof - 1;
			enum rShiftMask = 1;
		} else {
			enum rShiftAmount = 0;
			enum rShiftMask = 0;
		}
		enum extraShift = (X.sizeof - shift)/2;
	}
	struct rxsm {
		enum opBits = log2(X.sizeof * 8) - 1;
		enum shift = (I.sizeof - X.sizeof) * 8;
		enum mask = (1 << opBits) - 1;
	}
	struct xslrr {
		enum opBits = spareBits >= wantedOpBits ? wantedOpBits : spareBits;
		enum amplifier = wantedOpBits - opBits;
		enum mask = (1 << opBits) - 1;
		enum topSpare = spareBits;
		enum bottomSpare = spareBits - topSpare;
		enum xShift = (topSpare + X.sizeof * 8) / 2;
	}
}

private X xorshift(X, I)(I tmp, uint amt1, uint amt2) {
	tmp ^= tmp >> amt1;
	return cast(X)(tmp >> amt2);
}

/// XSH RR (xorshift high, random rotate) - decent performance, slightly better results
private X xshrr(X, I)(const I state) {
	alias constants = PCGConsts!(X, I).xshrr;
	static if (constants.opBits > 0) {
		auto rot = (state >> (I.sizeof * 8 - constants.opBits)) & constants.mask;
	} else {
		enum rot = 0;
	}
	uint amprot = cast(uint)((rot << constants.amplifier) & constants.mask);
	I tmp = state;
	return rotr(xorshift!X(tmp, constants.xShift, constants.bottomSpare), amprot);
}

/// XSH RS (xorshift high, random shift) - decent performance
private X xshrs(X, I)(const I state) {
	alias constants = PCGConsts!(X, I).xshrs;
	static if (constants.opBits > 0) {
		uint rshift = (state >> (I.sizeof * 8 - constants.opBits)) & constants.mask;
	} else {
		uint rshift = 0;
	}
	I tmp = state;
	return xorshift!X(tmp, constants.xShift, cast(uint)(constants.bottomSpare - constants.mask + rshift));
}

/// XSH (fixed xorshift, high) - don't use this for anything smaller than ulong
private X xsh(X, I)(const I state) {
	alias constants = PCGConsts!(X, I).xsh;

	I tmp = state;
	return xorshift!X(tmp, constants.xShift, constants.bottomSpare);
}

/// XSL (fixed xorshift, low) - don't use this for anything smaller than ulong
private X xsl(X, I)(const I state) {
	alias constants = PCGConsts!(X, I).xsl;

	I tmp = state;
	return xorshift!X(tmp, constants.xShift, constants.bottomSpare);
}

/// RXS (random xorshift)
private X rxs(X, I)(const I state) {
	alias constants = PCGConsts!(X, I).rxs;
	uint rshift = (state >> constants.rShiftAmount) & constants.rShiftMask;
	I tmp = state;
	return xorshift!X(tmp, cast(uint)(constants.shift + constants.extraShift - rshift), rshift);
}

/++
	RXS M XS (random xorshift, multiply, fixed xorshift)
	This has better statistical properties, but supposedly performs worse. This
	was not reproducible, however.
+/
private X rxsmxs(X, I)(const I state) {
	X result = rxsm!X(state);
	result ^= result >> ((2 * X.sizeof * 8 + 2) / 3);
	return result;
}

/// RXS M (random xorshift, multiply)
private X rxsm(X, I)(const I state) {
	alias constants = PCGConsts!(X, I).rxsm;
	I tmp = state;
	static if (constants.opBits > 0) {
		uint rshift = (tmp >> (I.sizeof * 8 - constants.opBits)) & constants.mask;
	} else {
		uint rshift = 0;
	}
	tmp ^= tmp >> (constants.opBits + rshift);
	tmp *= PCGMMultiplier!I;
	return cast(X)(tmp >> constants.shift);
}

/// DXSM (double xorshift, multiply) - newer, better performance for types 2x the size of the largest type the cpu can handle
private X dxsm(X, I)(const I state) {
	static assert(X.sizeof <= I.sizeof/2, "Output type must be half the size of the state type.");
	X hi = cast(X)(state >> ((I.sizeof - X.sizeof) * 8));
	X lo = cast(X)state;

	lo |= 1;
	hi ^= hi >> (X.sizeof * 8 / 2);
	hi *= PCGMMultiplier!I;
	hi ^= hi >> (3*(X.sizeof * 8 / 4));
	hi *= lo;
	return hi;
}
/// XSL RR (fixed xorshift, random rotate) - better performance for types 2x the size of the largest type the cpu can handle
private X xslrr(X, I)(const I state) {
	alias constants = PCGConsts!(X, I).xslrr;

	I tmp = state;
	static if (constants.opBits > 0) {
		uint rot = (tmp >> (I.sizeof * 8 - constants.opBits)) & constants.mask;
	} else {
		uint rot = 0;
	}
	uint amprot = (rot << constants.amplifier) & constants.mask;
	return rotr(xorshift!X(tmp, constants.xShift, constants.bottomSpare), amprot);
}

struct PCG(T, S, alias func, S multiplier = DefaultPCGMultiplier!S, S increment = DefaultPCGIncrement!S) {
	private S state;

	this(S val) @safe pure nothrow @nogc {
		seed(val);
	}
	void seed(S val) @safe pure nothrow @nogc {
		state = cast(S)(val + increment);
		popFront();
	}
	void popFront() @safe pure nothrow @nogc {
		state = cast(S)(state * multiplier + increment);
	}
	T front() const @safe pure nothrow @nogc @property {
		return func!T(state);
	}
	typeof(this) save() @safe pure nothrow @nogc {
		return this;
	}
	enum bool empty = false;
	enum bool isUniformRandom = true;
	enum T min = T.min;
	enum T max = T.max;
	const(S) toSiryulType()() const @safe {
		return state;
	}
	static PCG fromSiryulType()(S val) @safe {
		PCG result;
		result.state = val;
		return result;
	}
}

template DefaultPCGMultiplier(T) {
	static if (is(T == ubyte)) {
		enum DefaultPCGMultiplier = 141;
	} else static if (is(T == ushort)) {
		enum DefaultPCGMultiplier = 12829;
	} else static if (is(T == uint)) {
		enum DefaultPCGMultiplier = 747796405;
	} else static if (is(T == ulong)) {
		enum DefaultPCGMultiplier = 6364136223846793005;
	} else static if (is(T == ucent)) {
		//enum DefaultPCGMultiplier = 47026247687942121848144207491837523525;
	}
}

template DefaultPCGIncrement(T) {
	static if (is(T == ubyte)) {
		enum DefaultPCGIncrement = 77;
	} else static if (is(T == ushort)) {
		enum DefaultPCGIncrement = 47989;
	} else static if (is(T == uint)) {
		enum DefaultPCGIncrement = 2891336453;
	} else static if (is(T == ulong)) {
		enum DefaultPCGIncrement = 1442695040888963407;
	} else static if (is(T == ucent)) {
		//enum DefaultPCGIncrement = 117397592171526113268558934119004209487;
	}
}

private template PCGMMultiplier(T) {
	static if (is(T : ubyte)) {
		enum PCGMMultiplier = 217;
	} else static if (is(T : ushort)) {
		enum PCGMMultiplier = 62169;
	} else static if (is(T : uint)) {
		enum PCGMMultiplier = 277803737;
	} else static if (is(T : ulong)) {
		enum PCGMMultiplier = 12605985483714917081;
	//} else static if (is(T == ucent)) {
		//enum PCGMMultiplier = 327738287884841127335028083622016905945;
	}
}

version(arsd_random_unittest) {
	private alias AliasSeq(T...) = T;

	alias SupportedTypes = AliasSeq!(ubyte, ushort, uint, ulong);
	alias SupportedFunctions = AliasSeq!(xshrr, xshrs, xsh, xsl, rxs, rxsmxs, rxsm, xslrr);

	static foreach (ResultType; SupportedTypes) {
		static foreach (StateType; SupportedTypes) {
			static if (StateType.sizeof >= ResultType.sizeof) {
				static foreach (Function; SupportedFunctions) {
					mixin("alias PCG", int(StateType.sizeof * 8), int(ResultType.sizeof * 8), __traits(identifier, Function), " = PCG!(ResultType, StateType, Function);");
				}
			}
		}
	}
	alias PCG6432dxsm = PCG!(uint, ulong, dxsm);

	@safe unittest {
		import std.algorithm : reduce;
		import std.datetime.stopwatch : benchmark;
		import std.math : pow, sqrt;
		import std.random : isSeedable, Mt19937, uniform, uniform01, unpredictableSeed;
		import std.stdio : writefln, writeln;
		auto seed = unpredictableSeed;

		void testRNG(RNG, string name)(uint seed) {
			static if (isSeedable!(RNG, uint)) {
				auto rng = RNG(seed);
			} else static if (isSeedable!(RNG, ushort)) {
				auto rng = RNG(cast(ushort)seed);
			} else static if (isSeedable!(RNG, ubyte)) {
				auto rng = RNG(cast(ubyte)seed);
			}
			writefln!"--%s--"(name);
			double total = 0;
			ulong[ubyte] distribution;
			void test() {
				total += uniform01(rng);
				distribution.require(uniform!ubyte(rng), 0)++;
			}
			auto result = benchmark!(test)(1000000)[0];
			writefln!"Benchmark completed in %s"(result);
			writeln(total);
			double avg = reduce!((a, b) => a + b / distribution.length)(0.0f, distribution);
			auto var = reduce!((a, b) => a + pow(b - avg, 2) / distribution.length)(0.0f, distribution);
			auto sd = sqrt(var);
			writefln!"Average: %s, Standard deviation: %s"(avg, sd);
		}

		testRNG!(PCG168xshrr, "PCG168xshrr")(seed);
		testRNG!(PCG3216xshrr, "PCG3216xshrr")(seed);
		testRNG!(PCG6432xslrr, "PCG6432xslrr")(seed);
		testRNG!(PCG648rxsmxs, "PCG648rxsmxs")(seed);
		testRNG!(PCG6432dxsm, "PCG6432dxsm")(seed);
		testRNG!(Mt19937, "Mt19937")(seed);
		testRNG!(reasonableDefault, "reasonableDefault")(seed);
	}
}