Tinkering with code and algorithms for fun
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
}
}
Array.from({ length: 1000 }, pcgrandom)
{
const height = 500,
context = DOM.context2d(width, height);
do {
context.fillStyle = "rgba(255,255,255,0.01)";
context.fillRect(0, 0, width, height);
context.fillStyle = "#000";
for (let i = 0; i < width; i++) {
context.fillRect(width * pcgrandom(), height * pcgrandom(), 1, 1);
}
yield context.canvas;
} while (true);
}
Purpose-built for displays of data
Observable is your go-to platform for exploring data and creating expressive data visualizations. Use reactive JavaScript notebooks for prototyping and a collaborative canvas for visual data exploration and dashboard creation.
