#include "rng.h"

#include <limits.h>
#include <math.h>
#include <stdio.h>
#include <sys/time.h>

#ifndef M_PI
#define M_PI 3.14159265258979323846264
#endif

static uint32_t next(rng_t *rng)
{
	uint32_t x = rng->state;
	x ^= x << 13;
	x ^= x >> 17;
	x ^= x << 5;
	return rng->state = x;
}

static double canonical(rng_t *rng)
{
	return (double)next(rng) / (double)UINT32_MAX;
}

static double disc(rng_t *rng)
{
	return 2.0 * canonical(rng) - 1.0;
}

rng_t rng_init(unsigned seed)
{
	struct timeval tv;
	gettimeofday(&tv, nullptr);
	return (rng_t) { .state = tv.tv_usec + seed };
}

vec3_t rng_vec3(rng_t *rng)
{
	const vec3_t v = { disc(rng), disc(rng), disc(rng) };
	return vec3_unit(v);
}

vec3_t rng_gaussian_xy(rng_t *rng, double stddev)
{
	const double r1 = canonical(rng);
	double r2;
	do
		r2 = canonical(rng);
	while (r2 == 0);

	const double theta = 2.0 * M_PI * r1;
	const double mag = stddev * sqrt(-2.0 * log(r2));
	return (vec3_t) {
		.x = mag * cos(theta),
		.y = mag * sin(theta),
	};
}