/**
 * tolmacdef.h - header file containing various definitions, including the xoshiro256** PRNG algorithm.
 *
 * Copyright (C) 2025 https://optics-design.com
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * See the COPYING file for the full license text.
 */
#ifndef TOLDEF_H
#define TOLDEF_H

#include <stdint.h>
#include <uchar.h>
#include <time.h>
/* Type definitions for consistent sizing across platforms. Idea taken from https://nullprogram.com/blog/2023/10/08/ (archive link: )  */
typedef uint8_t   u8;
typedef char16_t  c16;
typedef int32_t   b32;
typedef int32_t   i32;
typedef uint32_t  u32;
typedef uint64_t  u64;
typedef float     f32;
typedef double    f64;
typedef uintptr_t uptr;
typedef char      byte;
typedef ptrdiff_t size;
typedef size_t    usize;
/* Utility macros */
#define countof(a)    (size)(sizeof(a) / sizeof(*(a)))
#define lengthof(s)   (countof(s) - 1)
// #define new(a, t, n)  (t *)alloc(a, sizeof(t), _Alignof(t), n) //From Nullprogram, not using currently
#define new(t, n)  (t *)malloc(n*sizeof(t))
#define new_arr(t,arr)  (t)malloc(sizeof(arr))
#endif /* TOLDEF_H */

/*
 * xoshiro256** PRNG algorithm by David Blackman and Sebastiano Vigna
 * Fast, high-quality random number generator with 256-bit state
 * Period: 2^256 - 1
 * References:
 * - https://prng.di.unimi.it/
 * - https://vigna.di.unimi.it/ftp/papers/ScrambledLinear.pdf
 *   Implemented by Claude
 */

// PRNG state structure
typedef struct {
    uint64_t s[4];
} xoshiro256_state;

// Static state for the default global generator
static xoshiro256_state xoshiro_global_state;

// Left rotation operation (helper function)
static inline uint64_t rotl(const uint64_t x, int k) {
    return (x << k) | (x >> (64 - k));
}

// Main xoshiro256** function: generates a random 64-bit value
static inline uint64_t xoshiro256_next(xoshiro256_state *state) {
    const uint64_t result = rotl(state->s[1] * 5, 7) * 9;
    const uint64_t t = state->s[1] << 17;

    state->s[2] ^= state->s[0];
    state->s[3] ^= state->s[1];
    state->s[1] ^= state->s[2];
    state->s[0] ^= state->s[3];

    state->s[2] ^= t;
    state->s[3] = rotl(state->s[3], 45);

    return result;
}

// Jump function - equivalent to 2^128 calls to xoshiro256_next
// Useful for generating non-overlapping sequences for parallel computations
static inline void xoshiro256_jump(xoshiro256_state *state) {
    static const uint64_t JUMP[] = { 0x180ec6d33cfd0aba, 0xd5a61266f0c9392c, 0xa9582618e03fc9aa, 0x39abdc4529b1661c };

    uint64_t s0 = 0;
    uint64_t s1 = 0;
    uint64_t s2 = 0;
    uint64_t s3 = 0;
    for(int i = 0; i < sizeof JUMP / sizeof *JUMP; i++) {
        for(int b = 0; b < 64; b++) {
            if (JUMP[i] & UINT64_C(1) << b) {
                s0 ^= state->s[0];
                s1 ^= state->s[1];
                s2 ^= state->s[2];
                s3 ^= state->s[3];
            }
            xoshiro256_next(state);
        }
    }
    
    state->s[0] = s0;
    state->s[1] = s1;
    state->s[2] = s2;
    state->s[3] = s3;
}

// Splitmix64 algorithm for initializing the state
static inline uint64_t splitmix64(uint64_t *x) {
    uint64_t z = (*x += 0x9e3779b97f4a7c15);
    z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
    z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
    return z ^ (z >> 31);
}

// Initialize xoshiro256** state with a 64-bit seed
static inline void xoshiro256_seed(xoshiro256_state *state, uint64_t seed) {
    uint64_t tmp = seed;
    state->s[0] = splitmix64(&tmp);
    state->s[1] = splitmix64(&tmp);
    state->s[2] = splitmix64(&tmp);
    state->s[3] = splitmix64(&tmp);
}

// Initialize the global state with current time or provided seed
static inline void rng_init(uint64_t seed) {
    if (seed == 0) {
        seed = (uint64_t)time(NULL);
    }
    xoshiro256_seed(&xoshiro_global_state, seed);
}

// Generate a random double between 0.0 and 1.0 (inclusive 0.0, exclusive 1.0)
static inline f64 rng_uniform_f64(void) {
    // Take the top 53 bits (the precision of a double) to get a value in [0,1)
    return (xoshiro256_next(&xoshiro_global_state) >> 11) * (1.0 / 9007199254740992.0);
}

// Generate a random int32 in the range [0, max)
static inline i32 rng_int32(i32 max) {
    // Avoid modulo bias by using rejection sampling
    const uint64_t threshold = (((uint64_t)-1) - ((uint64_t)-1) % max);
    uint64_t r;
    do {
        r = xoshiro256_next(&xoshiro_global_state);
    } while (r >= threshold);
    return (i32)(r % max);
}