mirror of https://github.com/adamdruppe/arsd.git
new module for good'n'quick random numbers from herringway
This commit is contained in:
Adam D. Ruppe
parent
09f209ed6d
commit
dd0d79fb51
|
|
@ -0,0 +1,395 @@
|
|||
/++
|
||||
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);
|
||||
}
|
||||
}
|
||||
Loading…
Reference in New Issue