pcgrandom = { /** Original C implementation by M.E. O'Neill that this is based on: // *Really* minimal PCG32 code / (c) 2014 M.E. O'Neill / pcg-random.org // Licensed under Apache License 2.0 (NO WARRANTY, etc. see website) typedef struct { uint64_t state; uint64_t inc; } pcg32_random_t; uint32_t pcg32_random_r(pcg32_random_t* rng) { uint64_t oldstate = rng->state; // Advance internal state rng->state = oldstate * 6364136223846793005ULL + (rng->inc|1); // Calculate output function (XSH RR), uses old state for max ILP uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u; uint32_t rot = oldstate >> 59u; return (xorshifted >> rot) | (xorshifted << ((-rot) & 31)); } */ // Also, based on the example code on the wiki page // it seems like 1442695040888963407u is a good default value // for increment. // https://en.wikipedia.org/wiki/Permuted_congruential_generator#Example_code // We need Uint64 values, so we resort to using a Uint16Array, // dividing the 64 bits into four chunks of 16 bits. // The reason for that is twofold: first, multiplying a number // that takes up n bits results in a number that takes up // twice as many bits to fit all possible outcomes (so multiplying // two 16 bit numbers results in a number 32 bits wide). // Second, bitshifting and bitmasking in JavaScript is limited to // 32 bits. Because we need all the bits of our multiplications, // and we plan to bitmask and bitshift the results, we need // to use 16 bit input numbers before multiplications. // For the curious, see also the following link for advice on // how to implement big integer addition: // https://www.chosenplaintext.ca/articles/radix-2-51-trick.html // (our use-case is slightly different, and we also implemnt // multiplication and bit-shifting) // So compared to the struct in the above code, the "state" // of pcg32_t is represented in a Uint16Array // (also, the reason I didn't use BigInt was because I thought // that faking shifting and bitmasking with division and modulo // was going to be more of a pain than doing this. Yes really, // and yes I probably was wrong about that). const pcg32_t = Uint16Array.from([ Math.random() * (1<<16), Math.random() * (1<<16), Math.random() * (1<<16), Math.random() * (1<<16) ]); return function(){ // uint64_t oldstate = rng->state; const state0 = pcg32_t[0]; // highest bits const state1 = pcg32_t[1]; const state2 = pcg32_t[2]; const state3 = pcg32_t[3]; // lowest bits // rng->state = oldstate * 6364136223846793005ULL + (rng->inc|1); // 6364136223846793005 in binary is // 101100001010001111101000010110101001100100101010111111100101101; // that splits into: // [ "101100001010001", "1111010000101101", "0100110010010101", "0111111100101101" ] // or in other words: [ 22609, 62509, 19605, 32557 ] const mult0 = 22609; const mult1 = 62509; const mult2 = 19605; const mult3 = 32557; // 1442695040888963407 in binary is // 1010000000101011110110111111011110111011001111000000100000000 // that splits into: // [ "1010000000101", "0111101101111110", "1111011101100111", "1000000100000000" ] // or in other words: [ 5125, 31614, 63335, 33024 ] const inc0 = 5125; const inc1 = 31614; const inc2 = 63335; const inc3 = 33024; // Remember how multiplying by 10 is easy? Just add the // number of zeroes. So that is true for 0b10 (or 2) as // well in binary numbers: 0b100 * 0b100 = 0b10000 // (4 * 4 = 16). Why do I mention this? Because we can // use that to reason about how multiplying our chunks // will behave. For example, state2 and mult2 both represent // 16 bits values with 16 zeroes behind it, and multiplying // them is guaranteed to be a value of 32 bits wide with // 32 zeroes behind it. That means that it has completely // shifted into the upper two chunks. From that we can // conclude that we don't have to multiply all chunks with // all constants, because some will be guaranteed to result // in values that are outside the 64 bits that our container // has. For example multiplying state0 and mult0 is // guaranteed to result in a number with (3 + 3)*16 = 96 bits // of zeroes behind it, which is far outside of our 64 bit // container. So we can skip doing that. // So let's only multiply all the chunks with the constants // that *do* result in values that are of consequence, and // skip the rest. The result is that we only multiply 10 // times instead of the 16 required if we multiplied every // combination of state and multiplication constant. const s0m3 = state0 * mult3; // lower 16 bits go in state0 const s1m2 = state1 * mult2; // lower 16 bits go in state0 const s1m3 = state1 * mult3; // lower 16 bits go in state1, upper in state0 const s2m1 = state2 * mult1; // lower 16 bits go in state0 const s2m2 = state2 * mult2; // lower 16 bits go in state1, upper in state0 const s2m3 = state2 * mult3; // lower 16 bits go in state2, upper in state1 const s3m0 = state3 * mult0; // lower 16 bits go in state0 const s3m1 = state3 * mult1; // lower 16 bits go in state1, upper in state0 const s3m2 = state3 * mult2; // lower 16 bits go in state2, upper in state1 const s3m3 = state3 * mult3; // lower 16 bits go in state3, upper in state2 // Now we add all chunks 16 bits that go into a particular // chunk of state: const new_state3 = (s3m3 & 0xFFFF) + inc3; // And now the final gotcha: the addition used to calculate // new_state3 may have a carry-over into new_state2, so that // has to be added to the chunk above it. Same for the // other chunks after that. const new_state2 = (s3m3 >>> 16) + (s3m2 & 0xFFFF) + (s2m3 & 0xFFFF) + inc2 + (new_state3 >>> 16); const new_state1 = (s3m2 >>> 16) + (s3m1 & 0xFFFF) + (s2m3 >>> 16) + (s2m2 & 0xFFFF) + (s1m3 & 0xFFFF) + inc1 + (new_state2 >>> 16); const new_state0 = (s3m1 >>> 16) + (s3m0 & 0xFFFF) + (s2m2 >>> 16) + (s2m1 & 0xFFFF) + (s1m3 >>> 16) + (s1m2 & 0xFFFF) + (s0m3 & 0xFFFF) + inc0 + (new_state1 >>> 16); // We let the assignment to a Uint16Array handle the truncation pcg32_t[0] = new_state0; pcg32_t[1] = new_state1; pcg32_t[2] = new_state2; pcg32_t[3] = new_state3; // If we inspect this line again: // uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u; // we can see that we can completely ignore state3, because the // bottom 27 bits are completely shifted away, and only have to // shift the upper three chunks. Similarly, >> 18 completely // shifts away state0 before masking it with itself, and so on. // What we'll end up with is: const xorshifted = ((state0 << 21) >>> 0) + // see below for >>> 0 explanation (((state0 >> 2) ^ state1) << 5) + (((state1 >> 2) ^ state2) >> 11); // uint32_t rot = oldstate >> 59u; // The lower three chunks are all shifted away, so that's 48 bits // we don't need to shift. So oldstate >> 59 becomes state0 >> (59 - 48), // which is: const rot = state0 >> 11; // in javascript, 0xFFFFFFFF | 0 === -1, or in other words: OR-masking the // most significant bit results in a signed 32 bit number. To undo the effect // of creating a negative number via OR-masking we logical shift to the // right by zero. return (((xorshifted >>> rot) | (xorshifted << ((-rot) & 31))) >>> 0) / 0x100000000; // TODO: Is there a more elegant approach than dividing by 0x100000000? // Using a Float64Array and a Uint8Array with the same backing buffer, // we can directly access the bits of a double. Sadly, just dumping // all 32 bits in the 52 fraction bits of a double won't work, since floats // don't work that way. // See also: // https://mumble.net/~campbell/2014/04/28/uniform-random-float } }