arsd/random.d

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

	It is designed to be reasonably good and fast, more for fun than for security.


	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
+/

// desired functions:
// https://phobos.dpldocs.info/source/std.random.d.html#L2119

/++
	Gets a random number from a uniform distribution including min and up to (but not including) max from the reasonable default generator.

	History:
		Added April 17, 2025
+/
int uniform(int min, int max) {
	return uniform(getReasonableDefaultGenerator(), min, max);
}

/// ditto
int uniform(Rng gen, int min, int max) {
	auto f = cast(uint) gen.next;

	// FIXME i think this is biased but also meh
	return f % (max - min) + min;
}

/// ditto
alias randomInteger = uniform;

/+
unittest {
	import arsd.core;
	writeln(uniform(-10, 0));
}
+/

/++
	Gets a random number between 0.0 and 1.0, including 0.0, but not including 1.0.

	History:
		Added April 18, 2025
+/
float randomFloat() {
	return randomFloat(getReasonableDefaultGenerator());
}

/// ditto
float randomFloat(Rng gen) {
	auto max = (1 << float.mant_dig) - 1;
	float n = uniform(gen, 0, max);
	return n / max;
}

// might do a long uniform and maybe double too? idk

/++
	Shuffles the contents of the array, in place. Assumes elements can be easily swapped.
	(the current implementation is an in-place Fisher-Yates algorithm)

	History:
		Added April 19, 2025
+/
void shuffle(T)(T[] array) {
	shuffle(getReasonableDefaultGenerator(), array);
}

/// ditto
void shuffle(T)(Rng gen, T[] array) {
	assert(array.length < int.max);

	foreach(index, item; array) {
		auto ridx = uniform(gen, cast(int) index, cast(int) array.length);
		if(ridx == index)
			continue;
		array[index] = array[ridx];
		array[ridx] = item;
	}
}

version(arsd_random_unittest)
unittest {
	auto array = [1,2,3,4,5,6,7,8,9,0];
	auto results = new int[](array.length);

	foreach(i; 0 .. 1_000_000) {
	shuffle(array);
	auto searchingFor = 9;
	foreach(where, item; array)
		if(searchingFor == item)
			results[where]++;
	}

	import arsd.core; writeln(results);
}

/++
	Returns an index of the weights, with the proportional odds given by the weights.

	So weightedChoice([1, 2, 1]) is twice as likely to return 1 as it is 0 or 2.

	History:
		Added April 19, 2025
+/
int weightedChoice(scope const int[] weights...) {
	return weightedChoice(getReasonableDefaultGenerator(), weights);
}

/// ditto
int weightedChoice(Rng gen, scope const int[] weights...) {
	int sum = 0;
	foreach(weight; weights)
		sum += weight;

	int val = uniform(gen, 0, sum);

	sum = 0;
	foreach(idx, weight; weights) {
		sum += weight;

		if(val < sum)
			return cast(int) idx;
	}

	assert(0);
}

/++
	Pick a random number off the normal (aka gaussian) distribution bell curve.

	Parameters:
		median = median
		standardDeviation = standard deviation
		min = minimum value to ever return
		max = one above the highest value to ever return; an exclusive endpoint

	History:
		Added April 18, 2025
+/
int bellCurve(int median, int standardDeviation, int min = int.min, int max = int.max) {
	return bellCurve(getReasonableDefaultGenerator(), median, standardDeviation, min, max);
}

/// ditto
int bellCurve(Rng gen, int median, int standardDeviation, int min = int.min, int max = int.max) {
	import core.stdc.math;

	auto mag =  standardDeviation * sqrt(-2.0 * log(randomFloat(gen)));
	int value = cast(int) (mag * cos(2 * 3.14159268f * randomFloat(gen)) + median);
	if(value < min)
		value = min;
	if(value >= max)
		value = max - 1;
	return value;
}
// bimodal distribution?
// maybe a pareto distribution too idk tho


