### Pseudorandom and -Gaussian

An extension to Pseudorandom number generator

I wanted to see if the random generator I use in my Turtles (by David Bau) has similar properties as the one Reinder uses.

When it seemed to have those properties I thought use them to extend the 'Random wrapper' I use to generate values that form a normal (Gaussian) distributed collection of values instead of a 'random collection'.

I used the Box-Muller transform to get from pseudorandom numbers to normal distributed values from 0 to 1, by code at the first listed URL below. There's also a (faster) method listed that approximates normal distributed values (less precise) using the central limit theorem (using ctrlLimSampleSize).

For completeness I included the random number generator by Randall Munroe (last URL).

I also added a way to skew values from any generator method.

stackoverflow.com/qu…-gaussian-bell-curve
en.wikipedia.org/wiki/box%e2%80%93muller_transform
en.wikipedia.org/wiki/central_limit_theorem
davidbau.com/archive…nd_quintillions.html
xkcd.com/221/

```// Forked from "Pseudorandom number generator" by reinder
// https://turtletoy.net/turtle/a2274fd1fe

// Pseudorandom number generator. Created by Reinder Nijhoff 2024 - @reindernijhoff
//
// https://turtletoy.net/turtle/a2274fd1fe
//

const turtle = new Turtle();

const generator = 2; //min=0 max=4 step=1 (Pseudorandom used by Reinder - see description, Pseudorandom used by Jurgen - see description, Normal distribution using Box-Muller transform - see description, Normal distribution using central limit theorem - see description, Random number generator by Randall Munroe - see description)

let seed = 1005640; // min=1, max=1000000, step=1
const bits = 20; // min=1, max=30, step=1
const samples = 150000; // min=1, max=500000, step=1

const mod = 1<<bits;
const bins = Math.min(1<<8, 1<<bits);
const ksAlpha = 0.05; // 95% confidence level
const acLag = 1; // min=1, max=128, step=1

const ctrlLimSampleSize = 6; //min=2 max=10 step=1

const skew = 0; //min=-.999 max=.999 step=.001

const values = [];

init();
if(generator > 0) { //only when not Reinder's random
R.seed(seed);
}
////////////////////////////////////////////////////////////////
// Pseudorandom number generator. Created by Reinder Nijhoff 2024
// https://turtletoy.net/turtle/a2274fd1fe
////////////////////////////////////////////////////////////////
function random() { // returns a number [0, 1[
let r = 1103515245 * (((seed+=12345) >> 1) ^ seed);
r = 1103515245 * (r ^ (r >> 3));
r = r ^ (r >> 16);
return (r % mod) / mod;
}

//https://www.xkcd.com/221/
function getRandomNumber()
{
return 4; // chosen by fair dice roll.
// guaranteed to be random.
}

const oseed = seed;

function walk(i) {
let r = 0;
switch(generator) {
case 0: r = R.skew(random(), skew); break;
case 1: r = R.skew(R.get(), skew); break;
case 2: r = R.getNormalDistributed(skew); break;
case 3: r = R.skew([ctrlLimSampleSize].map(n => Array.from({length: n}, e => R.get()).reduce((a,c) => a + c, 0)/n).pop(), skew); break;
case 4: r = (getRandomNumber() - 1) / 6; break; //transform dice roll to [0, 1>
}

values.push(r);

drawPoint( i/samples, r, transform(10, -80, 70, 70) );

if (i >= samples) {
const text = new Text(HersheySans1);
turtle.jump(-80,-5);
text.print(turtle, 'Histogram', .35);
turtle.jump(10,-5);
text.print(turtle, 'Noise', .35);

const histogram = createHistogram(values, bins);
drawHistogram(histogram, transform(-80, -80, 70, 70));

const mean = values.reduce((sum, value) => sum + value, 0) / values.length;
const stddev = Math.sqrt(values.map(x => Math.pow(x - mean, 2)).reduce((a, b) => a + b) / values.length);
const corr = autocorrelation(values, acLag);

const stdDevIntervalCounts = Array.from({length: 6}, (e,i) =>
values.filter(e => mean+((i-3)*stddev) <= e && e < mean+((i-2)*stddev)).length/values.length
);

const criticalValue = ksCriticalValue(values, ksAlpha);
const ksStatistic = ksTest(values);

turtle.jump(-80, 10);
text.print(turtle, `Mean: \${mean.toFixed(4)}
StdDev: \${stddev.toFixed(4)}
Values in StdDev intervals from mean (-3 to 3):
\${stdDevIntervalCounts.map(e => (100*e).toFixed(2)).join('% ')}%

