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

#include "rng.h"

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

float _rng_next_float(rng *rng) {
    uint64_t state = (uint64_t)rng->seed;
    
    state = (6364136223846793005ULL * state) + 1ULL;
    
    rng->seed = (int64_t)state;

    uint32_t high_bits = (uint32_t)(state >> 40);
    
    return (float)high_bits * (1.0f / 16777216.0f); 
}

rng *get_rng(int64_t seed) {
    rng *r = (rng *)malloc(sizeof(rng));
    if (r == NULL) {
        return NULL;
    }
    r->seed = seed;
    return r;
}

void free_rng(rng* rng) {
    free(rng);
}

float rng_uniform(rng *rng, float a, float b) {
    float r = _rng_next_float(rng);
    
    return a + r * (b - a);
}

float rng_gaussian(rng *rng, float m, float s) {
    static int has_spare = 0;
    static float spare;

    if (has_spare) {
        has_spare = 0;
        return m + s * spare;
    }

    has_spare = 1;
    float u1, u2;

    do {
        u1 = _rng_next_float(rng);
    } while (u1 == 0.0f);
    
    u2 = _rng_next_float(rng);

    float mag = s * sqrt(-2.0f * log(u1));
    float z1 = mag * cos(2.0f * M_PI * u2) + m;
    float z2 = mag * sin(2.0f * M_PI * u2);

    mag = sqrt(-2.0f * log(u1));
    z1 = mag * cos(2.0f * M_PI * u2);
    spare = mag * sin(2.0f * M_PI * u2); 

    return m + s * z1;
}