#include <u.h>
#include <libn.h>

// ----------------------------------------------------------------------------
// Internal structure

uint64 
rol64(uint64 x, int k)
{
	return (x << k) | (x >> (64 - k));
}

typedef struct Rng {
	uint64 s[4];
} Rng;

uint64
xoshiro256ss(Rng *state)
{
	uint64 *s     = state->s;
	uint64 result = rol64(s[1] * 5, 7) * 9;
	uint64 t      = s[1] << 17;

	s[2] ^= s[0];
	s[3] ^= s[1];
	s[1] ^= s[2];
	s[0] ^= s[3];

	s[2] ^= t;
	s[3] = rol64(s[3], 45);

	return result;
}

typedef struct Mix
{
	uint64 s;
} Mix;

uint64 
splitmix64(struct Mix *state) {
	uint64 result = state->s;

	state->s = result + 0x9E3779B97f4A7C15;
	result = (result ^ (result >> 30)) * 0xBF58476D1CE4E5B9;
	result = (result ^ (result >> 27)) * 0x94D049BB133111EB;
	return result ^ (result >> 31);
}

static Rng RNG;

// ----------------------------------------------------------------------------
// Exported functions

/* Initializes the global RNG */

error
rng·init(uint64 seed)
{
	Mix smstate = {seed};

    for (int i=0; i < 4; i++) 
        RNG.s[i] = splitmix64(&smstate);

    return 0;
}

/* Returns a random float64 between 0 and 1 */
double
rng·random(void)
{
    uint64 r = xoshiro256ss(&RNG); 
    return (double)r / (double)UINT64_MAX;
}

double
rng·exponential(double lambda)
{
    double f;

    f = rng·random();
    return -log(1 - f)/lambda;
}

double
rng·normal(void)
{
    double f;
    f = rng·random();

}

/* Returns true or false on success of trial */
bool
rng·bernoulli(double f)
{
    return rng·random() < f;
}

/* Returns a random integer between 0 and max
 * TODO: Modulo throws off uniformity 
 */
uint64
rng·randi(int max)
{
    uint64 r = xoshiro256ss(&RNG); 
    return r % max;
}

uint64
rng·poisson(double mean)
{
    uint64 n;
    double c;

    if(mean<10.0) {
        for(n=0, c=rng·exponential(1.0); c<mean; ++n, c+=rng·exponential(1.0))
            ;
        return n;
    }

    return 0;
}