The OpenD Programming Language

1 /++
2 	A random number generator that can work with [std.random] but does not have to.
3 
4 
5 	Authors:
6 		Forked from Herringway's pcg.d file:
7 		https://github.com/Herringway/unexpected/blob/main/pcg/source/unexpected/pcg.d
8 
9 		Modified by Adam D. Ruppe
10 	Copyright:
11 		Original version copyright Herringway, 2023
12 
13 	License: BSL-1.0
14 
15 		Boost Software License - Version 1.0 - August 17th, 2003
16 
17 		Permission is hereby granted, free of charge, to any person or organization
18 		obtaining a copy of the software and accompanying documentation covered by
19 		this license (the "Software") to use, reproduce, display, distribute,
20 		execute, and transmit the Software, and to prepare derivative works of the
21 		Software, and to permit third-parties to whom the Software is furnished to
22 		do so, all subject to the following:
23 
24 		The copyright notices in the Software and this entire statement, including
25 		the above license grant, this restriction and the following disclaimer,
26 		must be included in all copies of the Software, in whole or in part, and
27 		all derivative works of the Software, unless such copies or derivative
28 		works are solely in the form of machine-executable object code generated by
29 		a source language processor.
30 
31 		THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
32 		IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
33 		FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
34 		SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
35 		FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
36 		ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
37 		DEALINGS IN THE SOFTWARE.
38 +/
39 module arsd.random;
40 
41 /+
42 Herringway: adam_d_ruppe: when you get back, there're a few other things you might wanna consider for your std.random
43 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)
44 Herringway: as well as providing more sources of data to seed with, ike OS APIs n such
45 +/
46 
47 
48 alias reasonableDefault = PCG!(uint, ulong, xslrr);
49 
50 // desired functions:
51 // https://phobos.dpldocs.info/source/std.random.d.html#L2119
52 int uniform(int min, int max) { return 0; }
53 // might do a long uniform and maybe float? but idk divide that mebbe
54 void shuffle(T)(T[] array) {} // fisher-yates algorithm
55 int weightedChoice(scope const int[] weights...) { return 0; } // std.random.dice
56 // the normal / gaussian distribution
57 int bellCurve(int median, int standardDeviation) { return 0; }
58 // bimodal distribution?
59 // maybe a pareto distribution too idk tho
60 
61 
62 private V rotr(V)(V value, uint r) {
63 	return cast(V)(value >> r | value << (-r & (V.sizeof * 8 - 1)));
64 }
65 
66 private int log2(int d) {
67 	assert(__ctfe);
68 	if(d == 8) return 3;
69 	if(d == 16) return 4;
70 	if(d == 32) return 5;
71 	if(d == 64) return 6;
72 	if(d == 128) return 7;
73 	assert(0);
74 }
75 
76 struct PCGConsts(X, I) {
77 	enum spareBits = (I.sizeof - X.sizeof) * 8;
78 	enum wantedOpBits = log2(X.sizeof * 8);
79 	struct xshrr {
80 		enum opBits = spareBits >= wantedOpBits ? wantedOpBits : spareBits;
81 		enum amplifier = wantedOpBits - opBits;
82 		enum xShift = (opBits + X.sizeof * 8) / 2;
83 		enum mask = (1 << opBits) - 1;
84 		enum bottomSpare = spareBits - opBits;
85 	}
86 	struct xshrs {
87 		// there must be a simpler way to express this
88 		static if (spareBits - 5 >= 64) {
89 			enum opBits = 5;
90 		} else static if (spareBits - 4 >= 32) {
91 			enum opBits = 4;
92 		} else static if (spareBits - 3 >= 16) {
93 			enum opBits = 3;
94 		} else static if (spareBits - 2 >= 4) {
95 			enum opBits = 2;
96 		} else static if (spareBits - 1 >= 1) {
97 			enum opBits = 1;
98 		} else {
99 			enum opBits = 0;
100 		}
101 		enum xShift = opBits + ((X.sizeof * 8) + mask) / 2;
102 		enum mask = (1 << opBits) - 1;
103 		enum bottomSpare = spareBits - opBits;
104 	}
105 	struct xsh {
106 		enum topSpare = 0;
107 		enum bottomSpare = spareBits - topSpare;
108 		enum xShift = (topSpare + X.sizeof * 8) / 2;
109 	}
110 	struct xsl {
111 		enum topSpare = spareBits;
112 		enum bottomSpare = spareBits - topSpare;
113 		enum xShift = (topSpare + X.sizeof * 8) / 2;
114 	}
115 	struct rxs {
116 		enum shift = (I.sizeof - X.sizeof) * 8;
117 		// there must be a simpler way to express this
118 		static if (shift > 64 + 8) {
119 			enum rShiftAmount = I.sizeof - 6;
120 			enum rShiftMask = 63;
121 		} else static if (shift > 32 + 4) {
122 			enum rShiftAmount = I.sizeof - 5;
123 			enum rShiftMask = 31;
124 		} else static if (shift > 16 + 2) {
125 			enum rShiftAmount = I.sizeof - 4;
126 			enum rShiftMask = 15;
127 		} else static if (shift > 8 + 1) {
128 			enum rShiftAmount = I.sizeof - 3;
129 			enum rShiftMask = 7;
130 		} else static if (shift > 4 + 1) {
131 			enum rShiftAmount = I.sizeof - 2;
132 			enum rShiftMask = 3;
133 		} else static if (shift > 2 + 1) {
134 			enum rShiftAmount = I.sizeof - 1;
135 			enum rShiftMask = 1;
136 		} else {
137 			enum rShiftAmount = 0;
138 			enum rShiftMask = 0;
139 		}
140 		enum extraShift = (X.sizeof - shift)/2;
141 	}
142 	struct rxsm {
143 		enum opBits = log2(X.sizeof * 8) - 1;
144 		enum shift = (I.sizeof - X.sizeof) * 8;
145 		enum mask = (1 << opBits) - 1;
146 	}
147 	struct xslrr {
148 		enum opBits = spareBits >= wantedOpBits ? wantedOpBits : spareBits;
149 		enum amplifier = wantedOpBits - opBits;
150 		enum mask = (1 << opBits) - 1;
151 		enum topSpare = spareBits;
152 		enum bottomSpare = spareBits - topSpare;
153 		enum xShift = (topSpare + X.sizeof * 8) / 2;
154 	}
155 }
156 
157 private X xorshift(X, I)(I tmp, uint amt1, uint amt2) {
158 	tmp ^= tmp >> amt1;
159 	return cast(X)(tmp >> amt2);
160 }
161 
162 /// XSH RR (xorshift high, random rotate) - decent performance, slightly better results
163 private X xshrr(X, I)(const I state) {
164 	alias constants = PCGConsts!(X, I).xshrr;
165 	static if (constants.opBits > 0) {
166 		auto rot = (state >> (I.sizeof * 8 - constants.opBits)) & constants.mask;
167 	} else {
168 		enum rot = 0;
169 	}
170 	uint amprot = cast(uint)((rot << constants.amplifier) & constants.mask);
171 	I tmp = state;
172 	return rotr(xorshift!X(tmp, constants.xShift, constants.bottomSpare), amprot);
173 }
174 
175 /// XSH RS (xorshift high, random shift) - decent performance
176 private X xshrs(X, I)(const I state) {
177 	alias constants = PCGConsts!(X, I).xshrs;
178 	static if (constants.opBits > 0) {
179 		uint rshift = (state >> (I.sizeof * 8 - constants.opBits)) & constants.mask;
180 	} else {
181 		uint rshift = 0;
182 	}
183 	I tmp = state;
184 	return xorshift!X(tmp, constants.xShift, cast(uint)(constants.bottomSpare - constants.mask + rshift));
185 }
186 
187 /// XSH (fixed xorshift, high) - don't use this for anything smaller than ulong
188 private X xsh(X, I)(const I state) {
189 	alias constants = PCGConsts!(X, I).xsh;
190 
191 	I tmp = state;
192 	return xorshift!X(tmp, constants.xShift, constants.bottomSpare);
193 }
194 
195 /// XSL (fixed xorshift, low) - don't use this for anything smaller than ulong
196 private X xsl(X, I)(const I state) {
197 	alias constants = PCGConsts!(X, I).xsl;
198 
199 	I tmp = state;
200 	return xorshift!X(tmp, constants.xShift, constants.bottomSpare);
201 }
202 
203 /// RXS (random xorshift)
204 private X rxs(X, I)(const I state) {
205 	alias constants = PCGConsts!(X, I).rxs;
206 	uint rshift = (state >> constants.rShiftAmount) & constants.rShiftMask;
207 	I tmp = state;
208 	return xorshift!X(tmp, cast(uint)(constants.shift + constants.extraShift - rshift), rshift);
209 }
210 
211 /++
212 	RXS M XS (random xorshift, multiply, fixed xorshift)
213 	This has better statistical properties, but supposedly performs worse. This
214 	was not reproducible, however.
215 +/
216 private X rxsmxs(X, I)(const I state) {
217 	X result = rxsm!X(state);
218 	result ^= result >> ((2 * X.sizeof * 8 + 2) / 3);
219 	return result;
220 }
221 
222 /// RXS M (random xorshift, multiply)
223 private X rxsm(X, I)(const I state) {
224 	alias constants = PCGConsts!(X, I).rxsm;
225 	I tmp = state;
226 	static if (constants.opBits > 0) {
227 		uint rshift = (tmp >> (I.sizeof * 8 - constants.opBits)) & constants.mask;
228 	} else {
229 		uint rshift = 0;
230 	}
231 	tmp ^= tmp >> (constants.opBits + rshift);
232 	tmp *= PCGMMultiplier!I;
233 	return cast(X)(tmp >> constants.shift);
234 }
235 
236 /// DXSM (double xorshift, multiply) - newer, better performance for types 2x the size of the largest type the cpu can handle
237 private X dxsm(X, I)(const I state) {
238 	static assert(X.sizeof <= I.sizeof/2, "Output type must be half the size of the state type.");
239 	X hi = cast(X)(state >> ((I.sizeof - X.sizeof) * 8));
240 	X lo = cast(X)state;
241 
242 	lo |= 1;
243 	hi ^= hi >> (X.sizeof * 8 / 2);
244 	hi *= PCGMMultiplier!I;
245 	hi ^= hi >> (3*(X.sizeof * 8 / 4));
246 	hi *= lo;
247 	return hi;
248 }
249 /// XSL RR (fixed xorshift, random rotate) - better performance for types 2x the size of the largest type the cpu can handle
250 private X xslrr(X, I)(const I state) {
251 	alias constants = PCGConsts!(X, I).xslrr;
252 
253 	I tmp = state;
254 	static if (constants.opBits > 0) {
255 		uint rot = (tmp >> (I.sizeof * 8 - constants.opBits)) & constants.mask;
256 	} else {
257 		uint rot = 0;
258 	}
259 	uint amprot = (rot << constants.amplifier) & constants.mask;
260 	return rotr(xorshift!X(tmp, constants.xShift, constants.bottomSpare), amprot);
261 }
262 
263 struct PCG(T, S, alias func, S multiplier = DefaultPCGMultiplier!S, S increment = DefaultPCGIncrement!S) {
264 	private S state;
265 
266 	this(S val) @safe pure nothrow @nogc {
267 		seed(val);
268 	}
269 	void seed(S val) @safe pure nothrow @nogc {
270 		state = cast(S)(val + increment);
271 		popFront();
272 	}
273 	void popFront() @safe pure nothrow @nogc {
274 		state = cast(S)(state * multiplier + increment);
275 	}
276 	T front() const @safe pure nothrow @nogc @property {
277 		return func!T(state);
278 	}
279 	typeof(this) save() @safe pure nothrow @nogc {
280 		return this;
281 	}
282 	enum bool empty = false;
283 	enum bool isUniformRandom = true;
284 	enum T min = T.min;
285 	enum T max = T.max;
286 	const(S) toSiryulType()() const @safe {
287 		return state;
288 	}
289 	static PCG fromSiryulType()(S val) @safe {
290 		PCG result;
291 		result.state = val;
292 		return result;
293 	}
294 }
295 
296 template DefaultPCGMultiplier(T) {
297 	static if (is(T == ubyte)) {
298 		enum DefaultPCGMultiplier = 141;
299 	} else static if (is(T == ushort)) {
300 		enum DefaultPCGMultiplier = 12829;
301 	} else static if (is(T == uint)) {
302 		enum DefaultPCGMultiplier = 747796405;
303 	} else static if (is(T == ulong)) {
304 		enum DefaultPCGMultiplier = 6364136223846793005;
305 	} else static if (is(T == ucent)) {
306 		//enum DefaultPCGMultiplier = 47026247687942121848144207491837523525;
307 	}
308 }
309 
310 template DefaultPCGIncrement(T) {
311 	static if (is(T == ubyte)) {
312 		enum DefaultPCGIncrement = 77;
313 	} else static if (is(T == ushort)) {
314 		enum DefaultPCGIncrement = 47989;
315 	} else static if (is(T == uint)) {
316 		enum DefaultPCGIncrement = 2891336453;
317 	} else static if (is(T == ulong)) {
318 		enum DefaultPCGIncrement = 1442695040888963407;
319 	} else static if (is(T == ucent)) {
320 		//enum DefaultPCGIncrement = 117397592171526113268558934119004209487;
321 	}
322 }
323 
324 private template PCGMMultiplier(T) {
325 	static if (is(T : ubyte)) {
326 		enum PCGMMultiplier = 217;
327 	} else static if (is(T : ushort)) {
328 		enum PCGMMultiplier = 62169;
329 	} else static if (is(T : uint)) {
330 		enum PCGMMultiplier = 277803737;
331 	} else static if (is(T : ulong)) {
332 		enum PCGMMultiplier = 12605985483714917081;
333 	//} else static if (is(T == ucent)) {
334 		//enum PCGMMultiplier = 327738287884841127335028083622016905945;
335 	}
336 }
337 
338 version(arsd_random_unittest) {
339 	private alias AliasSeq(T...) = T;
340 
341 	alias SupportedTypes = AliasSeq!(ubyte, ushort, uint, ulong);
342 	alias SupportedFunctions = AliasSeq!(xshrr, xshrs, xsh, xsl, rxs, rxsmxs, rxsm, xslrr);
343 
344 	static foreach (ResultType; SupportedTypes) {
345 		static foreach (StateType; SupportedTypes) {
346 			static if (StateType.sizeof >= ResultType.sizeof) {
347 				static foreach (Function; SupportedFunctions) {
348 					mixin("alias PCG", int(StateType.sizeof * 8), int(ResultType.sizeof * 8), __traits(identifier, Function), " = PCG!(ResultType, StateType, Function);");
349 				}
350 			}
351 		}
352 	}
353 	alias PCG6432dxsm = PCG!(uint, ulong, dxsm);
354 
355 	@safe unittest {
356 		import std.algorithm : reduce;
357 		import std.datetime.stopwatch : benchmark;
358 		import std.math : pow, sqrt;
359 		import std.random : isSeedable, Mt19937, uniform, uniform01, unpredictableSeed;
360 		import std.stdio : writefln, writeln;
361 		auto seed = unpredictableSeed;
362 
363 		void testRNG(RNG, string name)(uint seed) {
364 			static if (isSeedable!(RNG, uint)) {
365 				auto rng = RNG(seed);
366 			} else static if (isSeedable!(RNG, ushort)) {
367 				auto rng = RNG(cast(ushort)seed);
368 			} else static if (isSeedable!(RNG, ubyte)) {
369 				auto rng = RNG(cast(ubyte)seed);
370 			}
371 			writefln!"--%s--"(name);
372 			double total = 0;
373 			ulong[ubyte] distribution;
374 			void test() {
375 				total += uniform01(rng);
376 				distribution.require(uniform!ubyte(rng), 0)++;
377 			}
378 			auto result = benchmark!(test)(1000000)[0];
379 			writefln!"Benchmark completed in %s"(result);
380 			writeln(total);
381 			double avg = reduce!((a, b) => a + b / distribution.length)(0.0f, distribution);
382 			auto var = reduce!((a, b) => a + pow(b - avg, 2) / distribution.length)(0.0f, distribution);
383 			auto sd = sqrt(var);
384 			writefln!"Average: %s, Standard deviation: %s"(avg, sd);
385 		}
386 
387 		testRNG!(PCG168xshrr, "PCG168xshrr")(seed);
388 		testRNG!(PCG3216xshrr, "PCG3216xshrr")(seed);
389 		testRNG!(PCG6432xslrr, "PCG6432xslrr")(seed);
390 		testRNG!(PCG648rxsmxs, "PCG648rxsmxs")(seed);
391 		testRNG!(PCG6432dxsm, "PCG6432dxsm")(seed);
392 		testRNG!(Mt19937, "Mt19937")(seed);
393 		testRNG!(reasonableDefault, "reasonableDefault")(seed);
394 	}
395 }