Skip to content

Instantly share code, notes, and snippets.

@skeeto
Last active November 19, 2022 16:24
Show Gist options
  • Save skeeto/2c4f3935f645ac647ec9d5d48ed49f41 to your computer and use it in GitHub Desktop.
Save skeeto/2c4f3935f645ac647ec9d5d48ed49f41 to your computer and use it in GitHub Desktop.
Approximate Random Gaussian Generator
// 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);
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment