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

const turtle = new Turtle();

let seed = 1000; // 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 values = [];


////////////////////////////////////////////////////////////////
// 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;
}

const oseed = seed;

function walk(i) {
    const r = random();
    // const r = (random() + random()) / 2;
    // const r = (random() + random() + random()) / 3;

    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 corr = autocorrelation(values, acLag);
        
        const criticalValue = ksCriticalValue(values, ksAlpha);
        const ksStatistic = ksTest(values);
        
        turtle.jump(-80, 10);
        text.print(turtle, `Mean: ${mean.toFixed(4)}
        
Autocorrelation: ${corr.toFixed(4)} (lag ${acLag})
K-S Statistic: ${ksStatistic.toFixed(4)} (${ksStatistic < criticalValue ? '<':'!! >='} ${criticalValue.toFixed(4)})

Seed: ${oseed}
Bits: ${bits}
Samples: ${samples}`, .6);
    }
    
    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;
    for (let i=0; i<bins; i++) {
        turtle.jump( t(i/bins, 1) );
        turtle.goto( t(i/bins, 1 - 0.5 * histogram[i] / mean));
    }
}

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 rotAdd([0, 10*scale], pos, h);
        }
    }

    return new Text();
}

const HersheySans1 = {"name":"HersheySans1","unitsPerEm":1000,"ascent":800,"descent":-200,"capHeight":500,"xHeight":300,"glyphs":"IADEDgDAIQBODE4MJOZODGj3AEBODIr9GAvF\/k4MAACEDcX+TgyK\/QBAAMAiALATmAgk5pgIxu4AQHASJOZwEsbuAEAAwCMA3Bk6ETjhmAiYCABAnBg44QQQmAgAQJgIPPHcGTzxAEBiB574nBie+ABAAMAkAJwYhA044YQN7AQAQHASOOFwEuwEAECcGNrpJhZk53ASJOaEDSTm2Alk52IH2uliB1DsmAjG7tgJ\/O9ODDzxsBOy8yYW6PRcFyj2nBie+JwYT\/wmFsX+cBIAAIQNAADYCcX+YgdP\/ABAAMAlAIgdiB0k5mIHAAAAQIQNJOYEEKToBBAQ68QOkO1ODMbu2AnG7mIHUOxiB9rpmAhk5xgLJOaEDSTmBBBk57ATpOhcF6ToEhtk54gdJOYAQJwYaPcmFp748BQU+\/AUiv1cFwAA3BkAAEgcxf6IHU\/8iB3U+RIbaPecGGj3AEAAwCYA\/h\/+Hzzx\/h\/878gexu6IHcbuSBz87xIbfPKcGJ74JhZP\/LATxf46EQAATgwAANgJxf6YCIr9YgcU+2IHnviYCCj22Ano9HAS\/O+wE8bu8BRQ7PAU2umwE2TnOhEk5sQOZOeEDdrphA1Q7MQO\/O86EbLzXBdP\/NwZxf5IHAAAyB4AAP4fxf7+H4r9AEAAwCcATgzYCaTomAhk59gJJOYYC2TnGAva6dgJUOyYCJDtAEAAwCgAOhE6ETjhxA64404MZOfYCVDsmAh88pgIaPfYCYr9Tgx2AsQOLAY6EZgIAEAAwCkAOhFiBzjh2Am4404MZOfEDlDsBBB88gQQaPfEDor9Tgx2AtgJLAZiB5gIAEAAwCoAsBOEDZDthA1P\/ABAYgc88bATnvgAQLATPPFiB574AEAAwCsA\/h+wE9rpsBMAAABAmAjo9Mge6PQAQADALADYCdgJFPuYCE\/8YgcU+5gI1PnYCRT72AmK\/WIHAAAAQADALQD+H5gI6PTIHuj0AEAAwC4A2AmYCNT5YgcU+5gIT\/zYCRT7mAjU+QBAAMAvABIbSBw44SwGmAgAQADAMACcGMQOJOYYC2TnmAgQ62IHPPFiB+j0mAgU+xgLxf7EDgAAOhEAAPAUxf5cFxT7nBjo9JwYPPFcFxDr8BRk5zoRJObEDiTmAEAAwDEAnBgYCxDrhA3a6ToRJOY6EQAAAEAAwDIAnBiYCFDsmAgQ69gJpOgYC2TnhA0k5nASJObwFGTnJhak6FwXEOtcF5DtJhb877ATsvNiBwAAnBgAAABAAMAzAJwY2Akk5lwXJOYEEPzvsBP87yYWPPFcF3zynBgo9pwYnvhcF0\/88BTF\/joRAACEDQAA2AnF\/pgIiv1iBxT7AEAAwDQAnBiwEyTmYgdo99wZaPcAQLATJOawEwAAAEAAwDUAnBgmFiTm2Akk5pgIPPHYCfzvhA3G7joRxu7wFPzvXBd88pwYKPacGJ74XBdP\/PAUxf46EQAAhA0AANgJxf6YCIr9YgcU+wBAAMA2AJwYXBfa6SYWZOdwEiTmBBAk5k4MZOfYCRDrmAg88ZgIaPfYCU\/8TgzF\/gQQAAA6EQAA8BTF\/lwXT\/ycGJ74nBho91wXsvPwFDzxOhH87wQQ\/O9ODDzx2Amy85gIaPcAQADANwCcGJwYJOZODAAAAEBiByTmnBgk5gBAAMA4AJwYhA0k5tgJZOeYCNrpmAhQ7NgJxu5ODPzvOhE88fAUfPJcF+j0nBho95wYFPtcF4r9JhbF\/nASAACEDQAA2AnF\/pgIiv1iBxT7Ygdo95gI6PQYC3zyxA488bAT\/O8mFsbuXBdQ7FwX2ukmFmTncBIk5oQNJOYAQADAOQCcGFwXxu4mFnzysBPo9AQQKPbEDij2GAvo9JgIfPJiB8buYgeQ7ZgI2ukYC2TnxA4k5gQQJOawE2TnJhba6VwXxu5cF+j0JhYU+7ATxf4EEAAAhA0AANgJxf6YCE\/8AEAAwDoA2AmYCDzxYgd88pgIsvPYCXzymAg88QBAmAjU+WIHFPuYCE\/82AkU+5gI1PkAQADAOwDYCZgIPPFiB3zymAiy89gJfPKYCDzxAEDYCRT7mAhP\/GIHFPuYCNT52AkU+9gJiv1iBwAAAEAAwDwAiB1IHNrpmAjo9EgcAAAAQADAPQD+H5gIPPHIHjzxAECYCJ74yB6e+ABAAMA+AIgdmAja6Ugc6PSYCAAAAEAAwD8AJhZiB1DsYgcQ65gIpOjYCWTnTgwk5joRJOawE2Tn8BSk6CYWEOsmFpDt8BT877ATPPHEDrLzxA5o9wBAxA6K\/YQNxf7EDgAABBDF\/sQOiv0AQADAQAA0IdwZ\/O+cGJDtJhZQ7HASUOwEEJDtxA7G7oQNfPKEDSj2xA6e+DoR1PnwFNT5XBee+JwYKPYAQHASUOwEEMbuxA588sQOKPYEEJ74OhHU+QBA3BlQ7JwYKPacGJ74EhvU+Ygd1Pn+H2j3NCGy8zQhPPH+H5DtyB4Q60gcpOjcGWTnJhYk5nASJObEDmTnTgyk6NgJEOuYCJDtYgc88WIH6PSYCJ742AkU+04Miv3EDsX+cBIAACYWAADcGcX+SByK\/YgdT\/wAQBIbUOzcGSj23Bme+BIb1PkAQADAQQAmFsQOJObsBAAAAEDEDiTmnBgAAABAmAho9\/AUaPcAQADAQgDcGZgIJOaYCAAAAECYCCTmsBMk5lwXZOecGKTo3BkQ69wZkO2cGPzvXBc88bATfPIAQJgIfPKwE3zyXBey85wY6PTcGWj33BkU+5wYiv1cF8X+sBMAAJgIAAAAQADAQwDcGdwZUOycGNrpJhZk57ATJObEDiTmTgxk59gJ2umYCFDsYgf872IHKPaYCNT52AlP\/E4Mxf7EDgAAsBMAACYWxf6cGE\/83BnU+QBAAMBEANwZmAgk5pgIAAAAQJgIJOY6ESTm8BRk51wX2umcGFDs3Bn879wZKPacGNT5XBdP\/PAUxf46EQAAmAgAAABAAMBFAFwXmAgk5pgIAAAAQJgIJOacGCTmAECYCHzycBJ88gBAmAgAAJwYAAAAQADARgAmFpgIJOaYCAAAAECYCCTmnBgk5gBAmAh88nASfPIAQADARwDcGdwZUOycGNrpJhZk57ATJObEDiTmTgxk59gJ2umYCFDsYgf872IHKPaYCNT52AlP\/E4Mxf7EDgAAsBMAACYWxf6cGE\/83BnU+dwZKPYAQLATKPbcGSj2AEAAwEgAEhuYCCTmmAgAAABA3Bkk5twZAAAAQJgIfPLcGXzyAEAAwEkA2AmYCCTmmAgAAABAAMBKALATcBIk5nAS1Pk6EYr9BBDF\/oQNAAAYCwAAmAjF\/mIHiv0sBtT5LAZo9wBAAMBLANwZmAgk5pgIAAAAQNwZJOaYCGj3AEDEDjzx3BkAAABAAMBMAPAUmAgk5pgIAAAAQJgIAABcFwAAAEAAwE0AiB2YCCTmmAgAAABAmAgk5nASAAAAQEgcJOZwEgAAAEBIHCTmSBwAAABAAMBOABIbmAgk5pgIAAAAQJgIJObcGQAAAEDcGSTm3BkAAABAAMBPABIbxA4k5k4MZOfYCdrpmAhQ7GIH\/O9iByj2mAjU+dgJT\/xODMX+xA4AALATAAAmFsX+nBhP\/NwZ1PkSGyj2Ehv879wZUOycGNrpJhZk57ATJObEDiTmAEAAwFAA3BmYCCTmmAgAAABAmAgk5rATJOZcF2TnnBik6NwZEOvcGcbunBg88VwXfPKwE7LzmAiy8wBAAMBRABIbxA4k5k4MZOfYCdrpmAhQ7GIH\/O9iByj2mAjU+dgJT\/xODMX+xA4AALATAAAmFsX+nBhP\/NwZ1PkSGyj2Ehv879wZUOycGNrpJhZk57ATJObEDiTmAEBwEhT73Bl2AgBAAMBSANwZmAgk5pgIAAAAQJgIJOawEyTmXBdk55wYpOjcGRDr3BmQ7ZwY\/O9cFzzxsBN88pgIfPIAQDoRfPLcGQAAAEAAwFMAnBicGNrpJhZk53ASJOaEDSTm2Alk52IH2uliB1DsmAjG7tgJ\/O9ODDzxsBOy8yYW6PRcFyj2nBie+JwYT\/wmFsX+cBIAAIQNAADYCcX+YgdP\/ABAAMBUALAThA0k5oQNAAAAQOwEJOYmFiTmAEAAwFUAEhuYCCTmmAie+NgJT\/xODMX+BBAAAHASAAAmFsX+nBhP\/NwZnvjcGSTmAEAAwFYAJhbsBCTmxA4AAABAnBgk5sQOAAAAQADAVwCIHSwGJOZODAAAAEBwEiTmTgwAAABAcBIk5pwYAAAAQMgeJOacGAAAAEAAwFgAnBhiByTmnBgAAABAnBgk5mIHAAAAQADAWQAmFuwEJObEDnzyxA4AAABAnBgk5sQOfPIAQADAWgCcGJwYJOZiBwAAAEBiByTmnBgk5gBAYgcAAJwYAAAAQADAWwA6EZgIOOGYCJgIAEDYCTjh2AmYCABAmAg44ToROOEAQJgImAg6EZgIAEAAwFwAOhGxAyTm8BSxAwBAAMBdADoRxA444cQOmAgAQAQQOOEEEJgIAEBiBzjhBBA44QBAYgeYCAQQmAgAQADAXgCwE4QNuOOxA+j0AECEDbjjXBfo9ABAAMBfACYWsQOYCNwZmAgAQADAYADYCdgJUOxiB8buYgc88ZgIfPLYCTzxmAj872IHPPEAQADAYQBcFyYWxu4mFgAAAEAmFnzysBP87zoRxu6EDcbuGAv875gIfPJiByj2Ygee+JgIT\/wYC8X+hA0AADoRAACwE8X+JhZP\/ABAAMBiAFwXmAgk5pgIAAAAQJgIfPIYC\/zvhA3G7joRxu6wE\/zvJhZ88lwXKPZcF574JhZP\/LATxf46EQAAhA0AABgLxf6YCE\/8AEAAwGMAJhYmFnzysBP87zoRxu6EDcbuGAv875gIfPJiByj2Ygee+JgIT\/wYC8X+hA0AADoRAACwE8X+JhZP\/ABAAMBkAFwXJhYk5iYWAAAAQCYWfPKwE\/zvOhHG7oQNxu4YC\/zvmAh88mIHKPZiB574mAhP\/BgLxf6EDQAAOhEAALATxf4mFk\/8AEAAwGUAJhZiByj2JhYo9iYWsvPwFDzxsBP87zoRxu6EDcbuGAv875gIfPJiByj2Ygee+JgIT\/wYC8X+hA0AADoRAACwE8X+JhZP\/ABAAMBmAMQOBBAk5oQNJOYYC2Tn2AkQ69gJAAAAQCwGxu7EDsbuAEAAwGcAXBcmFsbuJhZ2AvAULAawE2IHOhGYCIQNmAgYC2IHAEAmFnzysBP87zoRxu6EDcbuGAv875gIfPJiByj2Ygee+JgIT\/wYC8X+hA0AADoRAACwE8X+JhZP\/ABAAMBoAFwXmAgk5pgIAAAAQJgIsvNODPzvxA7G7nASxu7wFPzvJhay8yYWAAAAQADAaQDYCWIHJOaYCGTn2Akk5pgI7uRiByTmAECYCMbumAgAAABAAMBqAE4M2Akk5hgLZOdODCTmGAvu5NgJJOYAQBgLxu4YC7ED2AliB2IHmAjsBJgIAEAAwGsA8BSYCCTmmAgAAABA8BTG7pgIFPsAQIQNKPYmFgAAAEAAwGwA2AmYCCTmmAgAAABAAMBtAOokmAjG7pgIAAAAQJgIsvNODPzvxA7G7nASxu7wFPzvJhay8yYWAAAAQCYWsvPcGfzvSBzG7v4fxu50IvzvtCOy87QjAAAAQADAbgBcF5gIxu6YCAAAAECYCLLzTgz878QOxu5wEsbu8BT87yYWsvMmFgAAAEAAwG8AXBeEDcbuGAv875gIfPJiByj2Ygee+JgIT\/wYC8X+hA0AADoRAACwE8X+JhZP\/FwXnvhcFyj2JhZ88rAT\/O86EcbuhA3G7gBAAMBwAFwXmAjG7pgImAgAQJgIfPIYC\/zvhA3G7joRxu6wE\/zvJhZ88lwXKPZcF574JhZP\/LATxf46EQAAhA0AABgLxf6YCE\/8AEAAwHEAXBcmFsbuJhaYCABAJhZ88rAT\/O86EcbuhA3G7hgL\/O+YCHzyYgco9mIHnviYCE\/8GAvF\/oQNAAA6EQAAsBPF\/iYWT\/wAQADAcgAEEJgIxu6YCAAAAECYCCj22Al88k4M\/O\/EDsbucBLG7gBAAMBzAPAU8BR88rAT\/O8EEMbuTgzG7pgI\/O9iB3zymAjo9BgLKPY6EWj3sBOe+PAUFPvwFE\/8sBPF\/gQQAABODAAAmAjF\/mIHT\/wAQADAdADEDtgJJObYCRT7GAvF\/oQNAAAEEAAAAEAsBsbuxA7G7gBAAMB1AFwXmAjG7pgIFPvYCcX+TgwAAAQQAABwEsX+JhYU+wBAJhbG7iYWAAAAQADAdgCwEywGxu6EDQAAAEDwFMbuhA0AAABAAMB3ABIbYgfG7k4MAAAAQDoRxu5ODAAAAEA6EcbuJhYAAABAEhvG7iYWAAAAQADAeADwFGIHxu7wFAAAAEDwFMbuYgcAAABAAMB5ALATLAbG7oQNAAAAQPAUxu6EDQAAGAvsBJgIYgcsBpgI7ASYCABAAMB6APAU8BTG7mIHAAAAQGIHxu7wFMbuAEBiBwAA8BQAAABAAMB7ADoRxA444U4MeOIYC7jj2Akk5tgJpOgYCxDrTgxQ7IQNxu6EDTzxGAuy8wBATgx44hgL7uQYC2TnTgza6YQNEOvEDpDtxA7874QNfPKYCOj0hA1o98QO1PnEDk\/8hA3F\/k4MAAAYC3YCGAvsBE4MYgcAQBgLKPaEDZ74hA0U+04Miv0YC8X+2Ak7AdgJsQMYCywGTgxiB8QOmAgAQADAfADYCZgIOOGYCJgIAEAAwH0AOhHYCTjhTgx44oQNuOPEDiTmxA6k6IQNEOtODFDsGAvG7hgLPPGEDbLzAEBODHjihA3u5IQNZOdODNrpGAsQ69gJkO3YCfzvGAt88gQQ6PQYC2j32AnU+dgJT\/wYC8X+TgwAAIQNdgKEDewETgxiBwBAhA0o9hgLnvgYCxT7TgyK\/YQNxf7EDjsBxA6xA4QNLAZODGIH2AmYCABAAMB+AIgdYgee+GIHKPaYCHzyGAs88YQNPPEEEHzy8BQo9lwXaPfcGWj3SBwo9ogdsvMAQGIHKPaYCLLzGAt88oQNfPIEELLz8BRo91wXnvjcGZ74SBxo94gdsvOIHTzxAEAAwA=="};