From 0a64ff62830495ea21bdafe33b6166c3955d35fd Mon Sep 17 00:00:00 2001 From: admin Date: Fri, 25 Apr 2025 21:13:44 +0200 Subject: First Commit --- include/toldef.h | 147 ++++++++++++++++++++++++++++++++++++++++++++++++++++ include/tolmacdef.h | 147 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 294 insertions(+) create mode 100644 include/toldef.h create mode 100644 include/tolmacdef.h (limited to 'include') diff --git a/include/toldef.h b/include/toldef.h new file mode 100644 index 0000000..a486054 --- /dev/null +++ b/include/toldef.h @@ -0,0 +1,147 @@ +/** + * raymacdef.h - header file containing various definition for rayfile manipulations. + * + * 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 +#include +#include +/* 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); +} + diff --git a/include/tolmacdef.h b/include/tolmacdef.h new file mode 100644 index 0000000..d8a5543 --- /dev/null +++ b/include/tolmacdef.h @@ -0,0 +1,147 @@ +/** + * 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 +#include +#include +/* 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); +} + -- cgit v1.2.3