Autocorrelation: \${corr.toFixed(4)} (lag \${acLag})
K-S Statistic: \${ksStatistic.toFixed(4)} (\${ksStatistic < criticalValue ? '<':'!! >='} \${criticalValue.toFixed(4)})

Seed: \${oseed}
Bits: \${bits}
Samples: \${samples}`, .52);
}

return i < samples;
}

function transform(l, t, w, h) {
return (x, y) => [l+w*x, t+h*y];
}

function drawPoint(x, y, t) {
const d = 0.0002;
turtle.jump( t(x-d, y-d) );
turtle.goto( t(x+d, y+d) );
}

//
// Noise Statistics, thanks to chatGPT
//
function createHistogram(a, bins) {
const histogram = new Array(bins).fill(0);
for (let value of a) {
histogram[value * bins | 0]++;
}
return histogram;
}

function drawHistogram(histogram, t) {
const n = histogram.filter(v => v > 0).length;
const mean = histogram.reduce((sum, value) => sum + value, 0) / n;

const values = histogram.map(e => e/mean)
const max = values.reduce((a,c) => Math.max(a,c), 0);
const ratio = max < 2? 1: 2/max;

for (let i=0; i<bins; i++) {
turtle.jump( t(i/bins, 1) );
turtle.goto( t(i/bins, 1 - 0.5 * values[i] * ratio));
}
}

function autocorrelation(a, lag) {
const n = a.length;
const mean = a.reduce((sum, value) => sum + value, 0) / n;
let numerator = 0;
let denominator = 0;
for (let i = 0; i < n - lag; i++) {
numerator += (a[i] - mean) * (a[i + lag] - mean);
}
for (let i = 0; i < n; i++) {
denominator += (a[i] - mean) ** 2;
}
return numerator / denominator;
}

function ksCriticalValue(a, alpha) {
return Math.sqrt(-0.5 * Math.log(alpha / 2) / a.length);
}

function ksTest(a) {
a.sort((x, y) => x - y); // Sort the array in ascending order
const n = a.length;
let d = 0;

for (let i = 0; i < n; i++) {
const empiricalCDF = (i + 1) / n; // Empirical CDF
const uniformCDF = a[i]; // Uniform CDF (assuming values in [0, 1])
d = Math.max(d, Math.abs(empiricalCDF - uniformCDF)); // Update K-S statistic
}

return d;
}

////////////////////////////////////////////////////////////////
// Text utility code. Created by Reinder Nijhoff 2024
// https://turtletoy.net/turtle/0f84fd3ae4
////////////////////////////////////////////////////////////////

function Text(fontData) {
const decode = (data) => {
const b = atob(data), a = new Int16Array(b.length / 2), g = {};
for (let i = 0; i < b.length; i+=2) {
a[i / 2 | 0] = b.charCodeAt(i) | (b.charCodeAt(i + 1) << 8);
}
for (let i = 0;i < a.length; i++) {
let u = String.fromCharCode(a[i++]), x = a[i++], d = [], p = [];
for (;a[i] !== -16384; i++) a[i] === 16384 ? (d.push(p), p = []) : p.push(a[i]);
g[u] = { x, d };
}
return g;
}
const rotAdd = (a, b, h) => [Math.cos(h)*a[0] - Math.sin(h)*a[1] + b[0],
Math.cos(h)*a[1] + Math.sin(h)*a[0] + b[1]];

const data = Object.fromEntries(Object.entries(fontData).map(([key, value]) =>
[key, key === 'glyphs' ? decode(value) : value]));

class Text {
print (t, str, scale = 1) {
let pos = t ? [t.x(), t.y()] : [0,0], h = t ? t.h() * Math.PI * 2 / t.fullCircle() : 0,
o = pos, s = scale / data.unitsPerEm;
str.split('').map(c => {
if (c == '\n') {
pos = o = rotAdd([0, 14 * scale], o, h);
return;
}
const glyph = (data.glyphs[c] || data.glyphs[' ']), d = glyph.d;
d.forEach( (p, k) => {
t && t.up();
for (let i=0; i<p.length; i+=2) {
t && t.goto(rotAdd([p[i]*s, p[i+1]*s], pos, h));
t && t.down();
}
});
pos = rotAdd([glyph.x*s, 0], pos, h);
});
}
}

return new Text();
}

function init() {
class Random {
static #apply(seed) {
// Seedable random number generator by David Bau: http://davidbau.com/archives/2010/01/30/random_seeds_coded_hints_and_quintillions.html
!function(a,b,c,d,e,f,g,h,i){function j(a){var b,c=a.length,e=this,f=0,g=e.i=e.j=0,h=e.S=[];for(c||(a=[c++]);d>f;)h[f]=f++;for(f=0;d>f;f++)h[f]=h[g=s&g+a[f%c]+(b=h[f])],h[g]=b;(e.g=function(a){for(var b,c=0,f=e.i,g=e.j,h=e.S;a--;)b=h[f=s&f+1],c=c*d+h[s&(h[f]=h[g=s&g+b])+(h[g]=b)];return e.i=f,e.j=g,c})(d)}function k(a,b){var c,d=[],e=typeof a;if(b&&"object"==e)for(c in a)try{d.push(k(a[c],b-1))}catch(f){}return d.length?d:"string"==e?a:a+"\0"}function l(a,b){for(var c,d=a+"",e=0;e<d.length;)b[s&e]=s&(c^=19*b[s&e])+d.charCodeAt(e++);return n(b)}function m(c){try{return o?n(o.randomBytes(d)):(a.crypto.getRandomValues(c=new Uint8Array(d)),n(c))}catch(e){return[+new Date,a,(c=a.navigator)&&c.plugins,a.screen,n(b)]}}function n(a){return String.fromCharCode.apply(0,a)}var o,p=c.pow(d,e),q=c.pow(2,f),r=2*q,s=d-1,t=c["seed"+i]=function(a,f,g){var h=[];f=1==f?{entropy:!0}:f||{};var o=l(k(f.entropy?[a,n(b)]:null==a?m():a,3),h),s=new j(h);return l(n(s.S),b),(f.pass||g||function(a,b,d){return d?(c[i]=a,b):a})(function(){for(var a=s.g(e),b=p,c=0;q>a;)a=(a+c)*d,b*=d,c=s.g(1);for(;a>=r;)a/=2,b/=2,c>>>=1;return(a+c)/b},o,"global"in f?f.global:this==c)};if(l(c[i](),b),g&&g.exports){g.exports=t;try{o=require("crypto")}catch(u){}}else h&&h.amd&&h(function(){return t})}(this,[],Math,256,6,52,"object"==typeof module&&module,"function"==typeof define&&define,"random");
Math.seedrandom(seed);
}
static seedRandom() { this.#apply(new Date().getMilliseconds()); }
static seedDaily() { this.#apply(new Date().toDateString()); }
static seed(seed) { this.#apply(seed); }
static getInt(min, max) {
if(max == undefined) {
max = min + 1;
min = 0;
}
const [mi, ma] = [Math.min(min, max), Math.max(min, max)];
return (mi + Math.random() * (ma - mi)) | 0;
}
static get(min, max) {
if(min == undefined) {
return Math.random();
}
if(max == undefined) {
max = min;
min = 0;
}
const [mi, ma] = [Math.min(min, max), Math.max(min, max)];
return mi + Math.random() * (ma - mi);
}
static getAngle(l = 1) { return l * this.get(0, 2*Math.PI); }

// Standard Normal variate using Box-Muller transform.
static getGaussian(mean=.5, stdev=.1) {
const u = 1 - this.get(); // Converting [0,1) to (0,1]
const v = this.get();
const z = ( -2.0 * Math.log( u ) )**.5 * Math.cos( 2.0 * Math.PI * v );
// Transform to the desired mean and standard deviation:
return z * stdev + mean;
}
static skew(value, skew = 0) { //skew values (from 0 to 1) by a skew from -1 to 1, respectively right and left skewed (resp more values to left or to right), 0 is not skewed
const v = Math.pow(value, 1-Math.abs(skew));
return skew > 0? v: 1 - v;
}
static getNormalDistributed(skew = 0) { //skew values (from 0 to 1) by a skew from -1 to 1, respectively right and left skewed (resp more values to left or to right), 0 is not skewed
let v = -1;
while(v < 0 || 1 <= v) { v = this.getGaussian(.5, .1) };
return this.skew(v, skew);
}
}
this.R = Random;
}```