31:T682,md`# Permuted congruential generator A computationally cheap and nicely behaving pseudorandom number generator. It’s a bit tricky in Javascript because there is no native support for 64-bit-integer arithmetic. This is a fork/rewrite of the code in [Job van der Zwan](../@jobleonard)’s [notebook](../@jobleonard/permuted-congruential-generator), shortened and simplified a bit and with some typos fixed. (Also a whole lot fewer code comments, for better or worse.) All credit to Job for porting the C code to Javascript, but any lingering errors below are my fault. > Aside: A stripped down version of Job’s code (with swapped endianness) is twice as fast as mine in Safari. I can’t immediately figure out *why* it would be faster, since it *seems* like that code is doing more arithmetic. It’s substantially *slower* in Chrome/Firefox. Maybe numbers that explicitly look like 16-bit integers get specifically optimized in Safari? > > I added seemingly unnecessary parentheses to the state updating part in both versions, because adding numbers in pairs improves speed by another few percent (better [ILP](https://en.wikipedia.org/wiki/Instruction-level_parallelism) presumably). I also added some seemingly unnecessary \`|0\` to my version, because it seems to improve type inference and also boosts speed by a few percent. For more information on the algorithm, see the [website](https://www.pcg-random.org) Melissa O’Neill set up for the PRNG, the scholarly [paper](https://www.pcg-random.org/paper.html) explaining the details in the context of other PRNGs, or the video recording below of a nice lecture she gave at Stanford about her research.`32:T581,pcgrandom = { // Note that the index order [0, 1, 2, 3] is little-endian const eps = Math.pow(2, -32), m0 = 0x7F2D, m1 = 0x4C95, m2 = 0xF42D, m3 = 0x5851, // 6364136223846793005 a0 = 0x814F, a1 = 0xF767, a2 = 0x7B7E, a3 = 0x1405, // 1442695040888963407 state = Uint16Array.from({length: 4}, () => Math.random() * 0x10000); // random seed return function pcgrandom() { // Advance internal state const s0 = state[0], s1 = state[1], s2 = state[2], s3 = state[3], new0 = a0 + m0*s0 |0, new1 = a1 + m0*s1 + (m1*s0 + (new0 >>> 16)) |0, new2 = a2 + m0*s2 + m1*s1 + (m2*s0 + (new1 >>> 16)) |0, new3 = a3 + m0*s3 + (m1*s2 + m2*s1) + (m3*s0 + (new2 >>> 16)); state[0] = new0, state[1] = new1, state[2] = new2; state[3] = new3; // Calculate output function (XSH RR), uses old state const xorshifted = (s3 << 21) + ((s3 >> 2 ^ s2) << 5) + ((s2 >> 2 ^ s1) >> 11), out_int32 = xorshifted >>> (s3 >> 11) | xorshifted << (-(s3 >> 11) & 31); return eps * (out_int32 >>> 0); } } // In the code above both new2 and new3 can overflow 32 bits, and the // resulting carry bits are lost, but we don’t need to care because // those bits are larger than our 64-bit number anyway. // // new1 can't overflow 32 bits because m0 and m1 are small: // a1 + m0*s1 + m1*s0 + carry // <= (2 + 0x7F2D + 0x4C95) * 0xFFFF == 0xCBC3343C33:T524,pcgrandom_job = { // Note that the index order [0, 1, 2, 3] is little-endian const eps = Math.pow(2, -32), lo = 0xFFFF, m0 = 0x7F2D, m1 = 0x4C95, m2 = 0xF42D, m3 = 0x5851, // 6364136223846793005 a0 = 0x814F, a1 = 0xF767, a2 = 0x7B7E, a3 = 0x1405, // 1442695040888963407 state = Uint16Array.from({length: 4}, () => Math.random() * 0x10000); // random seed return function pcgrandom(){ // Advance internal state const s0 = state[0], s1 = state[1], s2 = state[2], s3 = state[3], new0 = a0 + (s0*m0 & lo), new1 = a1 + (new0 >>> 16) + ((s0*m1 & lo) + (s0*m0 >>> 16)) + (s1*m0 & lo), new2 = a2 + (new1 >>> 16) + ((s0*m2 & lo) + (s0*m1 >>> 16)) + ((s1*m1 & lo) + (s1*m0 >>> 16)) + (s2*m0 & lo), new3 = a3 + (new2 >>> 16) + ((s0*m3 & lo) + (s0*m2 >>> 16)) + ((s1*m2 & lo) + (s1*m1 >>> 16)) + ((s2*m1 & lo) + (s2*m0 >>> 16)) + (s3*m0 & lo); state[0] = new0; state[1] = new1; state[2] = new2; state[3] = new3; // Calculate output function (XSH RR), uses old state const xorshifted = (s3 << 21) + ((s3 >> 2 ^ s2) << 5) + ((s2 >> 2 ^ s1) >> 11), out_int32 = xorshifted >>> (s3 >> 11) | xorshifted << (-(s3 >> 11) & 31); return eps * (out_int32 >>> 0); } }34:[] 2f:{"title":"Permuted congruential generator","cells":[{"id":0,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":140,"value":"html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":220,"value":"md`\nOriginal C code from the [PCG-random website](https://www.pcg-random.org/):\n\n\\`\\`\\`c\n// *Really* minimal PCG32 code / (c) 2014 M.E. O'Neill / pcg-random.org\n// Licensed under Apache License 2.0 (NO WARRANTY, etc. see website)\n\ntypedef struct { uint64_t state; uint64_t inc; } pcg32_random_t;\n\nuint32_t pcg32_random_r(pcg32_random_t* rng)\n{\n uint64_t oldstate = rng->state;\n // Advance internal state\n rng->state = oldstate * 6364136223846793005ULL + (rng->inc|1);\n // Calculate output function (XSH RR), uses old state for max ILP\n uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;\n uint32_t rot = oldstate >> 59u;\n return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));\n}\n\\`\\`\\`\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":90,"value":"$32","pinned":true,"hidden":false,"mode":"ojs"},{"id":258,"value":"{\n const n = 1e7;\n const A = new Array(n).fill(0);\n \n let time = +new Date;\n for (let i = 0; i < n; i++) A[i] = pcgrandom();\n time = (+new Date - time) / 1000;\n \n return md`${(n / 1e6 / time).toFixed(1)} million pseudorandom numbers per second`\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":304,"value":"{\n const n = 1e7;\n const A = new Array(n).fill(0);\n \n let time = +new Date;\n for (let i = 0; i < n; i++) A[i] = pcgrandom_job();\n time = (+new Date - time) / 1000;\n \n return md`${(n / 1e6 / time).toFixed(1)} million pseudorandom numbers per second`\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":300,"value":"$33","pinned":true,"hidden":false,"mode":"ojs"},{"id":6,"value":"Array.from({ length: 1000 }, pcgrandom)","pinned":true,"hidden":false,"mode":"ojs"},{"id":12,"value":"{\n const\n height = 500,\n context = DOM.context2d(width, height);\n\n animate(context.canvas, invalidation, () => {\n context.fillStyle = \"rgba(255,255,255,0.01)\";\n context.fillRect(0, 0, width, height);\n context.fillStyle = \"#000\";\n for (let i = 0; i < width; i++) {\n context.fillRect(width * pcgrandom(), height * pcgrandom(), 1, 1);\n }\n });\n return context.canvas;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":457,"value":"import {animate} from '@jrus/misc'","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W34","stdlib":"1","titling":"auto"}