Add guf_rand with 32-bit operations (using xoshiro128** 1.1)
This commit is contained in:
jun
parent
f911933588
commit
22ebd8aa05
@ -29,8 +29,8 @@ typedef enum guf_cpy_opt {
|
||||
#define GUF_SWAP(TYPE, val_a, val_b) do {TYPE guf_swap_tmp = val_a; val_a = val_b; val_b = guf_swap_tmp;} while (0);
|
||||
#define GUF_STATIC_BUF_SIZE(BUF) (sizeof((BUF)) / (sizeof((BUF)[0])))
|
||||
|
||||
#define GUF_MIN(X, Y) ((X) <= (Y) ? (X) : (Y))
|
||||
#define GUF_MAX(X, Y) ((X) >= (Y) ? (X) : (Y))
|
||||
#define GUF_MIN(X, Y) ((X) < (Y) ? (X) : (Y))
|
||||
#define GUF_MAX(X, Y) ((X) > (Y) ? (X) : (Y))
|
||||
#define GUF_CLAMP(X, MIN, MAX) GUF_MAX(GUF_MIN((X), (MAX)), (MIN))
|
||||
|
||||
// The GUF_CAT/GUF_TOK_CAT indirection is necessary because the ## operation alone does not evaluate the macro arguments.
|
||||
|
||||
@ -1,9 +1,6 @@
|
||||
#ifndef GUF_HASH_H
|
||||
#define GUF_HASH_H
|
||||
#include "guf_common.h"
|
||||
#include "guf_assert.h"
|
||||
|
||||
// #define GUF_USE_32_BIT_HASH /* Define GUF_USE_32_BIT_HASH to make guflib use 32 bit hashes */
|
||||
|
||||
/*
|
||||
FNV-1a (32-bit and 64-bit) hash functions.
|
||||
@ -11,50 +8,67 @@
|
||||
cf. http://www.isthe.com/chongo/tech/comp/fnv/ (last retrieved: 2023-11-30)
|
||||
*/
|
||||
|
||||
#define GUF_HASH32_INIT 2166136261ul
|
||||
#define GUF_HASH64_INIT 14695981039346656037ull
|
||||
#ifdef GUF_IMPL_STATIC
|
||||
#define GUF_FN_KEYWORDS static
|
||||
#else
|
||||
#define GUF_FN_KEYWORDS
|
||||
#endif
|
||||
|
||||
static inline uint32_t guf_hash32(const void *data, ptrdiff_t num_bytes, uint32_t hash)
|
||||
#define GUF_HASH32_INIT UINT32_C(2166136261)
|
||||
#define GUF_HASH64_INIT UINT64_C(14695981039346656037)
|
||||
|
||||
GUF_FN_KEYWORDS uint32_t guf_hash32(const void *data, ptrdiff_t num_bytes, uint32_t hash); // FNV-1a (32 bit)
|
||||
GUF_FN_KEYWORDS uint64_t guf_hash64(const void *data, ptrdiff_t num_bytes, uint64_t hash); // FNV-1a (64 bit)
|
||||
|
||||
#ifdef GUF_HASH_32_BIT
|
||||
typedef uint32_t guf_hash_size_t;
|
||||
#define GUF_HASH_INIT GUF_HASH32_INIT
|
||||
#define GUF_HASH_MAX UINT32_MAX
|
||||
static inline guf_hash_size_t guf_hash(const void *data, ptrdiff_t num_bytes, uint32_t hash) {
|
||||
return guf_hash32(data, num_bytes, hash);
|
||||
}
|
||||
#else
|
||||
typedef uint64_t guf_hash_size_t;
|
||||
#define GUF_HASH_INIT GUF_HASH64_INIT
|
||||
#define GUF_HASH_MAX UINT64_MAX
|
||||
static inline guf_hash_size_t guf_hash(const void *data, ptrdiff_t num_bytes, uint64_t hash) {
|
||||
return guf_hash64(data, num_bytes, hash);
|
||||
}
|
||||
#endif
|
||||
|
||||
#endif
|
||||
|
||||
#if defined(GUF_IMPL) || defined(GUF_IMPL_STATIC)
|
||||
#include "guf_assert.h"
|
||||
|
||||
GUF_FN_KEYWORDS uint32_t guf_hash32(const void *data, ptrdiff_t num_bytes, uint32_t hash)
|
||||
{
|
||||
GUF_ASSERT_RELEASE(data);
|
||||
GUF_ASSERT_RELEASE(num_bytes > 0);
|
||||
const unsigned char *data_bytes = (const unsigned char*)data; // This does not break strict-aliasing rules I think...
|
||||
const uint32_t FNV_32_PRIME = 16777619ul;
|
||||
for (size_t i = 0; i < (size_t)num_bytes; ++i) {
|
||||
for (ptrdiff_t i = 0; i < num_bytes; ++i) {
|
||||
hash ^= data_bytes[i];
|
||||
hash *= FNV_32_PRIME;
|
||||
}
|
||||
return hash;
|
||||
}
|
||||
|
||||
static inline uint64_t guf_hash64(const void *data, ptrdiff_t num_bytes, uint64_t hash)
|
||||
GUF_FN_KEYWORDS uint64_t guf_hash64(const void *data, ptrdiff_t num_bytes, uint64_t hash)
|
||||
{
|
||||
GUF_ASSERT_RELEASE(data);
|
||||
GUF_ASSERT_RELEASE(num_bytes > 0);
|
||||
const unsigned char *data_bytes = (const unsigned char*)data; // This does not break strict-aliasing rules I think...
|
||||
const uint64_t FNV_64_PRIME = 1099511628211ull;
|
||||
for (size_t i = 0; i < (size_t)num_bytes; ++i) {
|
||||
for (ptrdiff_t i = 0; i < num_bytes; ++i) {
|
||||
hash ^= data_bytes[i];
|
||||
hash *= FNV_64_PRIME;
|
||||
}
|
||||
return hash;
|
||||
}
|
||||
|
||||
#ifdef GUF_HASH_USE_32_BIT
|
||||
typedef uint32_t guf_hash_size_t;
|
||||
#define GUF_HASH_INIT GUF_HASH32_INIT
|
||||
#undef GUF_IMPL
|
||||
#undef GUF_IMPL_STATIC
|
||||
#endif /* endif GUF_IMPL/GUF_IMPL_STATIC */
|
||||
|
||||
static inline guf_hash_size_t guf_hash(const void *data, ptrdiff_t num_bytes, uint32_t hash) {
|
||||
return guf_hash32(data, num_bytes, hash);
|
||||
}
|
||||
#define GUF_HASH_MAX UINT32_MAX
|
||||
#else
|
||||
typedef uint64_t guf_hash_size_t;
|
||||
#define GUF_HASH_INIT GUF_HASH64_INIT
|
||||
static inline guf_hash_size_t guf_hash(const void *data, ptrdiff_t num_bytes, uint64_t hash) {
|
||||
return guf_hash64(data, num_bytes, hash);
|
||||
}
|
||||
#define GUF_HASH_MAX UINT64_MAX
|
||||
#endif
|
||||
|
||||
#endif
|
||||
#undef GUF_FN_KEYWORDS
|
||||
|
||||
@ -1,45 +0,0 @@
|
||||
#ifndef GUF_INT_H
|
||||
#define GUF_INT_H
|
||||
#include "guf_common.h"
|
||||
#include "guf_assert.h"
|
||||
|
||||
// #define GUF_DECLARE_MIN_MAX_CLAMP(int_type, int_type_name)\
|
||||
// static inline int_type GUF_CAT(guf_min_, int_type_name)(int_type a, int_type b);\
|
||||
// static inline int_type GUF_CAT(guf_max_, int_type_name)(int_type a, int_type b);\
|
||||
// static inline int_type GUF_CAT(guf_clamp_, int_type_name)(int_type x, int_type min, int_type max);
|
||||
|
||||
#define GUF_DEFINE_MIN_MAX_CLAMP(int_type, int_type_name)\
|
||||
static inline int_type GUF_CAT(guf_min_, int_type_name)(int_type a, int_type b) {return a >= b ? a : b;}\
|
||||
static inline int_type GUF_CAT(guf_max_, int_type_name)(int_type a, int_type b) {return a >= b ? a : b;}\
|
||||
static inline int_type GUF_CAT(guf_clamp_, int_type_name)(int_type x, int_type min, int_type max) {if (x < min) {return min;} if (x > max) {return max;} return x;}
|
||||
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(char, char)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(int, int)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(int8_t, i8)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(int16_t, i16)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(int32_t, i32)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(int64_t, i64)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(ptrdiff_t, ptrdiff_t)
|
||||
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(unsigned char, uchar)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(unsigned, unsigned)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(uint8_t, u8)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(uint16_t, u16)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(uint32_t, u32)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(uint64_t, u64)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(size_t, size_t)
|
||||
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(float, f32)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(double, f64)
|
||||
|
||||
// static inline int guf_abs_int(int x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > INT_MIN); return -x;} // I would not drink that...
|
||||
// static inline long guf_abs_long(long x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > LONG_MIN); return -x;}
|
||||
// static inline long long guf_abs_long_long(long long x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > LLONG_MIN); return -x;}
|
||||
// static inline int8_t guf_abs_i8 (int8_t x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > INT8_MIN); return -x;}
|
||||
// static inline int16_t guf_abs_i16(int16_t x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > INT16_MIN); return -x;}
|
||||
// static inline int32_t guf_abs_i32(int32_t x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > INT32_MIN); return -x;}
|
||||
// static inline int64_t guf_abs_i64(int64_t x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > INT64_MIN); return -x;}
|
||||
// static inline ptrdiff_t guf_abs_ptrdiff(ptrdiff_t x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > PTRDIFF_MIN); return -x;}
|
||||
|
||||
#undef GUF_DEFINE_MIN_MAX_CLAMP
|
||||
#endif
|
||||
125
src/guf_math.h
125
src/guf_math.h
@ -1,71 +1,104 @@
|
||||
#ifndef GUF_MATH_H
|
||||
#define GUF_MATH_H
|
||||
#include "guf_common.h"
|
||||
#include "guf_assert.h"
|
||||
|
||||
#define GUF_PI 3.14159265358979323846264338327950288
|
||||
|
||||
// Rotate left.
|
||||
static inline uint64_t guf_rotl_u64(uint64_t x, int k)
|
||||
{
|
||||
return (x << k) | (x >> (64 - k));
|
||||
}
|
||||
static inline uint64_t guf_rotl_u64(uint64_t x, int k) {return (x << k) | (x >> (64 - k));}
|
||||
static inline uint32_t guf_rotl_u32(uint32_t x, int k) {return (x << k) | (x >> (32 - k));}
|
||||
|
||||
static inline uint32_t guf_rotl_u32(uint32_t x, int k) {
|
||||
return (x << k) | (x >> (32 - k));
|
||||
}
|
||||
#define GUF_DEFINE_MIN_MAX_CLAMP(int_type, int_type_name)\
|
||||
static inline int_type GUF_CAT(guf_min_, int_type_name)(int_type a, int_type b) {return a < b ? a : b;}\
|
||||
static inline int_type GUF_CAT(guf_max_, int_type_name)(int_type a, int_type b) {return a > b ? a : b;}\
|
||||
static inline int_type GUF_CAT(guf_clamp_, int_type_name)(int_type x, int_type min, int_type max) {if (x < min) {return min;} if (x > max) {return max;} return x;}
|
||||
|
||||
static inline float guf_clamp_f32(float x, float min, float max)
|
||||
{
|
||||
if (x < min) return min;
|
||||
if (x > max) return max;
|
||||
return x;
|
||||
}
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(char, char)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(int, int)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(long, long)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(long long, long_long)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(int8_t, i8)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(int16_t, i16)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(int32_t, i32)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(int64_t, i64)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(ptrdiff_t, ptrdiff_t)
|
||||
|
||||
static inline double guf_clamp_f64(double x, double min, double max)
|
||||
{
|
||||
if (x < min) return min;
|
||||
if (x > max) return max;
|
||||
return x;
|
||||
}
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(unsigned char, uchar)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(unsigned, unsigned)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(unsigned long, ulong)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(unsigned long long, ulong_long)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(uint8_t, u8)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(uint16_t, u16)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(uint32_t, u32)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(uint64_t, u64)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(size_t, size_t)
|
||||
|
||||
static inline float guf_lerp_f32(float a, float b, float alpha)
|
||||
{
|
||||
return (1 - alpha) * a + alpha * b;
|
||||
}
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(float, f32)
|
||||
GUF_DEFINE_MIN_MAX_CLAMP(double, f64)
|
||||
|
||||
static inline double guf_lerp_f64(double a, double b, double alpha)
|
||||
{
|
||||
return (1 - alpha) * a + alpha * b;
|
||||
}
|
||||
|
||||
static inline float guf_smoothstep_f32(float edge0, float edge1, float x)
|
||||
{
|
||||
GUF_ASSERT(edge0 != edge1);
|
||||
x = guf_clamp_f32((x - edge0) / (edge1 - edge0), 0, 1);
|
||||
return x * x * (3.f - 2.f * x);
|
||||
}
|
||||
|
||||
static inline float guf_smootherstep_f32(float edge0, float edge1, float x)
|
||||
{
|
||||
GUF_ASSERT(edge0 != edge1);
|
||||
x = guf_clamp_f32((x - edge0) / (edge1 - edge0), 0, 1);
|
||||
return x * x * x * (x * (6.f * x - 15.f) + 10.f);
|
||||
}
|
||||
#undef GUF_DEFINE_MIN_MAX_CLAMP
|
||||
|
||||
static inline int guf_abs_int(int x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > INT_MIN); return -x;} // I would not drink that...
|
||||
static inline long guf_abs_long(long x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > LONG_MIN); return -x;}
|
||||
static inline long long guf_abs_long_long(long long x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > LLONG_MIN); return -x;}
|
||||
static inline int8_t guf_abs_i8 (int8_t x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > INT8_MIN); return -x;}
|
||||
static inline int16_t guf_abs_i16(int16_t x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > INT16_MIN); return -x;}
|
||||
static inline int32_t guf_abs_i32(int32_t x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > INT32_MIN); return -x;}
|
||||
static inline int64_t guf_abs_i64(int64_t x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > INT64_MIN); return -x;}
|
||||
static inline ptrdiff_t guf_abs_ptrdiff(ptrdiff_t x) {if (x >= 0) {return x;} GUF_ASSERT_RELEASE(x > PTRDIFF_MIN); return -x;}
|
||||
|
||||
static inline uint32_t guf_uabs_i32(int32_t x) {if (x >= 0) {return x;} else if (x == INT32_MIN) {return (uint32_t)INT32_MAX + 1;} else {return -x;}}
|
||||
static inline uint64_t guf_uabs_i64(int64_t x) {if (x >= 0) {return x;} else if (x == INT64_MIN) {return (uint64_t)INT64_MAX + 1;} else {return -x;}}
|
||||
static inline unsigned char guf_uabs_char(char x) {if (x >= 0) {return x;} else if (x == CHAR_MIN) {return (unsigned char)CHAR_MAX + 1;} else {return -x;}}
|
||||
static inline unsigned guf_uabs_int(int x) {if (x >= 0) {return x;} else if (x == INT_MIN) {return (unsigned)INT_MAX + 1;} else {return -x;}}
|
||||
static inline uint8_t guf_uabs_i8(int8_t x) {if (x >= 0) {return x;} else if (x == INT8_MIN) {return (uint8_t)INT8_MAX + 1;} else {return -x;}}
|
||||
static inline uint16_t guf_uabs_i16(int16_t x) {if (x >= 0) {return x;} else if (x == INT16_MIN) {return (uint16_t)INT16_MAX + 1;} else {return -x;}}
|
||||
static inline uint32_t guf_uabs_i32(int32_t x) {if (x >= 0) {return x;} else if (x == INT32_MIN) {return (uint32_t)INT32_MAX + 1;} else {return -x;}}
|
||||
static inline uint64_t guf_uabs_i64(int64_t x) {if (x >= 0) {return x;} else if (x == INT64_MIN) {return (uint64_t)INT64_MAX + 1;} else {return -x;}}
|
||||
|
||||
static inline unsigned char guf_absdiff_char(char a, char b) {return a > b ? (unsigned char)a - (unsigned char)b : (unsigned char)b - (unsigned char)a;}
|
||||
static inline unsigned guf_absdiff_int(int a, int b) {return a > b ? (unsigned)a - (unsigned)b : (unsigned)b - (unsigned)a;}
|
||||
static inline uint8_t guf_absdiff_i8(int8_t a, int8_t b) {return a > b ? (uint8_t)a - (uint8_t)b : (uint8_t)b - (uint8_t)a;}
|
||||
static inline uint16_t guf_absdiff_i16(int16_t a, int16_t b) {return a > b ? (uint16_t)a - (uint16_t)b : (uint16_t)b - (uint16_t)a;}
|
||||
static inline uint32_t guf_absdiff_i32(int32_t a, int32_t b) {return a > b ? (uint32_t)a - (uint32_t)b : (uint32_t)b - (uint32_t)a;}
|
||||
static inline uint64_t guf_absdiff_i64(int64_t a, int64_t b) {return a > b ? (uint64_t)a - (uint64_t)b : (uint64_t)b - (uint64_t)a;}
|
||||
|
||||
// An alternative lerp would be a + alpha * (b - a) (advantage: would be weakly monotonic, disadvantage: would not guarantee a for alpha = 0 and b for alpha = 1)
|
||||
static inline float guf_lerp_f32(float a, float b, float alpha) {return (1 - alpha) * a + alpha * b;}
|
||||
static inline double guf_lerp_f64(double a, double b, double alpha) {return (1 - alpha) * a + alpha * b;}
|
||||
|
||||
/*
|
||||
smoothstep interpolation, cf. https://en.wikipedia.org/wiki/Smoothstep (last-retrieved 2025-02-18)
|
||||
*/
|
||||
|
||||
static inline float guf_smoothstep_f32(float edge0, float edge1, float x)
|
||||
{
|
||||
if (edge0 == edge1) { // Prevent division by zero.
|
||||
return 1;
|
||||
}
|
||||
x = guf_clamp_f32((x - edge0) / (edge1 - edge0), 0, 1); // Bring in range [0, 1]
|
||||
return x * x * (3.f - 2.f * x);
|
||||
}
|
||||
static inline float guf_smootherstep_f32(float edge0, float edge1, float x)
|
||||
{
|
||||
if (edge0 == edge1) { // Prevent division by zero.
|
||||
return 1;
|
||||
}
|
||||
x = guf_clamp_f32((x - edge0) / (edge1 - edge0), 0, 1); // Bring in range [0, 1]
|
||||
return x * x * x * (x * (6.f * x - 15.f) + 10.f);
|
||||
}
|
||||
|
||||
static inline double guf_smoothstep_f64(double edge0, double edge1, double x)
|
||||
{
|
||||
if (edge0 == edge1) { // Prevent division by zero.
|
||||
return 1;
|
||||
}
|
||||
x = guf_clamp_f64((x - edge0) / (edge1 - edge0), 0, 1); // Bring in range [0, 1]
|
||||
return x * x * (3.0 - 2.0 * x);
|
||||
}
|
||||
static inline double guf_smootherstep_f64(double edge0, double edge1, double x)
|
||||
{
|
||||
if (edge0 == edge1) { // Prevent division by zero.
|
||||
return 1;
|
||||
}
|
||||
x = guf_clamp_f64((x - edge0) / (edge1 - edge0), 0, 1); // Bring in range [0, 1]
|
||||
return x * x * x * (x * (6.0 * x - 15.0) + 10.0);
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
194
src/guf_rand.h
194
src/guf_rand.h
@ -1,16 +1,5 @@
|
||||
#ifndef GUF_RAND_H
|
||||
#define GUF_RAND_H
|
||||
#include "guf_common.h"
|
||||
#include "guf_assert.h"
|
||||
#include "guf_math.h"
|
||||
#include <math.h>
|
||||
#include <float.h>
|
||||
|
||||
#define GUF_RAND_MAX UINT64_MAX
|
||||
|
||||
typedef struct guf_randstate { // State for xoshiro256** 1.0
|
||||
uint64_t s[4];
|
||||
} guf_randstate;
|
||||
|
||||
#ifdef GUF_IMPL_STATIC
|
||||
#define GUF_FN_KEYWORDS static
|
||||
@ -18,25 +7,40 @@ typedef struct guf_randstate { // State for xoshiro256** 1.0
|
||||
#define GUF_FN_KEYWORDS
|
||||
#endif
|
||||
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
#define GUF_RAND_MAX UINT32_MAX
|
||||
typedef struct guf_randstate { // State for xoshiro128** 1.1
|
||||
uint32_t s[4];
|
||||
} guf_randstate;
|
||||
#else
|
||||
#define GUF_RAND_MAX UINT64_MAX
|
||||
typedef struct guf_randstate { // State for xoshiro256** 1.0
|
||||
uint64_t s[4];
|
||||
} guf_randstate;
|
||||
#endif
|
||||
|
||||
GUF_FN_KEYWORDS uint64_t guf_rand_splitmix64(uint64_t *state);
|
||||
|
||||
GUF_FN_KEYWORDS void guf_randstate_init(guf_randstate *state, uint64_t seed);
|
||||
void guf_randstate_jump(guf_randstate *state); // Advance the state; equivalent to 2^128 calls to guf_rand_u64(state)
|
||||
|
||||
// uniform distributions using xoshiro256** 1.0
|
||||
GUF_FN_KEYWORDS uint64_t guf_rand_u64(guf_randstate *state); // [0, GUF_RAND_MAX]
|
||||
// uniform distributions
|
||||
GUF_FN_KEYWORDS uint32_t guf_rand_u32(guf_randstate *state); // [0, UINT32_MAX]
|
||||
GUF_FN_KEYWORDS uint64_t guf_rand_u64(guf_randstate *state); // [0, UINT64_MAX]
|
||||
GUF_FN_KEYWORDS double guf_rand_f64(guf_randstate *state); // [0.0, 1.0)
|
||||
GUF_FN_KEYWORDS float guf_rand_f32(guf_randstate *state); // [0.f, 1.f)
|
||||
|
||||
// return true with a probability of p, false with a probability of (1 - p)
|
||||
GUF_FN_KEYWORDS bool guf_rand_bernoulli_trial(guf_randstate *state, double p);
|
||||
GUF_FN_KEYWORDS bool guf_rand_bernoulli_trial_f32(guf_randstate *state, float p);
|
||||
GUF_FN_KEYWORDS bool guf_rand_bernoulli_trial_f64(guf_randstate *state, double p);
|
||||
GUF_FN_KEYWORDS bool guf_rand_flip(guf_randstate *state); // Fair coin flip (bernoulli trial with p == 0.5)
|
||||
|
||||
GUF_FN_KEYWORDS double guf_randrange_f64(guf_randstate *state, double min, double end); // [min, end)
|
||||
GUF_FN_KEYWORDS float guf_randrange_f32(guf_randstate *state, float min, float end); // [min, end)
|
||||
GUF_FN_KEYWORDS int32_t guf_randrange_i32(guf_randstate *state, int32_t min, int32_t max); // [min, max]
|
||||
GUF_FN_KEYWORDS uint32_t guf_randrange_u32(guf_randstate *state, uint32_t min, uint32_t max); //[min, max]
|
||||
GUF_FN_KEYWORDS int64_t guf_randrange_i64(guf_randstate *state, int64_t min, int64_t max); // [min, max]
|
||||
GUF_FN_KEYWORDS double guf_randrange_f64(guf_randstate *state, double min, double end); // [min, end)
|
||||
GUF_FN_KEYWORDS float guf_randrange_f32(guf_randstate *state, float min, float end); // [min, end)
|
||||
|
||||
GUF_FN_KEYWORDS int32_t guf_randrange_i32(guf_randstate *state, int32_t min, int32_t max); // [min, max]
|
||||
GUF_FN_KEYWORDS uint32_t guf_randrange_u32(guf_randstate *state, uint32_t min, uint32_t max); // [min, max]
|
||||
GUF_FN_KEYWORDS int64_t guf_randrange_i64(guf_randstate *state, int64_t min, int64_t max); // [min, max]
|
||||
|
||||
// normal distributions
|
||||
GUF_FN_KEYWORDS void guf_rand_normal_sample_f64(guf_randstate *state, double mean, double std_dev, double *result, ptrdiff_t n);
|
||||
@ -44,7 +48,14 @@ GUF_FN_KEYWORDS void guf_rand_normal_sample_f32(guf_randstate *state, float mean
|
||||
GUF_FN_KEYWORDS double guf_rand_normal_sample_one_f64(guf_randstate *state, double mean, double std_dev);
|
||||
GUF_FN_KEYWORDS float guf_rand_normal_sample_one_f32(guf_randstate *state, float mean, float std_dev);
|
||||
|
||||
#endif
|
||||
|
||||
#if defined(GUF_IMPL) || defined(GUF_IMPL_STATIC)
|
||||
#include <math.h>
|
||||
#include <float.h>
|
||||
#include "guf_common.h"
|
||||
#include "guf_assert.h"
|
||||
#include "guf_math.h"
|
||||
|
||||
/*
|
||||
splitmix64 (public domain) written in 2015 by Sebastiano Vigna (vigna@acm.org)
|
||||
@ -62,11 +73,21 @@ GUF_FN_KEYWORDS uint64_t guf_rand_splitmix64(uint64_t *state)
|
||||
GUF_FN_KEYWORDS void guf_randstate_init(guf_randstate *state, uint64_t seed)
|
||||
{
|
||||
GUF_ASSERT_RELEASE(state);
|
||||
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
for (size_t i = 0; i < GUF_STATIC_BUF_SIZE(state->s); ++i) {
|
||||
state->s[i] = (uint32_t)(guf_rand_splitmix64(&seed) >> 32);
|
||||
}
|
||||
if (!state->s[0] && !state->s[1] && !state->s[2] && !state->s[3]) { // State must not be only zeroes:
|
||||
state->s[0] = 0x9e3779b9; // arbitrary constant != 0
|
||||
seed = 0x9e3779b97f4a7c15;
|
||||
for (size_t i = 1; i < GUF_STATIC_BUF_SIZE(state->s); ++i) {
|
||||
state->s[i] = (uint32_t)(guf_rand_splitmix64(&seed) >> 32);
|
||||
}
|
||||
}
|
||||
#else
|
||||
for (size_t i = 0; i < GUF_STATIC_BUF_SIZE(state->s); ++i) {
|
||||
state->s[i] = guf_rand_splitmix64(&seed);
|
||||
}
|
||||
|
||||
if (!state->s[0] && !state->s[1] && !state->s[2] && !state->s[3]) { // State must not be only zeroes:
|
||||
state->s[0] = 0x9e3779b97f4a7c15; // arbitrary constant != 0
|
||||
seed = state->s[0];
|
||||
@ -74,51 +95,97 @@ GUF_FN_KEYWORDS void guf_randstate_init(guf_randstate *state, uint64_t seed)
|
||||
state->s[i] = guf_rand_splitmix64(&seed);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
/*
|
||||
xoshiro256** 1.0 (public domain) written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org)
|
||||
cf. https://prng.di.unimi.it/xoshiro256starstar.c (last-retrieved 2025-02-11)
|
||||
*/
|
||||
GUF_FN_KEYWORDS uint32_t guf_rand_u32(guf_randstate *state)
|
||||
{
|
||||
GUF_ASSERT(state);
|
||||
GUF_ASSERT(state->s[0] || state->s[1] || state->s[2] || state->s[3]);
|
||||
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
/*
|
||||
xoshiro128** 1.1 (public domain) written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org)
|
||||
cf. https://prng.di.unimi.it/xoshiro128starstar.c (last-retrieved 2025-02-11)
|
||||
*/
|
||||
const uint32_t result = guf_rotl_u32(state->s[1] * 5, 7) * 9;
|
||||
const uint32_t t = state->s[1] << 9;
|
||||
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] = guf_rotl_u32(state->s[3], 11);
|
||||
return result;
|
||||
#else
|
||||
return (uint32_t)(guf_rand_u64(state) >> 32);
|
||||
#endif
|
||||
}
|
||||
|
||||
GUF_FN_KEYWORDS uint64_t guf_rand_u64(guf_randstate *state)
|
||||
{
|
||||
GUF_ASSERT(state);
|
||||
GUF_ASSERT(state->s[0] || state->s[1] || state->s[2] || state->s[3]);
|
||||
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
const uint32_t lower_bits = guf_rand_u32(state);
|
||||
const uint32_t upper_bits = guf_rand_u32(state);
|
||||
return ((uint64_t)upper_bits << 32) | (uint64_t)lower_bits; // TODO: not sure if that's a good idea...
|
||||
#else
|
||||
/*
|
||||
xoshiro256** 1.0 (public domain) written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org)
|
||||
cf. https://prng.di.unimi.it/xoshiro256starstar.c (last-retrieved 2025-02-11)
|
||||
*/
|
||||
const uint64_t result = guf_rotl_u64(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] = guf_rotl_u64(state->s[3], 45);
|
||||
|
||||
return result;
|
||||
#endif
|
||||
}
|
||||
|
||||
/*
|
||||
Equivalent to 2^128 calls to guf_rand_u64(); it can be used to generate 2^128
|
||||
Equivalent to 2^128 calls to guf_rand() (or 2^64 calls if GUF_RAND_32_BIT); it can be used to generate 2^128 (or 2^64)
|
||||
non-overlapping subsequences for parallel computations.
|
||||
*/
|
||||
void guf_randstate_jump(guf_randstate *state)
|
||||
{
|
||||
GUF_ASSERT(state);
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
static const uint32_t JUMP[] = { 0x8764000b, 0xf542d2d3, 0x6fa035c3, 0x77f2db5b };
|
||||
uint32_t s0 = 0;
|
||||
uint32_t s1 = 0;
|
||||
uint32_t s2 = 0;
|
||||
uint32_t s3 = 0;
|
||||
for(size_t i = 0; i < sizeof JUMP / sizeof *JUMP; ++i) {
|
||||
for(int b = 0; b < 32; ++b) {
|
||||
if (JUMP[i] & UINT32_C(1) << b) {
|
||||
s0 ^= state->s[0];
|
||||
s1 ^= state->s[1];
|
||||
s2 ^= state->s[2];
|
||||
s3 ^= state->s[3];
|
||||
}
|
||||
guf_rand_u32(state);
|
||||
}
|
||||
}
|
||||
state->s[0] = s0;
|
||||
state->s[1] = s1;
|
||||
state->s[2] = s2;
|
||||
state->s[3] = s3;
|
||||
#else
|
||||
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 (size_t 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];
|
||||
if (JUMP[i] & UINT64_C(1) << b) {
|
||||
s0 ^= state->s[0];
|
||||
s1 ^= state->s[1];
|
||||
s2 ^= state->s[2];
|
||||
s3 ^= state->s[3];
|
||||
@ -130,6 +197,7 @@ void guf_randstate_jump(guf_randstate *state)
|
||||
state->s[1] = s1;
|
||||
state->s[2] = s2;
|
||||
state->s[3] = s3;
|
||||
#endif
|
||||
}
|
||||
|
||||
// Generate double in the unit interval [0, 1)
|
||||
@ -142,10 +210,20 @@ GUF_FN_KEYWORDS double guf_rand_f64(guf_randstate *state)
|
||||
// Generate float in the unit interval [0, 1)
|
||||
GUF_FN_KEYWORDS float guf_rand_f32(guf_randstate *state)
|
||||
{
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
return (guf_rand_u32(state) >> 8) * 0x1.0p-24f; // 8 == 32 - 24; (float has a 24-bit mantissa/significand)
|
||||
#else
|
||||
return (guf_rand_u64(state) >> 40) * 0x1.0p-24f; // 40 == 64 - 24; (float has a 24-bit mantissa/significand)
|
||||
#endif
|
||||
}
|
||||
|
||||
GUF_FN_KEYWORDS bool guf_rand_bernoulli_trial(guf_randstate *state, double p)
|
||||
GUF_FN_KEYWORDS bool guf_rand_bernoulli_trial_f32(guf_randstate *state, float p)
|
||||
{
|
||||
p = guf_clamp_f32(p, 0, 1);
|
||||
return guf_rand_f32(state) < p; // never true for p = 0, always true for p = 1 since guf_rand_f64 is in range [0, 1)
|
||||
}
|
||||
|
||||
GUF_FN_KEYWORDS bool guf_rand_bernoulli_trial_f64(guf_randstate *state, double p)
|
||||
{
|
||||
p = guf_clamp_f64(p, 0, 1);
|
||||
return guf_rand_f64(state) < p; // never true for p = 0, always true for p = 1 since guf_rand_f64 is in range [0, 1)
|
||||
@ -153,17 +231,25 @@ GUF_FN_KEYWORDS bool guf_rand_bernoulli_trial(guf_randstate *state, double p)
|
||||
|
||||
GUF_FN_KEYWORDS bool guf_rand_flip(guf_randstate *state)
|
||||
{
|
||||
return guf_rand_bernoulli_trial(state, 0.5);
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
return guf_rand_bernoulli_trial_f32(state, 0.5f);
|
||||
#else
|
||||
return guf_rand_bernoulli_trial_f64(state, 0.5);
|
||||
#endif
|
||||
}
|
||||
|
||||
// returns uniformly-distributed random double in range [min, end) (or min if min == end)
|
||||
GUF_FN_KEYWORDS double guf_randrange_f64(guf_randstate *state, double min, double end)
|
||||
{
|
||||
if (min == (double)INFINITY) {
|
||||
min = DBL_MAX;
|
||||
} else if (min == (double)-INFINITY) {
|
||||
min = -DBL_MAX;
|
||||
}
|
||||
if (end == (double)INFINITY) {
|
||||
end = DBL_MAX;
|
||||
}
|
||||
if (min == (double)-INFINITY) {
|
||||
min = -DBL_MAX;
|
||||
} else if (end == (double)-INFINITY) {
|
||||
end = -DBL_MAX;
|
||||
}
|
||||
GUF_ASSERT_RELEASE(end >= min);
|
||||
return guf_rand_f64(state) * (end - min) + min;
|
||||
@ -172,11 +258,15 @@ GUF_FN_KEYWORDS double guf_randrange_f64(guf_randstate *state, double min, doubl
|
||||
// returns uniformly-distributed random float in range [min, end) (or min if min == end)
|
||||
GUF_FN_KEYWORDS float guf_randrange_f32(guf_randstate *state, float min, float end)
|
||||
{
|
||||
if (min == INFINITY) {
|
||||
min = FLT_MAX;
|
||||
} else if (min == -INFINITY) {
|
||||
min = -FLT_MAX;
|
||||
}
|
||||
if (end == INFINITY) {
|
||||
end = FLT_MAX;
|
||||
}
|
||||
if (min == -INFINITY) {
|
||||
min = -FLT_MAX;
|
||||
} else if (end == -INFINITY) {
|
||||
end = -FLT_MAX;
|
||||
}
|
||||
GUF_ASSERT_RELEASE(end >= min);
|
||||
return guf_rand_f32(state) * (end - min) + min;
|
||||
@ -189,8 +279,7 @@ GUF_FN_KEYWORDS int32_t guf_randrange_i32(guf_randstate *state, int32_t min, int
|
||||
if (min == max) {
|
||||
return min;
|
||||
}
|
||||
const double delta = (int64_t)max - (int64_t)min; // Cast to int64_t to avoid overflow.
|
||||
|
||||
const double delta = (double)max - (double)min;
|
||||
// cf. https://developer.mozilla.org/en-US/docs/Web/JavaScript/Reference/Global_Objects/Math/random (last-retrieved 2025-02-12)
|
||||
const double result = floor(guf_rand_f64(state) * (delta + 1.0) + min);
|
||||
GUF_ASSERT(result >= min && result <= max);
|
||||
@ -203,14 +292,12 @@ GUF_FN_KEYWORDS uint32_t guf_randrange_u32(guf_randstate *state, uint32_t min, u
|
||||
if (min == max) {
|
||||
return min;
|
||||
}
|
||||
const double delta = max - min; // Cannot overflow here.
|
||||
|
||||
const double delta = (double)max - (double)min;
|
||||
const double result = floor(guf_rand_f64(state) * (delta + 1.0) + min);
|
||||
GUF_ASSERT(result >= min && result <= max);
|
||||
return (uint32_t)result;
|
||||
}
|
||||
|
||||
|
||||
// returns uniformly-distributed random int64_t in range [min, max] (max is inclusive as opposed to the f32/f64 versions)
|
||||
GUF_FN_KEYWORDS int64_t guf_randrange_i64(guf_randstate *state, int64_t min, int64_t max)
|
||||
{
|
||||
@ -219,12 +306,10 @@ GUF_FN_KEYWORDS int64_t guf_randrange_i64(guf_randstate *state, int64_t min, int
|
||||
return min;
|
||||
}
|
||||
|
||||
// rand_max is 2^63 - 1 for rand_max_shift == 1
|
||||
const unsigned rand_max_shift = 1;
|
||||
const uint64_t rand_max = GUF_RAND_MAX >> rand_max_shift; // 2^63 - 1
|
||||
const uint64_t rand_max_i64 = UINT64_MAX >> 1; // 2^63 - 1 (== INT64_MAX)
|
||||
|
||||
const uint64_t delta = guf_absdiff_i64(max, min);
|
||||
if (delta > rand_max) {
|
||||
if (delta > rand_max_i64) {
|
||||
guf_panic(GUF_ERR_INT_OVERFLOW, GUF_ERR_MSG("in function guf_randrange_i64: interval [min, max] larger than INT64_MAX"));
|
||||
return -1;
|
||||
}
|
||||
@ -234,7 +319,7 @@ GUF_FN_KEYWORDS int64_t guf_randrange_i64(guf_randstate *state, int64_t min, int
|
||||
cf. https://c-faq.com/lib/randrange.html (last-retrieved 2025-02-11)
|
||||
https://stackoverflow.com/a/6852396 (last-retrieved 2025-02-11)
|
||||
*/
|
||||
const uint64_t num_rand_vals = rand_max + 1u; // 2^63
|
||||
const uint64_t num_rand_vals = rand_max_i64 + 1u; // 2^63
|
||||
const uint64_t num_bins = (delta + 1u);
|
||||
|
||||
const uint64_t bin_size = num_rand_vals / num_bins; // bin_size = floor(num_rand_vals / num_bins)
|
||||
@ -247,7 +332,7 @@ GUF_FN_KEYWORDS int64_t guf_randrange_i64(guf_randstate *state, int64_t min, int
|
||||
*/
|
||||
uint64_t step;
|
||||
do {
|
||||
step = guf_rand_u64(state) >> rand_max_shift;
|
||||
step = guf_rand_u64(state) >> 1; // [0, 2^63 - 1]
|
||||
} while (step >= limit);
|
||||
step = step / bin_size;
|
||||
|
||||
@ -319,7 +404,6 @@ GUF_FN_KEYWORDS float guf_rand_normal_sample_one_f32(guf_randstate *state, float
|
||||
|
||||
#undef GUF_IMPL
|
||||
#undef GUF_IMPL_STATIC
|
||||
#endif
|
||||
#endif /* endif GUF_IMPL/GUF_IMPL_STATIC */
|
||||
|
||||
#undef GUF_FN_KEYWORDS
|
||||
#endif
|
||||
|
||||
239
src/guf_rand32.h
239
src/guf_rand32.h
@ -1,239 +0,0 @@
|
||||
#ifndef GUF_RAND32_H
|
||||
#define GUF_RAND32_H
|
||||
#include "guf_common.h"
|
||||
#include "guf_assert.h"
|
||||
#include "guf_math.h"
|
||||
#include <math.h>
|
||||
#include <float.h>
|
||||
|
||||
#define GUF_RAND32_MAX UINT32_MAX
|
||||
|
||||
typedef struct guf_randstate32 { // State for xoshiro128** 1.1
|
||||
uint32_t s[4];
|
||||
} guf_randstate32;
|
||||
|
||||
#ifdef GUF_IMPL_STATIC
|
||||
#define GUF_FN_KEYWORDS static
|
||||
#else
|
||||
#define GUF_FN_KEYWORDS
|
||||
#endif
|
||||
|
||||
GUF_FN_KEYWORDS uint64_t guf_rand32_splitmix64(uint64_t *state);
|
||||
|
||||
GUF_FN_KEYWORDS void guf_randstate32_init(guf_randstate32 *state, uint64_t seed);
|
||||
void guf_randstate32_jump(guf_randstate32 *state); // Advance the state; equivalent to 2^128 calls to guf_rand32_u32(state)
|
||||
|
||||
// uniform distributions using xoshiro128** 1.1
|
||||
GUF_FN_KEYWORDS uint32_t guf_rand32_u32(guf_randstate32 *state); // [0, GUF_RAND_MAX]
|
||||
GUF_FN_KEYWORDS float guf_rand32_f32(guf_randstate32 *state); // [0.f, 1.f)
|
||||
|
||||
// return true with a probability of p, false with a probability of (1 - p)
|
||||
GUF_FN_KEYWORDS bool guf_rand32_bernoulli_trial(guf_randstate32 *state, float p);
|
||||
GUF_FN_KEYWORDS bool guf_rand32_flip(guf_randstate32 *state); // Fair coin flip (bernoulli trial with p == 0.5)
|
||||
|
||||
GUF_FN_KEYWORDS float guf_rand32range_f32(guf_randstate32 *state, float min, float end); // [min, end)
|
||||
GUF_FN_KEYWORDS int32_t guf_rand32range_i32(guf_randstate32 *state, int32_t min, int32_t max); // [min, max]
|
||||
|
||||
// normal distributions
|
||||
GUF_FN_KEYWORDS void guf_rand32_normal_sample_f32(guf_randstate32 *state, float mean, float std_dev, float *result, ptrdiff_t n);
|
||||
GUF_FN_KEYWORDS float guf_rand32_normal_sample_one_f32(guf_randstate32 *state, float mean, float std_dev);
|
||||
|
||||
#if defined(GUF_IMPL) || defined(GUF_IMPL_STATIC)
|
||||
|
||||
/*
|
||||
splitmix64 (public domain) written in 2015 by Sebastiano Vigna (vigna@acm.org)
|
||||
cf. https://prng.di.unimi.it/splitmix64.c (last-retrieved 2025-02-11)
|
||||
*/
|
||||
GUF_FN_KEYWORDS uint64_t guf_rand32_splitmix64(uint64_t *state)
|
||||
{
|
||||
GUF_ASSERT(state);
|
||||
uint64_t z = ((*state) += 0x9e3779b97f4a7c15);
|
||||
z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9;
|
||||
z = (z ^ (z >> 27)) * 0x94d049bb133111eb;
|
||||
return z ^ (z >> 31);
|
||||
}
|
||||
|
||||
GUF_FN_KEYWORDS void guf_randstate32_init(guf_randstate32 *state, uint64_t seed)
|
||||
{
|
||||
GUF_ASSERT_RELEASE(state);
|
||||
|
||||
uint64_t split = guf_rand32_splitmix64(&seed);
|
||||
state->s[0] = (uint32_t)split; // lower 32-bits
|
||||
state->s[1] = (uint32_t)(split >> 32); // upper 32-bits
|
||||
split = guf_rand32_splitmix64(&seed);
|
||||
state->s[2] = (uint32_t)split; // lower 32-bits
|
||||
state->s[3] = (uint32_t)(split >> 32); // upper 32-bits
|
||||
|
||||
if (!state->s[0] && !state->s[1] && !state->s[2] && !state->s[3]) { // State must not be only zeroes:
|
||||
state->s[0] = 0x9e3779b9; // arbitrary constant != 0
|
||||
seed = 0x9e3779b97f4a7c15;
|
||||
split = guf_rand32_splitmix64(&seed);
|
||||
state->s[0] = (uint32_t)split; // lower 32-bits
|
||||
state->s[1] = (uint32_t)(split >> 32); // upper 32-bits
|
||||
split = guf_rand32_splitmix64(&seed);
|
||||
state->s[2] = (uint32_t)split; // lower 32-bits
|
||||
state->s[3] = (uint32_t)(split >> 32); // upper 32-bits
|
||||
}
|
||||
}
|
||||
|
||||
/*
|
||||
xoshiro128** 1.1 (public domain) written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org)
|
||||
cf. https://prng.di.unimi.it/xoshiro128starstar.c (last-retrieved 2025-02-11)
|
||||
*/
|
||||
|
||||
GUF_FN_KEYWORDS uint32_t guf_rand32_u32(guf_randstate32 *state)
|
||||
{
|
||||
GUF_ASSERT(state);
|
||||
GUF_ASSERT(state->s[0] || state->s[1] || state->s[2] || state->s[3]);
|
||||
|
||||
const uint32_t result = guf_rotl_u32(state->s[1] * 5, 7) * 9;
|
||||
|
||||
const uint32_t t = state->s[1] << 9;
|
||||
|
||||
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] = guf_rotl_u32(state->s[3], 11);
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
/*
|
||||
Equivalent to 2^128 calls to guf_rand_u64(); it can be used to generate 2^128
|
||||
non-overlapping subsequences for parallel computations.
|
||||
*/
|
||||
void guf_randstate32_jump(guf_randstate32 *state)
|
||||
{
|
||||
GUF_ASSERT(state);
|
||||
static const uint32_t JUMP[] = { 0x8764000b, 0xf542d2d3, 0x6fa035c3, 0x77f2db5b };
|
||||
uint32_t s0 = 0;
|
||||
uint32_t s1 = 0;
|
||||
uint32_t s2 = 0;
|
||||
uint32_t s3 = 0;
|
||||
for(size_t i = 0; i < sizeof JUMP / sizeof *JUMP; i++) {
|
||||
for(int b = 0; b < 32; b++) {
|
||||
if (JUMP[i] & UINT32_C(1) << b) {
|
||||
s0 ^= state->s[0];
|
||||
s1 ^= state->s[1];
|
||||
s2 ^= state->s[2];
|
||||
s3 ^= state->s[3];
|
||||
}
|
||||
guf_rand32_u32(state);
|
||||
}
|
||||
}
|
||||
state->s[0] = s0;
|
||||
state->s[1] = s1;
|
||||
state->s[2] = s2;
|
||||
state->s[3] = s3;
|
||||
}
|
||||
|
||||
// Generate float in the unit interval [0, 1)
|
||||
GUF_FN_KEYWORDS float guf_rand32_f32(guf_randstate32 *state)
|
||||
{
|
||||
return (guf_rand32_u32(state) >> 8) * 0x1.0p-24f; // 8 == 32 - 24; (float has a 24-bit mantissa/significand)
|
||||
}
|
||||
|
||||
GUF_FN_KEYWORDS bool guf_rand32_bernoulli_trial(guf_randstate32 *state, float p)
|
||||
{
|
||||
p = guf_clamp_f32(p, 0, 1);
|
||||
return guf_rand32_f32(state) < p; // never true for p = 0, always true for p = 1 since guf_rand32_f32 is in range [0, 1)
|
||||
}
|
||||
|
||||
GUF_FN_KEYWORDS bool guf_rand32_flip(guf_randstate32 *state)
|
||||
{
|
||||
return guf_rand32_bernoulli_trial(state, 0.5f);
|
||||
}
|
||||
|
||||
// returns uniformly-distributed random float in range [min, end) (or min if min == end)
|
||||
GUF_FN_KEYWORDS float guf_rand32range_f32(guf_randstate32 *state, float min, float end)
|
||||
{
|
||||
if (end == INFINITY) {
|
||||
end = FLT_MAX;
|
||||
}
|
||||
if (min == -INFINITY) {
|
||||
min = -FLT_MAX;
|
||||
}
|
||||
GUF_ASSERT_RELEASE(end >= min);
|
||||
return guf_rand32_f32(state) * (end - min) + min;
|
||||
}
|
||||
|
||||
// returns uniformly-distributed random int32_t in range [min, max] (max is inclusive as opposed to the f32/f64 versions)
|
||||
GUF_FN_KEYWORDS int32_t guf_rand32range_i32(guf_randstate32 *state, int32_t min, int32_t max)
|
||||
{
|
||||
GUF_ASSERT_RELEASE(max >= min);
|
||||
if (min == max) {
|
||||
return min;
|
||||
}
|
||||
|
||||
// rand_max is 2^32 - 1 for rand_max_shift == 1
|
||||
const unsigned rand_max_shift = 1;
|
||||
const uint32_t rand_max = GUF_RAND32_MAX >> rand_max_shift; // 2^32 - 1
|
||||
|
||||
const uint32_t delta = guf_absdiff_i32(max, min);
|
||||
if (delta > rand_max) {
|
||||
guf_panic(GUF_ERR_INT_OVERFLOW, GUF_ERR_MSG("in function guf_randrange32_i32: interval [min, max] larger than INT32_MAX"));
|
||||
return -1;
|
||||
}
|
||||
|
||||
/*
|
||||
cf. https://c-faq.com/lib/randrange.html (last-retrieved 2025-02-11)
|
||||
https://stackoverflow.com/a/6852396 (last-retrieved 2025-02-11)
|
||||
*/
|
||||
const uint32_t num_rand_vals = rand_max + 1u; // 2^31
|
||||
const uint32_t num_bins = (delta + 1u);
|
||||
|
||||
const uint32_t bin_size = num_rand_vals / num_bins; // bin_size = floor(num_rand_vals / num_bins)
|
||||
const uint32_t limit = num_rand_vals - (num_rand_vals % num_bins); // limit == bin_size * num_bins
|
||||
GUF_ASSERT(limit == bin_size * num_bins);
|
||||
uint32_t step;
|
||||
do {
|
||||
step = guf_rand32_u32(state) >> rand_max_shift;
|
||||
} while (step >= limit);
|
||||
step = step / bin_size;
|
||||
|
||||
const int32_t rnd = min + step;
|
||||
GUF_ASSERT(rnd >= min && rnd <= max);
|
||||
return rnd;
|
||||
}
|
||||
|
||||
// Box-Müller-transform transcribed from wikipedia, cf. https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform (last-retrieved 2025-02-12)
|
||||
|
||||
GUF_FN_KEYWORDS void guf_rand32_normal_sample_f32(guf_randstate32 *state, float mean, float std_dev, float *result, ptrdiff_t n)
|
||||
{
|
||||
GUF_ASSERT_RELEASE(result);
|
||||
GUF_ASSERT_RELEASE(n >= 0);
|
||||
const float TAU = 2.f * (float)GUF_PI;
|
||||
|
||||
ptrdiff_t i = 0;
|
||||
while (i < n) {
|
||||
float u1, u2;
|
||||
do {
|
||||
u1 = guf_rand32_f32(state);
|
||||
} while (u1 == 0);
|
||||
u2 = guf_rand32_f32(state);
|
||||
|
||||
const float mag = std_dev * sqrtf(-2.f * logf(u1));
|
||||
result[i++] = mag * cosf(TAU * u2) + mean;
|
||||
if (i < n) {
|
||||
result[i++] = mag * sinf(TAU * u2) + mean;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
GUF_FN_KEYWORDS float guf_rand32_normal_sample_one_f32(guf_randstate32 *state, float mean, float std_dev)
|
||||
{
|
||||
float result;
|
||||
guf_rand32_normal_sample_f32(state, mean, std_dev, &result, 1);
|
||||
return result;
|
||||
}
|
||||
|
||||
#undef GUF_IMPL
|
||||
#undef GUF_IMPL_STATIC
|
||||
#endif
|
||||
|
||||
#undef GUF_FN_KEYWORDS
|
||||
#endif
|
||||
@ -43,12 +43,10 @@
|
||||
#define GUF_IMPL_STATIC
|
||||
#include "guf_dbuf.h"
|
||||
|
||||
#define GUF_RAND_32_BIT
|
||||
#define GUF_IMPL_STATIC
|
||||
#include "guf_rand.h"
|
||||
|
||||
#define GUF_IMPL_STATIC
|
||||
#include "guf_rand32.h"
|
||||
|
||||
int main(void)
|
||||
{
|
||||
printf("libguf test: compiled with C %ld\n", __STDC_VERSION__);
|
||||
@ -165,18 +163,13 @@ int main(void)
|
||||
}
|
||||
|
||||
dbuf_int_free(&integers, NULL);
|
||||
|
||||
printf("\n");
|
||||
guf_randstate rng;
|
||||
guf_randstate_init(&rng, time(NULL));
|
||||
|
||||
guf_randstate32 rng32;
|
||||
guf_randstate32_init(&rng32, time(NULL));
|
||||
|
||||
printf("\n");
|
||||
int heads = 0, tails = 0;
|
||||
int throws = 10;
|
||||
for (i = 0; i < throws; ++i) {
|
||||
bool is_head = guf_rand32_flip(&rng32);
|
||||
bool is_head = guf_rand_flip(&rng);
|
||||
if (is_head) {
|
||||
puts("head");
|
||||
++heads;
|
||||
@ -189,13 +182,12 @@ int main(void)
|
||||
|
||||
int result[256];
|
||||
memset(result, 0, sizeof result);
|
||||
for (int n = 0; n < 24000; ++n) {
|
||||
double r = round(guf_rand_normal_sample_one_f64(&rng, 100, 10));
|
||||
if (r >= 0 && r < GUF_STATIC_BUF_SIZE(result)) {
|
||||
result[(int)r] += 1;
|
||||
}
|
||||
for (int n = 0; n < 32000; ++n) {
|
||||
float r = roundf(guf_rand_normal_sample_one_f32(&rng, 100, 15));
|
||||
r = guf_clamp_f32(r, 0, 255);
|
||||
result[(int)r] += 1;
|
||||
}
|
||||
for (size_t n = 50; n <= 150; ++n) {
|
||||
for (size_t n = 60; n <= 140; ++n) {
|
||||
printf("%zu:\t", n);
|
||||
for (int j = 0; j < result[n] / 8; ++j) {
|
||||
putc('#', stdout);
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user