59 lines
996 B
C
59 lines
996 B
C
#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 disc(rng_t *rng)
|
|
{
|
|
return 2.0 * rng_canon(rng) - 1.0;
|
|
}
|
|
|
|
rng_t rng_init(unsigned seed)
|
|
{
|
|
struct timeval tv;
|
|
gettimeofday(&tv, nullptr);
|
|
return (rng_t) { .state = tv.tv_usec + seed };
|
|
}
|
|
|
|
double rng_canon(rng_t *rng)
|
|
{
|
|
return (double)next(rng) / (double)UINT32_MAX;
|
|
}
|
|
|
|
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 = rng_canon(rng);
|
|
double r2;
|
|
do
|
|
r2 = rng_canon(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),
|
|
};
|
|
}
|