version(arsd_random_unittest)
unittest {
	int[21] results;
	foreach(i; 0 .. 1_000_00) {
		//results[uniform(0, cast(int) results.length)] += 1;
		//results[bellCurve(10, 3, 0, cast(int) results.length)] += 1;

		results[weightedChoice([0, 2, 1, 0, 2, 6, 0, 6, 6])] += 1;
	}
	import std.stdio; writeln(results);

	// foreach(i; 0 .. 10) writeln(bellCurve(100, 10));
}

/++
	A simple generic interface to a random number generator.
+/
interface Rng {
	/++
		Seeds the generator, calling the delegate zero (if it is a true rng), one, or more times to get all the state it needs.
	+/
	void seed(scope ulong delegate() getEntropy);
	/++
		Get the next number in the sequence. Some may not actually use all 64 bits of the return type.
	+/
	ulong next();

	/+
	/++
		Saves a copy of the current generator state to a fresh object.

		See_Also:
			[saveState], which returns an array of bytes you can save to a file (or whatever)
	+/
	Rng save() const;
	+/
}

/+
interface RestorableRng {
	/++
		Saves the rng state to an array.

		To restore state, you must first construct an object of the same type, then call `restoreState`
		on that fresh object. If you get the wrong type, it won't work right (and may or may not throw an exception).
	+/
	ubyte[] saveState() const;

	/// ditto
	void restoreState(in ubyte[]);
}
+/

class RngFromRange(R) : Rng {
	private R r;

	void seed(scope ulong delegate() getEntropy) {
		r = R(getEntropy());
	}
	ulong next() {
		auto f = r.front;
		r.popFront;
		return f;
	}
	Rng save() const {
		auto t = new RngFromRange();
		t.r = this.r.save;
		return t;
	}
}

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

/++
	Gets a "reasonable default" random number generator, one good enough
	for my casual use. This is the object used by the other functions when
	you don't explicitly use your own generator.

	It will be automatically seeded from the operating system random number
	pool if you don't pass one of your own.

	History:
		Added April 17, 2025
+/
Rng getReasonableDefaultGenerator(lazy ulong seed) @trusted {
	static Rng generator;
	if(generator is null) {
		generator = new RngFromRange!reasonableDefault();
		generator.seed(&seed);
	}
	return generator;
}

/// ditto
Rng getReasonableDefaultGenerator() {
	return getReasonableDefaultGenerator(unpredictableSeed());
}

private ulong unpredictableSeed() {
	ulong r;
	osRandom(r);
	return r;
}

version (none) {
} else version (linux) {
	private bool osRandom(out ulong result) @trusted {
		import core.sys.posix.unistd;
		import core.sys.posix.fcntl;
		int fd = open("/dev/urandom", O_RDONLY);
		if(fd == -1)
			return false;
		auto ret = read(fd, &result, typeof(result).sizeof);
		if(ret != typeof(result).sizeof) {
			close(fd);
			return false;
		}
		close(fd);
		return true;
	}

} else version (Windows) {
	pragma(lib, "Bcrypt.lib");

	private bool osRandom(out ulong result) @trusted {
		import core.sys.windows.windef : PUCHAR, ULONG;
		import core.sys.windows.ntdef : NT_SUCCESS;
		import core.sys.windows.bcrypt : BCryptGenRandom, BCRYPT_USE_SYSTEM_PREFERRED_RNG;

		const gotRandom = BCryptGenRandom(
			null,
			cast(PUCHAR) &result,
			ULONG(typeof(result).sizeof),
			BCRYPT_USE_SYSTEM_PREFERRED_RNG,
		);

		return NT_SUCCESS(gotRandom);
	}
} else version (all) {
	private bool osRandom(out ulong result) @trusted {
		import std.random;
		result = std.random.unpredictableSeed;
		return false;
	}
}

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 const {
		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);
	}
}