Approximate Random Gaussian Generator

randn.c

// Approximate Random Gaussian Generator
//   $ cc -O3 -fopenmp -o randn randn.c -lm
//   $ ./randn
// Ref: https://old.reddit.com/r/algorithms/comments/yyz59u
// This is free and unencumbered software released into the public domain.
#include <math.h>
#include <stdint.h>
#include <stdio.h>

struct kahan { double sum, err; };

static struct kahan kahan(struct kahan k, double x)
{
    double y = x - k.err;
    double t = k.sum + y;
    k.err = t - k.sum - y;
    k.sum = t;
    return k;
}

static uint64_t sfc64(uint64_t s[4])
{
    uint64_t r = s[0] + s[1] + s[3]++;
    s[0] = (s[1] >> 11) ^ s[1];
    s[1] = (s[2] <<  3) + s[2];
    s[2] = r + (s[2]<<24 | s[2] >>40);
    return r;
}

static double randn(uint64_t s[4])
{
    uint64_t u = sfc64(s);
    double x = __builtin_popcount(u>>32);
    x += (uint32_t)u * (1 / 4294967296.);
    x -= 16.5;
    x *= 0.3517262290563295;
    return x;
}

int main(void)
{
    #pragma omp parallel for
    for (int i = 0; i < 20; i++) {
        long n = 100000000;
        struct kahan sum1={0}, sum2={0};
        uint64_t s[4] = {i+1, i+2, i+3, i+4};

        sfc64(s); sfc64(s); sfc64(s); sfc64(s);  // mix
        for (long i = 0; i < n; i++) {
            double x = randn(s);
            sum1 = kahan(sum1, x);
            sum2 = kahan(sum2, x*x);
        }

        double mean   = sum1.sum/n;
        double stddev = sqrt(sum2.sum/n - mean*mean);
        printf("%-+23.17g %-+.17g\n", mean, stddev);
    }
}