/**
 * @file rand.004.c
 * @ingroup experimental
 * Random number generation using custom LCG as bitstream.
 * @date 01/04/2025
 */

#include <stdlib.h>
#include <stdint.h>
#include <limits.h>
#include <assert.h>
#include <math.h>
#include <float.h>
#include <stdio.h>

_Static_assert((-1 & 3) == 3, "Not 2's complement");
_Static_assert(UINT_MAX == (unsigned int) INT_MAX - INT_MIN, "Bad integer");

//
// Utility.
//

#define BITS(x) (sizeof(x) * CHAR_BIT)

int bitlength(unsigned int x)
{
    return x ? BITS(x) - __builtin_clz(x) : 0;
}

//
// Private.
//

#define GENERATOR_BITS 32

static _Thread_local uint_fast64_t generator_state = UINT64_C(0x1234ABCD330E);

static void generator_seed(uint_fast64_t seed)
{
    generator_state = seed;
    generator_state %= (uint_fast64_t) 1 << 48;
}

static uint_fast32_t generator_next(void)
{
    generator_state *= (uint_fast64_t) 25214903917;
    generator_state += (uint_fast64_t) 11;
    generator_state %= (uint_fast64_t) 1 << 48;
    return generator_state >> 16;
}

static uint_fast32_t takebits(uint_fast32_t x, int available, int take)
{
    assert(take >= 1);
    assert(take <= available);
    return (x >> (available - take)) & (UINT32_C(-1) >> (BITS(x) - take));
}

//
// Public.
//

void randseed(uint_fast64_t seed)
{
    generator_seed(seed);
}

uint_fast64_t randbits(int n)
{
    assert(n <= 64);
    static _Thread_local uint_fast32_t reservoir = 0;
    static _Thread_local int available = 0;

    uint_fast64_t r = 0;

    while (n > 0)
    {
        if (available == 0)
        {
            reservoir = generator_next();
            available = GENERATOR_BITS;
        }

        int take = (available < n) ? available : n;
        r = (r << take) | takebits(reservoir, available, take);
        available -= take;
        n -= take;
    }
    return r;
}

unsigned int randuint(unsigned int max)
{
    if (max == 0)
        return 0;

    int n = bitlength(max);

    unsigned int bit_max = -1U >> (BITS(max) - n);
    unsigned int r = randbits(n);

    if (max == bit_max)
        return r;

    unsigned int mod = max + 1;
    unsigned int min = (bit_max + 1 - mod) % mod;

    while (r < min)
        r = randbits(n);
    return r % mod;
}

double randreal(void)
{
    double r = randbits(DBL_MANT_DIG) / exp2(DBL_MANT_DIG);
    assert(r < 1.0);
    return r;
}

_Bool randbool(double p)
{
    assert(0.0 <= p && p <= 1.0);
    return randreal() < p;
}

//
// Test.
//

int uniform_int_distribution_uint(int a, int b)
{
    assert(a <= b);
    return randuint((unsigned int) b - a) + a;
}

int uniform_int_distribution_real(int a, int b)
{
    assert(a <= b);
    return randreal() * (1.0 + b - a) + a;
}

int bernoulli_distribution_50_50(int a, int b)
{
    assert(a == 0 && b == 1);
    return randbool(0.5);
}

unsigned long random_device()
{
    unsigned long r = -1;
    FILE *fp = fopen("/dev/urandom", "r");
    if (fp)
    {
        fread(&r, sizeof r, 1, fp);
        fclose(fp);
    }
    return r;
}

double chi_square(int n, int k, int (*f)(int, int))
{
    int *h = calloc(k, sizeof *h);
    assert(h != 0 || k == 0);
    for (int j = 0; j < n; j++)
    {
        int i = f(0, k-1);
        assert(0 <= i && i < k);
        h[i] += 1;
    }
    double expect = (double) n / k;
    double x2 = 0;
    for (int i = 0; i < k; i++)
        x2 += pow(h[i] - expect, 2) / expect;
    free(h);
    return x2;
}

void test(int m, int n, int df, double x2, int (*f)(int, int),
          const char *header)
{
    puts(header);
    for (int i = 0; i < m; i++)
    {
        double r = chi_square(n, df, f);
        printf(" significant: %-5s [%6.2f][%6.2f]\n",
               r >= x2 ? "true" : "false", r, x2);
    }
}

int main(void)
{
    unsigned long seed = random_device();
    printf("seed: %lu\n", seed);
    randseed(seed);

    int m = 5;
    int n = 10000;

    // PV=0.05 for all tests.

    test(m, n, 100, 124.34, uniform_int_distribution_uint,
         "uniform_int_distribution_uint");
    test(m, n, 100, 124.34, uniform_int_distribution_real,
         "uniform_int_distribution_real");
    test(m, n, 2, 5.99, bernoulli_distribution_50_50,
         "bernoulli_distribution_50_50");
    return 0;
}