1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
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 <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);
}
|
