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 }