approximating gaussian white noise with equal-distribution random numbers

jprng.h

/**
 * code from http://burtleburtle.net/bob/rand/smallprng.html
 * "I wrote this PRNG. I place it in the public domain."
 */

#pragma once

#include <stdint.h>

#define JPRNG_RAND_MAX UINT32_MAX

struct jprng_ctx {
	uint32_t a, b, c, d;
};

#define JPRNG_ROT(x,k) (((x)<<(k))|((x)>>(32-(k))))
inline static uint32_t
jprng_rand(struct jprng_ctx *ctx)
{
	uint32_t e;
	
	e      = ctx->a - JPRNG_ROT(ctx->b, 27);
	ctx->a = ctx->b ^ JPRNG_ROT(ctx->c, 17);
	ctx->b = ctx->c + ctx->d;
	ctx->c = ctx->d + e;
	ctx->d = e + ctx->a;

	return ctx->d;
}
#undef JPRNG_ROT

inline static void
jprng_srand(struct jprng_ctx *ctx, uint32_t seed)
{
	uint32_t i;

	ctx->a = 0xf1ea5eed;
	ctx->b = ctx->c = ctx->d = seed;

	for (i = 0; i < 20; ++i)
		jprng_rand(ctx);
}

prngtest.c

/**
 * This is free and unencumbered software released into the public domain.
 *
 * Anyone is free to copy, modify, publish, use, compile, sell, or
 * distribute this software, either in source code form or as a compiled
 * binary, for any purpose, commercial or non-commercial, and by any
 * means.
 *
 * In jurisdictions that recognize copyright laws, the author or authors
 * of this software dedicate any and all copyright interest in the
 * software to the public domain. We make this dedication for the benefit
 * of the public at large and to the detriment of our heirs and
 * successors. We intend this dedication to be an overt act of
 * relinquishment in perpetuity of all present and future rights to this
 * software under copyright law.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
 * IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
 * OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
 * OTHER DEALINGS IN THE SOFTWARE.
 *
 * For more information, please refer to <http://unlicense.org/>
 */

#include <stdlib.h>
#include <stdint.h>
#include <stdio.h>
#include <math.h>

#include "jprng.h"

#define RUNS 1000000
#define RANDOM_SAMPLES 6

#define ARRAY_LENGTH(x) (sizeof(x) / sizeof(*x))

struct jprng_ctx ctx;
int bins[20] = {[0 ... 19] = 0};

static float
randf(struct jprng_ctx *ctx)
{
	uint32_t r = jprng_rand(ctx);
	return (double) r / ((double) JPRNG_RAND_MAX);
}

static float
number_pls(struct jprng_ctx *ctx)
{
	float r;
	int i;

	for (r = 0, i = 0; i < RANDOM_SAMPLES; i++)
		r += randf(ctx);

	r /= (float) i;

	return r;
}

int
main(int argc, char **argv)
{
	float f;
	int i;

	jprng_srand(&ctx, 0xFEEDBEEF);

	for (i = 0; i < RUNS; i++) {
		f = number_pls(&ctx);
		bins[lrintf(floorf(ARRAY_LENGTH(bins) * f))]++;
	}

	f = RUNS / (float) ARRAY_LENGTH(bins);

	for (i = 0; i < ARRAY_LENGTH(bins); i++)
		printf(" %2d: %f\n", i, bins[i] / f);

	return 0;
}