Add separate 32/64 bit versions to guf_rand
This commit is contained in:
parent
05f995e855
commit
d062784425
@ -28,14 +28,31 @@
|
||||
|
||||
#if GUF_PLATFORM_BITS <= 32
|
||||
#define GUF_HASH_32_BIT
|
||||
#define GUF_RAND_32_BIT
|
||||
#endif
|
||||
|
||||
#if (defined(__STDC_VERSION__) && __STDC_VERSION__ >= 201112L)
|
||||
#if defined(__STDC_VERSION__)
|
||||
#if __STDC_VERSION__ >= 199901L
|
||||
#define GUF_STDC_AT_LEAST_C99
|
||||
#else
|
||||
#error "libguf only supports C99 and above"
|
||||
#endif
|
||||
|
||||
#if __STDC_VERSION__ >= 201112L
|
||||
#define GUF_STDC_AT_LEAST_C11
|
||||
#endif
|
||||
#if __STDC_VERSION__ >= 201710L
|
||||
#define GUF_STDC_AT_LEAST_C17
|
||||
#endif
|
||||
#if __STDC_VERSION__ >= 202311L
|
||||
#define GUF_STDC_AT_LEAST_C23
|
||||
#endif
|
||||
#endif
|
||||
|
||||
#if (defined(__cplusplus) && __cplusplus >= 201103L)
|
||||
#if defined(__cplusplus)
|
||||
#if __cplusplus >= 201103L
|
||||
#define GUF_STDCPP_AT_LEAST_CPP11
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/*
|
||||
|
||||
@ -19,10 +19,14 @@
|
||||
*/
|
||||
|
||||
#define GUF_HASH32_INIT UINT32_C(2166136261)
|
||||
#define GUF_HASH64_INIT UINT64_C(14695981039346656037)
|
||||
#ifdef UINT64_MAX
|
||||
#define GUF_HASH64_INIT UINT64_C(14695981039346656037)
|
||||
#endif
|
||||
|
||||
GUF_HASH_KWRDS uint32_t guf_hash32(const void *data, ptrdiff_t num_bytes, uint32_t hash); // FNV-1a (32 bit)
|
||||
GUF_HASH_KWRDS uint64_t guf_hash64(const void *data, ptrdiff_t num_bytes, uint64_t hash); // FNV-1a (64 bit)
|
||||
#ifdef UINT64_MAX
|
||||
GUF_HASH_KWRDS uint64_t guf_hash64(const void *data, ptrdiff_t num_bytes, uint64_t hash); // FNV-1a (64 bit)
|
||||
#endif
|
||||
|
||||
#ifdef GUF_HASH_32_BIT
|
||||
typedef uint32_t guf_hash_size_t;
|
||||
@ -32,6 +36,9 @@ GUF_HASH_KWRDS uint64_t guf_hash64(const void *data, ptrdiff_t num_bytes, uint64
|
||||
return guf_hash32(data, num_bytes, hash);
|
||||
}
|
||||
#else
|
||||
#ifndef UINT64_MAX
|
||||
#error "guf_hash.h: Platform does not support uint64_t (define GUF_HASH_32_BIT to fix)"
|
||||
#endif
|
||||
typedef uint64_t guf_hash_size_t;
|
||||
#define GUF_HASH_INIT GUF_HASH64_INIT
|
||||
#define GUF_HASH_MAX UINT64_MAX
|
||||
@ -59,8 +66,9 @@ GUF_HASH_KWRDS uint32_t guf_hash32(const void *data, ptrdiff_t num_bytes, uint32
|
||||
return hash;
|
||||
}
|
||||
|
||||
GUF_HASH_KWRDS uint64_t guf_hash64(const void *data, ptrdiff_t num_bytes, uint64_t hash)
|
||||
{
|
||||
#ifdef UINT64_MAX
|
||||
GUF_HASH_KWRDS 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...
|
||||
@ -70,7 +78,8 @@ GUF_HASH_KWRDS uint64_t guf_hash64(const void *data, ptrdiff_t num_bytes, uint64
|
||||
hash *= FNV_64_PRIME;
|
||||
}
|
||||
return hash;
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
#undef GUF_HASH_IMPL
|
||||
#undef GUF_HASH_IMPL_STATIC
|
||||
|
||||
@ -124,6 +124,17 @@ static inline bool guf_sub_is_overflow_i32(int32_t a, int32_t b)
|
||||
(b > 0 && a < INT32_MIN + b); // a - b underflow
|
||||
}
|
||||
|
||||
static inline bool guf_add_is_overflow_i64(int64_t a, int64_t b)
|
||||
{
|
||||
return (b > 0 && a > INT64_MAX - b) || // a + b overflow
|
||||
(b < 0 && a < INT64_MIN - b); // a + b underflow
|
||||
}
|
||||
static inline bool guf_sub_is_overflow_i64(int64_t a, int64_t b)
|
||||
{
|
||||
return (b < 0 && a > INT64_MAX + b) || // a - b overflow
|
||||
(b > 0 && a < INT64_MIN + b); // a - b underflow
|
||||
}
|
||||
|
||||
static inline bool guf_size_calc_safe(ptrdiff_t count, ptrdiff_t sizeof_elem, ptrdiff_t *result)
|
||||
{
|
||||
if (count < 0 || sizeof_elem <= 0) {
|
||||
|
||||
624
src/guf_rand.h
624
src/guf_rand.h
@ -1,6 +1,5 @@
|
||||
/*
|
||||
is parametrized: yes
|
||||
TODO: - maybe allow 64- and 32-bit implementations to co-exist...
|
||||
*/
|
||||
|
||||
#if defined(GUF_RAND_IMPL_STATIC)
|
||||
@ -14,41 +13,117 @@
|
||||
|
||||
#include "guf_common.h"
|
||||
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
#define GUF_RAND_MAX UINT32_MAX
|
||||
typedef struct guf_randstate { // State for xoshiro128** 1.1
|
||||
typedef struct guf_rand32_state {
|
||||
uint32_t s[4];
|
||||
} guf_randstate;
|
||||
#else
|
||||
#define GUF_RAND_MAX UINT64_MAX
|
||||
typedef struct guf_randstate { // State for xoshiro256** 1.0
|
||||
} guf_rand32_state;
|
||||
|
||||
#ifdef UINT64_MAX
|
||||
typedef struct guf_rand64_state {
|
||||
uint64_t s[4];
|
||||
} guf_randstate;
|
||||
} guf_rand64_state;
|
||||
#endif
|
||||
|
||||
GUF_RAND_KWRDS uint64_t guf_rand_splitmix64(uint64_t *state);
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
#define GUF_RAND_MAX UINT32_MAX
|
||||
typedef guf_rand32_state guf_randstate;
|
||||
typedef uint32_t guf_rand_seed_t;
|
||||
#else
|
||||
#ifndef UINT64_MAX
|
||||
#error "guf_rand.h: Platform does not support uint64_t (define GUF_RAND_32_BIT to fix)"
|
||||
#endif
|
||||
#define GUF_RAND_MAX UINT64_MAX
|
||||
typedef guf_rand64_state guf_randstate;
|
||||
typedef uint64_t guf_rand_seed_t;
|
||||
#endif
|
||||
|
||||
#ifdef UINT64_MAX
|
||||
GUF_RAND_KWRDS uint64_t guf_rand_splitmix64(uint64_t *state);
|
||||
#endif
|
||||
GUF_RAND_KWRDS uint32_t guf_rand_splitmix32(uint32_t *state);
|
||||
|
||||
GUF_RAND_KWRDS void guf_randstate_init(guf_randstate *state, uint64_t seed);
|
||||
GUF_RAND_KWRDS void guf_randstate_jump(guf_randstate *state); // Advance the state; equivalent to 2^128 calls to guf_rand_u64(state)
|
||||
GUF_RAND_KWRDS void guf_randstate_init(guf_randstate *state, guf_rand_seed_t seed);
|
||||
#ifdef UINT64_MAX
|
||||
GUF_RAND_KWRDS void guf_rand64_state_init(guf_rand64_state *state, uint64_t seed);
|
||||
#endif
|
||||
GUF_RAND_KWRDS void guf_rand32_state_init(guf_rand32_state *state, uint32_t seed);
|
||||
|
||||
GUF_RAND_KWRDS void guf_randstate_jump(guf_randstate *state); // Advance the state; equivalent to 2^128 (or 2^64) calls to guf_rand(state)
|
||||
#ifdef UINT64_MAX
|
||||
GUF_RAND_KWRDS void guf_rand64_state_jump(guf_rand64_state *state); // equivalent to 2^128 calls to guf_rand64_u64(state)
|
||||
#endif
|
||||
GUF_RAND_KWRDS void guf_rand32_state_jump(guf_rand32_state *state); // equivalent to 2^64 calls to guf_rand32_u32(state)
|
||||
|
||||
// Uniform distributions:
|
||||
|
||||
// uniform distributions
|
||||
GUF_RAND_KWRDS uint32_t guf_rand_u32(guf_randstate *state); // [0, UINT32_MAX]
|
||||
GUF_RAND_KWRDS uint64_t guf_rand_u64(guf_randstate *state); // [0, UINT64_MAX]
|
||||
GUF_RAND_KWRDS double guf_rand_f64(guf_randstate *state); // [0.0, 1.0)
|
||||
#ifdef UINT64_MAX
|
||||
// NOTE: slow on 32-bit platforms.
|
||||
GUF_RAND_KWRDS uint32_t guf_rand64_u32(guf_rand64_state *state);
|
||||
#endif
|
||||
GUF_RAND_KWRDS uint32_t guf_rand32_u32(guf_rand32_state *state);
|
||||
|
||||
#ifdef UINT64_MAX
|
||||
// NOTE: Slow on 32-bit platforms.
|
||||
GUF_RAND_KWRDS uint64_t guf_rand_u64(guf_randstate *state); // [0, UINT64_MAX]
|
||||
GUF_RAND_KWRDS uint64_t guf_rand32_u64(guf_rand32_state *state);
|
||||
GUF_RAND_KWRDS uint64_t guf_rand64_u64(guf_rand64_state *state);
|
||||
#else
|
||||
GUF_RAND_KWRDS uint_least64_t guf_rand_u64(guf_randstate *state) // NOTE: Slow on 32-bit platforms.
|
||||
GUF_RAND_KWRDS uint_least64_t guf_rand32_u64(guf_rand32_state *state); // NOTE: Slow on 32-bit platforms.
|
||||
#endif
|
||||
|
||||
GUF_RAND_KWRDS double guf_rand_f64(guf_randstate *state); // [0.0, 1.0) (NOTE: slow on 32-bit platforms.)
|
||||
#ifdef UINT64_MAX
|
||||
// NOTE: slow on 32-bit platforms.
|
||||
GUF_RAND_KWRDS double guf_rand64_f64(guf_rand64_state *state);
|
||||
#endif
|
||||
GUF_RAND_KWRDS double guf_rand32_f64(guf_rand32_state *state); // NOTE: slow on 32-bit platforms.
|
||||
|
||||
GUF_RAND_KWRDS float guf_rand_f32(guf_randstate *state); // [0.f, 1.f)
|
||||
#ifdef UINT64_MAX
|
||||
// NOTE: slow on 32-bit platforms.
|
||||
GUF_RAND_KWRDS float guf_rand64_f32(guf_rand64_state *state);
|
||||
#endif
|
||||
GUF_RAND_KWRDS float guf_rand32_f32(guf_rand32_state *state);
|
||||
|
||||
// return true with a probability of p, false with a probability of (1 - p)
|
||||
GUF_RAND_KWRDS bool guf_rand_bernoulli_trial_f32(guf_randstate *state, float p);
|
||||
GUF_RAND_KWRDS bool guf_rand_bernoulli_trial_f64(guf_randstate *state, double p);
|
||||
GUF_RAND_KWRDS bool guf_rand_bernoulli_trial_f64(guf_randstate *state, double p); // NOTE: slow on 32-bit platforms.
|
||||
GUF_RAND_KWRDS bool guf_rand_flip(guf_randstate *state); // Fair coin flip (bernoulli trial with p == 0.5)
|
||||
#ifdef UINT64_MAX
|
||||
// NOTE: slow on 32-bit platforms.
|
||||
GUF_RAND_KWRDS bool guf_rand64_bernoulli_trial_f32(guf_rand64_state *state, float p);
|
||||
GUF_RAND_KWRDS bool guf_rand64_bernoulli_trial_f64(guf_rand64_state *state, double p);
|
||||
GUF_RAND_KWRDS bool guf_rand64_flip(guf_rand64_state *state);
|
||||
#endif
|
||||
GUF_RAND_KWRDS bool guf_rand32_bernoulli_trial_f32(guf_rand32_state *state, float p);
|
||||
GUF_RAND_KWRDS bool guf_rand32_bernoulli_trial_f64(guf_rand32_state *state, double p); // NOTE: slow on 32-bit platforms.
|
||||
GUF_RAND_KWRDS bool guf_rand32_flip(guf_rand32_state *state);
|
||||
|
||||
GUF_RAND_KWRDS double guf_randrange_f64(guf_randstate *state, double min, double end); // [min, end)
|
||||
GUF_RAND_KWRDS double guf_randrange_f64(guf_randstate *state, double min, double end); // [min, end) (NOTE: slow on 32-bit platforms.)
|
||||
GUF_RAND_KWRDS float guf_randrange_f32(guf_randstate *state, float min, float end); // [min, end)
|
||||
#ifdef UINT64_MAX
|
||||
// (NOTE: slow on 32-bit platforms.)
|
||||
GUF_RAND_KWRDS double guf_rand64_range_f64(guf_rand64_state *state, double min, double end);
|
||||
GUF_RAND_KWRDS float guf_rand64_range_f32(guf_rand64_state *state, float min, float end);
|
||||
#endif
|
||||
GUF_RAND_KWRDS double guf_rand32_range_f64(guf_rand32_state *state, double min, double end); // (NOTE: slow on 32-bit platforms.)
|
||||
GUF_RAND_KWRDS float guf_rand32_range_f32(guf_rand32_state *state, float min, float end);
|
||||
|
||||
GUF_RAND_KWRDS int32_t guf_randrange_i32(guf_randstate *state, int32_t min, int32_t max); // [min, max]
|
||||
GUF_RAND_KWRDS uint32_t guf_randrange_u32(guf_randstate *state, uint32_t min, uint32_t max); // [min, max]
|
||||
GUF_RAND_KWRDS int64_t guf_randrange_i64(guf_randstate *state, int64_t min, int64_t max); // [min, max]
|
||||
GUF_RAND_KWRDS uint32_t guf_randrange_u32(guf_randstate *state, uint32_t min, uint32_t max); // [min, max] (NOTE: slow on 32-bit platforms.)
|
||||
#ifdef UINT64_MAX
|
||||
// (NOTE: slow on 32-bit platforms.)
|
||||
GUF_RAND_KWRDS int32_t guf_rand64_range_i32(guf_rand64_state *state, int32_t min, int32_t max);
|
||||
GUF_RAND_KWRDS uint32_t guf_rand64_range_u32(guf_rand64_state *state, uint32_t min, uint32_t max);
|
||||
#endif
|
||||
GUF_RAND_KWRDS int32_t guf_rand32_range_i32(guf_rand32_state *state, int32_t min, int32_t max);
|
||||
GUF_RAND_KWRDS uint32_t guf_rand32_range_u32(guf_rand32_state *state, uint32_t min, uint32_t max); // (NOTE: slow on 32-bit platforms.)
|
||||
|
||||
#ifdef UINT64_MAX
|
||||
// (NOTE: slow on 32-bit platforms.)
|
||||
GUF_RAND_KWRDS int64_t guf_randrange_i64(guf_randstate *state, int64_t min, int64_t max); // [min, max] (NOTE: slow on 32-bit platforms.)
|
||||
#endif
|
||||
|
||||
// normal distributions
|
||||
GUF_RAND_KWRDS void guf_rand_normal_sample_f64(guf_randstate *state, double mean, double std_dev, double *result, ptrdiff_t n);
|
||||
@ -65,18 +140,20 @@ GUF_RAND_KWRDS float guf_rand_normal_sample_one_f32(guf_randstate *state, float
|
||||
#include "guf_assert.h"
|
||||
#include "guf_math.h"
|
||||
|
||||
/*
|
||||
#ifdef UINT64_MAX
|
||||
/*
|
||||
splitmix64 written in 2015 by Sebastiano Vigna (vigna@acm.org) (released as public domain)
|
||||
cf. https://prng.di.unimi.it/splitmix64.c (last-retrieved 2025-02-11)
|
||||
*/
|
||||
GUF_RAND_KWRDS uint64_t guf_rand_splitmix64(uint64_t *state)
|
||||
{
|
||||
*/
|
||||
GUF_RAND_KWRDS uint64_t guf_rand_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);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
/*
|
||||
splitmix32 written in 2016 by Kaito Udagawa (released under CC0 <http://creativecommons.org/publicdomain/zero/1.0/>)
|
||||
@ -91,23 +168,23 @@ GUF_RAND_KWRDS uint32_t guf_rand_splitmix32(uint32_t *state)
|
||||
return z ^ (z >> 16);
|
||||
}
|
||||
|
||||
|
||||
GUF_RAND_KWRDS void guf_randstate_init(guf_randstate *state, uint64_t seed)
|
||||
GUF_RAND_KWRDS void guf_rand32_state_init(guf_rand32_state *state, uint32_t seed)
|
||||
{
|
||||
GUF_ASSERT_RELEASE(state);
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
for (size_t i = 0; i < GUF_ARR_SIZE(state->s); ++i) {
|
||||
state->s[i] = (uint32_t)(guf_rand_splitmix64(&seed) >> 32);
|
||||
state->s[i] = guf_rand_splitmix32(&seed);
|
||||
}
|
||||
|
||||
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;
|
||||
seed = state->s[0];
|
||||
for (size_t i = 1; i < GUF_ARR_SIZE(state->s); ++i) {
|
||||
state->s[i] = (uint32_t)(guf_rand_splitmix64(&seed) >> 32);
|
||||
state->s[i] = guf_rand_splitmix32(&seed);
|
||||
}
|
||||
}
|
||||
#else
|
||||
}
|
||||
|
||||
#ifdef UINT64_MAX
|
||||
GUF_RAND_KWRDS void guf_rand64_state_init(guf_rand64_state *state, uint64_t seed)
|
||||
{
|
||||
for (size_t i = 0; i < GUF_ARR_SIZE(state->s); ++i) {
|
||||
state->s[i] = guf_rand_splitmix64(&seed);
|
||||
}
|
||||
@ -118,15 +195,22 @@ GUF_RAND_KWRDS void guf_randstate_init(guf_randstate *state, uint64_t seed)
|
||||
state->s[i] = guf_rand_splitmix64(&seed);
|
||||
}
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
GUF_RAND_KWRDS void guf_randstate_init(guf_randstate *state, guf_rand_seed_t seed)
|
||||
{
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
guf_rand32_state_init(state, seed);
|
||||
#else
|
||||
guf_rand64_state_init(state, seed);
|
||||
#endif
|
||||
}
|
||||
|
||||
GUF_RAND_KWRDS uint32_t guf_rand_u32(guf_randstate *state)
|
||||
GUF_RAND_KWRDS uint32_t guf_rand32_u32(guf_rand32_state *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)
|
||||
@ -140,21 +224,32 @@ GUF_RAND_KWRDS uint32_t guf_rand_u32(guf_randstate *state)
|
||||
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_RAND_KWRDS uint64_t guf_rand_u64(guf_randstate *state)
|
||||
#ifdef UINT64_MAX
|
||||
GUF_RAND_KWRDS uint32_t guf_rand64_u32(guf_rand64_state *state)
|
||||
{
|
||||
return (uint32_t)(guf_rand64_u64(state) >> 32);
|
||||
}
|
||||
#endif
|
||||
|
||||
GUF_RAND_KWRDS 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
|
||||
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...
|
||||
return guf_rand32_u32(state);
|
||||
#else
|
||||
return guf_rand64_u32(state);
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
#ifdef UINT64_MAX
|
||||
GUF_RAND_KWRDS uint64_t guf_rand64_u64(guf_rand64_state *state)
|
||||
{
|
||||
GUF_ASSERT(state);
|
||||
GUF_ASSERT(state->s[0] || state->s[1] || state->s[2] || state->s[3]);
|
||||
/*
|
||||
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)
|
||||
@ -168,38 +263,78 @@ GUF_RAND_KWRDS uint64_t guf_rand_u64(guf_randstate *state)
|
||||
state->s[2] ^= t;
|
||||
state->s[3] = guf_rotl_u64(state->s[3], 45);
|
||||
return result;
|
||||
}
|
||||
#endif
|
||||
|
||||
#ifdef UINT64_MAX
|
||||
GUF_RAND_KWRDS uint64_t guf_rand32_u64(guf_rand32_state *state)
|
||||
#else
|
||||
GUF_RAND_KWRDS uint_least64_t guf_rand32_u64(guf_rand32_state *state)
|
||||
#endif
|
||||
{
|
||||
GUF_ASSERT(state);
|
||||
GUF_ASSERT(state->s[0] || state->s[1] || state->s[2] || state->s[3]);
|
||||
const uint32_t lower_bits = guf_rand32_u32(state);
|
||||
const uint32_t upper_bits = guf_rand32_u32(state);
|
||||
#ifdef UINT64_MAX
|
||||
return ((uint64_t)upper_bits << 32) | (uint64_t)lower_bits; // TODO: not sure if that's a good idea...
|
||||
#else
|
||||
return ((uint_least64_t)upper_bits << 32) | (uint_least64_t)lower_bits; // TODO: not sure if that's a good idea...
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
#ifdef UINT64_MAX
|
||||
GUF_RAND_KWRDS uint64_t guf_rand_u64(guf_randstate *state)
|
||||
#else
|
||||
GUF_RAND_KWRDS uint_least64_t guf_rand_u64(guf_randstate *state)
|
||||
#endif
|
||||
{
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
return guf_rand32_u64(state);
|
||||
#else
|
||||
return guf_rand64_u64(state);
|
||||
#endif
|
||||
}
|
||||
|
||||
/*
|
||||
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)
|
||||
Equivalent to 2^64 calls to guf_rand32_u32; can be used to generate 2^64
|
||||
non-overlapping subsequences for parallel computations.
|
||||
*/
|
||||
GUF_RAND_KWRDS void guf_randstate_jump(guf_randstate *state)
|
||||
GUF_RAND_KWRDS void guf_rand32_state_jump(guf_rand32_state *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) {
|
||||
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);
|
||||
guf_rand32_u32(state);
|
||||
}
|
||||
}
|
||||
state->s[0] = s0;
|
||||
state->s[1] = s1;
|
||||
state->s[2] = s2;
|
||||
state->s[3] = s3;
|
||||
#else
|
||||
}
|
||||
|
||||
|
||||
|
||||
#ifdef UINT64_MAX
|
||||
/*
|
||||
Equivalent to 2^128 calls to guf_rand64_u64(); can be used to generate 2^128
|
||||
non-overlapping subsequences for parallel computations.
|
||||
*/
|
||||
GUF_RAND_KWRDS void guf_rand64_state_jump(guf_rand64_state *state)
|
||||
{
|
||||
static const uint64_t JUMP[] = { 0x180ec6d33cfd0aba, 0xd5a61266f0c9392c, 0xa9582618e03fc9aa, 0x39abdc4529b1661c };
|
||||
uint64_t s0 = 0;
|
||||
uint64_t s1 = 0;
|
||||
@ -213,16 +348,44 @@ GUF_RAND_KWRDS void guf_randstate_jump(guf_randstate *state)
|
||||
s2 ^= state->s[2];
|
||||
s3 ^= state->s[3];
|
||||
}
|
||||
guf_rand_u64(state);
|
||||
guf_rand64_u64(state);
|
||||
}
|
||||
}
|
||||
state->s[0] = s0;
|
||||
state->s[1] = s1;
|
||||
state->s[2] = s2;
|
||||
state->s[3] = s3;
|
||||
}
|
||||
#endif
|
||||
|
||||
/*
|
||||
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.
|
||||
*/
|
||||
GUF_RAND_KWRDS void guf_randstate_jump(guf_randstate *state)
|
||||
{
|
||||
GUF_ASSERT(state);
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
guf_rand32_state_jump(state);
|
||||
#else
|
||||
guf_rand64_state_jump(state);
|
||||
#endif
|
||||
}
|
||||
|
||||
#ifdef UINT64_MAX
|
||||
GUF_RAND_KWRDS double guf_rand64_f64(guf_rand64_state *state)
|
||||
{
|
||||
// cf. https://prng.di.unimi.it/ and https://dotat.at/@/2023-06-23-random-double.html (last-retrieved 2025-02-11)
|
||||
return (guf_rand64_u64(state) >> 11) * 0x1.0p-53; // 11 == 64 - 53 (double has a 53-bit mantissa/significand)
|
||||
}
|
||||
#endif
|
||||
|
||||
GUF_RAND_KWRDS double guf_rand32_f64(guf_rand32_state *state)
|
||||
{
|
||||
// cf. https://prng.di.unimi.it/ and https://dotat.at/@/2023-06-23-random-double.html (last-retrieved 2025-02-11)
|
||||
return (guf_rand32_u64(state) >> 11) * 0x1.0p-53; // 11 == 64 - 53 (double has a 53-bit mantissa/significand)
|
||||
}
|
||||
|
||||
// Generate double in the unit interval [0, 1)
|
||||
GUF_RAND_KWRDS double guf_rand_f64(guf_randstate *state)
|
||||
{
|
||||
@ -230,28 +393,80 @@ GUF_RAND_KWRDS double guf_rand_f64(guf_randstate *state)
|
||||
return (guf_rand_u64(state) >> 11) * 0x1.0p-53; // 11 == 64 - 53 (double has a 53-bit mantissa/significand)
|
||||
}
|
||||
|
||||
#ifdef UINT64_MAX
|
||||
GUF_RAND_KWRDS float guf_rand64_f32(guf_rand64_state *state)
|
||||
{
|
||||
return (guf_rand64_u64(state) >> 40) * 0x1.0p-24f; // 40 == 64 - 24; (float has a 24-bit mantissa/significand)
|
||||
}
|
||||
#endif
|
||||
|
||||
GUF_RAND_KWRDS float guf_rand32_f32(guf_rand32_state *state)
|
||||
{
|
||||
return (guf_rand32_u32(state) >> 8) * 0x1.0p-24f; // 8 == 32 - 24; (float has a 24-bit mantissa/significand)
|
||||
}
|
||||
|
||||
// Generate float in the unit interval [0, 1)
|
||||
GUF_RAND_KWRDS 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)
|
||||
return guf_rand32_f32(state);
|
||||
#else
|
||||
return (guf_rand_u64(state) >> 40) * 0x1.0p-24f; // 40 == 64 - 24; (float has a 24-bit mantissa/significand)
|
||||
return guf_rand64_f32(state);
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
GUF_RAND_KWRDS bool guf_rand32_bernoulli_trial_f32(guf_rand32_state *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_rand_f64 is in range [0, 1)
|
||||
}
|
||||
#ifdef UINT64_MAX
|
||||
GUF_RAND_KWRDS bool guf_rand64_bernoulli_trial_f32(guf_rand64_state *state, float p)
|
||||
{
|
||||
p = guf_clamp_f32(p, 0, 1);
|
||||
return guf_rand64_f32(state) < p; // never true for p = 0, always true for p = 1 since guf_rand_f64 is in range [0, 1)
|
||||
}
|
||||
#endif
|
||||
|
||||
GUF_RAND_KWRDS 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_RAND_KWRDS bool guf_rand32_bernoulli_trial_f64(guf_rand32_state *state, double p)
|
||||
{
|
||||
p = guf_clamp_f64(p, 0, 1);
|
||||
return guf_rand32_f64(state) < p; // never true for p = 0, always true for p = 1 since guf_rand_f64 is in range [0, 1)
|
||||
}
|
||||
|
||||
#ifdef UINT64_MAX
|
||||
GUF_RAND_KWRDS bool guf_rand64_bernoulli_trial_f64(guf_rand64_state *state, double p)
|
||||
{
|
||||
p = guf_clamp_f64(p, 0, 1);
|
||||
return guf_rand64_f64(state) < p; // never true for p = 0, always true for p = 1 since guf_rand_f64 is in range [0, 1)
|
||||
}
|
||||
#endif
|
||||
|
||||
GUF_RAND_KWRDS 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)
|
||||
}
|
||||
|
||||
GUF_RAND_KWRDS bool guf_rand32_flip(guf_rand32_state *state)
|
||||
{
|
||||
return guf_rand32_bernoulli_trial_f32(state, 0.5f);
|
||||
}
|
||||
|
||||
#ifdef UINT64_MAX
|
||||
GUF_RAND_KWRDS bool guf_rand64_flip(guf_rand64_state *state)
|
||||
{
|
||||
return guf_rand64_bernoulli_trial_f64(state, 0.5);
|
||||
}
|
||||
#endif
|
||||
|
||||
GUF_RAND_KWRDS bool guf_rand_flip(guf_randstate *state)
|
||||
{
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
@ -261,8 +476,7 @@ GUF_RAND_KWRDS bool guf_rand_flip(guf_randstate *state)
|
||||
#endif
|
||||
}
|
||||
|
||||
// returns uniformly-distributed random double in range [min, end) (or min if min == end)
|
||||
GUF_RAND_KWRDS double guf_randrange_f64(guf_randstate *state, double min, double end)
|
||||
GUF_RAND_KWRDS double guf_rand32_range_f64(guf_rand32_state *state, double min, double end)
|
||||
{
|
||||
if (min == (double)INFINITY) {
|
||||
min = DBL_MAX;
|
||||
@ -275,11 +489,39 @@ GUF_RAND_KWRDS double guf_randrange_f64(guf_randstate *state, double min, double
|
||||
end = -DBL_MAX;
|
||||
}
|
||||
GUF_ASSERT_RELEASE(end >= min);
|
||||
return guf_rand_f64(state) * (end - min) + min;
|
||||
return guf_rand32_f64(state) * (end - min) + min;
|
||||
}
|
||||
|
||||
#ifdef UINT64_MAX
|
||||
GUF_RAND_KWRDS double guf_rand64_range_f64(guf_rand64_state *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;
|
||||
} else if (end == (double)-INFINITY) {
|
||||
end = -DBL_MAX;
|
||||
}
|
||||
GUF_ASSERT_RELEASE(end >= min);
|
||||
return guf_rand64_f64(state) * (end - min) + min;
|
||||
}
|
||||
#endif
|
||||
|
||||
// returns uniformly-distributed random double in range [min, end) (or min if min == end)
|
||||
GUF_RAND_KWRDS double guf_randrange_f64(guf_randstate *state, double min, double end)
|
||||
{
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
return guf_rand32_range_f64(state, min, end);
|
||||
#else
|
||||
return guf_rand64_range_f64(state, min, end);
|
||||
#endif
|
||||
}
|
||||
|
||||
// returns uniformly-distributed random float in range [min, end) (or min if min == end)
|
||||
GUF_RAND_KWRDS float guf_randrange_f32(guf_randstate *state, float min, float end)
|
||||
GUF_RAND_KWRDS float guf_rand32_range_f32(guf_rand32_state *state, float min, float end)
|
||||
{
|
||||
if (min == INFINITY) {
|
||||
min = FLT_MAX;
|
||||
@ -292,38 +534,145 @@ GUF_RAND_KWRDS float guf_randrange_f32(guf_randstate *state, float min, float en
|
||||
end = -FLT_MAX;
|
||||
}
|
||||
GUF_ASSERT_RELEASE(end >= min);
|
||||
return guf_rand_f32(state) * (end - min) + 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_RAND_KWRDS int32_t guf_randrange_i32(guf_randstate *state, int32_t min, int32_t max)
|
||||
#ifdef UINT64_MAX
|
||||
// returns uniformly-distributed random float in range [min, end) (or min if min == end)
|
||||
GUF_RAND_KWRDS float guf_rand64_range_f32(guf_rand64_state *state, float min, float end)
|
||||
{
|
||||
if (min == INFINITY) {
|
||||
min = FLT_MAX;
|
||||
} else if (min == -INFINITY) {
|
||||
min = -FLT_MAX;
|
||||
}
|
||||
if (end == INFINITY) {
|
||||
end = FLT_MAX;
|
||||
} else if (end == -INFINITY) {
|
||||
end = -FLT_MAX;
|
||||
}
|
||||
GUF_ASSERT_RELEASE(end >= min);
|
||||
return guf_rand64_f32(state) * (end - min) + min;
|
||||
}
|
||||
#endif
|
||||
|
||||
// returns uniformly-distributed random float in range [min, end) (or min if min == end)
|
||||
GUF_RAND_KWRDS float guf_randrange_f32(guf_randstate *state, float min, float end)
|
||||
{
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
return guf_rand32_range_f32(state, min, end);
|
||||
#else
|
||||
return guf_rand64_range_f32(state, min, end);
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
#ifdef UINT64_MAX
|
||||
// returns uniformly-distributed random int32_t in range [min, max] (max is inclusive as opposed to the f32/f64 versions)
|
||||
GUF_RAND_KWRDS int32_t guf_rand64_range_i32(guf_rand64_state *state, int32_t min, int32_t max)
|
||||
{
|
||||
GUF_ASSERT_RELEASE(max >= min);
|
||||
if (min == max) {
|
||||
return min;
|
||||
}
|
||||
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);
|
||||
const double result = floor(guf_rand64_f64(state) * (delta + 1.0) + min);
|
||||
GUF_ASSERT(result >= min && result <= max);
|
||||
return (int32_t)result;
|
||||
}
|
||||
#endif
|
||||
|
||||
// returns uniformly-distributed random int32_t in range [min, max] (max is inclusive as opposed to the f32/f64 versions)
|
||||
GUF_RAND_KWRDS int32_t guf_rand32_range_i32(guf_rand32_state *state, int32_t min, int32_t max)
|
||||
{
|
||||
GUF_ASSERT_RELEASE(max >= min);
|
||||
if (min == max) {
|
||||
return min;
|
||||
}
|
||||
|
||||
const uint32_t rand_max_i32 = UINT32_MAX >> 1; // 2^31 - 1 (== INT32_MAX)
|
||||
|
||||
const uint32_t delta = guf_absdiff_i32(max, min);
|
||||
if (delta > rand_max_i32) {
|
||||
guf_panic(GUF_ERR_INT_OVERFLOW, GUF_ERR_MSG("in function guf_rand32_range_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_i32 + 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);
|
||||
/*
|
||||
since (num_rand_vals % num_bins) is at most 2^30 + 1 (I think...), the minimum limit is 2^31 - (2^30 + 1),
|
||||
which means in the worst case, the chance of having to iterate (i.e. step >= limit)
|
||||
is 1 - (2^31 - (2^30 + 1)) / 2^31 == 0.5
|
||||
*/
|
||||
uint32_t step;
|
||||
do {
|
||||
step = guf_rand32_u32(state) >> 1; // [0, 2^31 - 1]
|
||||
} while (step >= limit);
|
||||
step = step / bin_size;
|
||||
|
||||
GUF_ASSERT(!guf_add_is_overflow_i32(min, step));
|
||||
const int32_t rnd = min + (int32_t)step;
|
||||
GUF_ASSERT(rnd >= min && rnd <= max);
|
||||
return rnd;
|
||||
}
|
||||
|
||||
GUF_RAND_KWRDS uint32_t guf_randrange_u32(guf_randstate *state, uint32_t min, uint32_t max)
|
||||
// returns uniformly-distributed random int32_t in range [min, max] (max is inclusive as opposed to the f32/f64 versions)
|
||||
GUF_RAND_KWRDS int32_t guf_randrange_i32(guf_randstate *state, int32_t min, int32_t max)
|
||||
{
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
return guf_rand32_range_i32(state, min, max);
|
||||
#else
|
||||
return guf_rand64_range_i32(state, min, max);
|
||||
#endif
|
||||
}
|
||||
|
||||
GUF_RAND_KWRDS uint32_t guf_rand32_range_u32(guf_rand32_state *state, uint32_t min, uint32_t max)
|
||||
{
|
||||
GUF_ASSERT_RELEASE(max >= min);
|
||||
if (min == max) {
|
||||
return min;
|
||||
}
|
||||
const double delta = (double)max - (double)min;
|
||||
const double result = floor(guf_rand_f64(state) * (delta + 1.0) + min);
|
||||
const double result = floor(guf_rand32_f64(state) * (delta + 1.0) + min); // TODO: slow for 32-bit platforms...
|
||||
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_RAND_KWRDS int64_t guf_randrange_i64(guf_randstate *state, int64_t min, int64_t max)
|
||||
#ifdef UINT64_MAX
|
||||
GUF_RAND_KWRDS uint32_t guf_rand64_range_u32(guf_rand64_state *state, uint32_t min, uint32_t max)
|
||||
{
|
||||
GUF_ASSERT_RELEASE(max >= min);
|
||||
if (min == max) {
|
||||
return min;
|
||||
}
|
||||
const double delta = (double)max - (double)min;
|
||||
const double result = floor(guf_rand64_f64(state) * (delta + 1.0) + min);
|
||||
GUF_ASSERT(result >= min && result <= max);
|
||||
return (uint32_t)result;
|
||||
}
|
||||
#endif
|
||||
|
||||
GUF_RAND_KWRDS uint32_t guf_randrange_u32(guf_randstate *state, uint32_t min, uint32_t max)
|
||||
{
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
return guf_rand32_range_u32(state, min, max);
|
||||
#else
|
||||
return guf_rand64_range_u32(state, min, max);
|
||||
#endif
|
||||
}
|
||||
|
||||
#ifdef UINT64_MAX
|
||||
GUF_RAND_KWRDS int64_t guf_rand64_range_i64(guf_rand64_state *state, int64_t min, int64_t max)
|
||||
{
|
||||
GUF_ASSERT_RELEASE(max >= min);
|
||||
if (min == max) {
|
||||
return min;
|
||||
@ -355,19 +704,28 @@ GUF_RAND_KWRDS int64_t guf_randrange_i64(guf_randstate *state, int64_t min, int6
|
||||
*/
|
||||
uint64_t step;
|
||||
do {
|
||||
step = guf_rand_u64(state) >> 1; // [0, 2^63 - 1]
|
||||
step = guf_rand64_u64(state) >> 1; // [0, 2^63 - 1]
|
||||
} while (step >= limit);
|
||||
step = step / bin_size;
|
||||
|
||||
const int64_t rnd = min + step;
|
||||
GUF_ASSERT(!guf_add_is_overflow_i64(min, step));
|
||||
const int64_t rnd = min + (int64_t)step;
|
||||
GUF_ASSERT(rnd >= min && rnd <= max);
|
||||
return rnd;
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
#if !defined(GUF_RAND_32_BIT) && defined(UINT64_MAX)
|
||||
// returns uniformly-distributed random int64_t in range [min, max] (max is inclusive as opposed to the f32/f64 versions)
|
||||
GUF_RAND_KWRDS int64_t guf_randrange_i64(guf_randstate *state, int64_t min, int64_t max)
|
||||
{
|
||||
return guf_rand64_range_i64(state, min, max);
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
// Box-Müller-transform transcribed from wikipedia, cf. https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform (last-retrieved 2025-02-12)
|
||||
|
||||
GUF_RAND_KWRDS void guf_rand_normal_sample_f64(guf_randstate *state, double mean, double std_dev, double *result, ptrdiff_t n)
|
||||
GUF_RAND_KWRDS void guf_rand32_normal_sample_f64(guf_rand32_state *state, double mean, double std_dev, double *result, ptrdiff_t n)
|
||||
{
|
||||
GUF_ASSERT_RELEASE(result);
|
||||
GUF_ASSERT_RELEASE(n >= 0);
|
||||
@ -377,9 +735,9 @@ GUF_RAND_KWRDS void guf_rand_normal_sample_f64(guf_randstate *state, double mean
|
||||
while (i < n) {
|
||||
double u1, u2;
|
||||
do {
|
||||
u1 = guf_rand_f64(state);
|
||||
u1 = guf_rand32_f64(state);
|
||||
} while (u1 == 0);
|
||||
u2 = guf_rand_f64(state);
|
||||
u2 = guf_rand32_f64(state);
|
||||
|
||||
const double mag = std_dev * sqrt(-2.0 * log(u1));
|
||||
result[i++] = mag * cos(TAU * u2) + mean;
|
||||
@ -389,7 +747,43 @@ GUF_RAND_KWRDS void guf_rand_normal_sample_f64(guf_randstate *state, double mean
|
||||
}
|
||||
}
|
||||
|
||||
GUF_RAND_KWRDS void guf_rand_normal_sample_f32(guf_randstate *state, float mean, float std_dev, float *result, ptrdiff_t n)
|
||||
#ifdef UINT64_MAX
|
||||
// Box-Müller-transform transcribed from wikipedia, cf. https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform (last-retrieved 2025-02-12)
|
||||
GUF_RAND_KWRDS void guf_rand64_normal_sample_f64(guf_rand64_state *state, double mean, double std_dev, double *result, ptrdiff_t n)
|
||||
{
|
||||
GUF_ASSERT_RELEASE(result);
|
||||
GUF_ASSERT_RELEASE(n >= 0);
|
||||
const double TAU = 2.0 * GUF_PI;
|
||||
|
||||
ptrdiff_t i = 0;
|
||||
while (i < n) {
|
||||
double u1, u2;
|
||||
do {
|
||||
u1 = guf_rand64_f64(state);
|
||||
} while (u1 == 0);
|
||||
u2 = guf_rand64_f64(state);
|
||||
|
||||
const double mag = std_dev * sqrt(-2.0 * log(u1));
|
||||
result[i++] = mag * cos(TAU * u2) + mean;
|
||||
if (i < n) {
|
||||
result[i++] = mag * sin(TAU * u2) + mean;
|
||||
}
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
// Box-Müller-transform transcribed from wikipedia, cf. https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform (last-retrieved 2025-02-12)
|
||||
GUF_RAND_KWRDS void guf_rand_normal_sample_f64(guf_randstate *state, double mean, double std_dev, double *result, ptrdiff_t n)
|
||||
{
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
guf_rand32_normal_sample_f64(state, mean, std_dev, result, n);
|
||||
#else
|
||||
guf_rand64_normal_sample_f64(state, mean, std_dev, result, n);
|
||||
#endif
|
||||
}
|
||||
|
||||
GUF_RAND_KWRDS void guf_rand32_normal_sample_f32(guf_rand32_state *state, float mean, float std_dev, float *result, ptrdiff_t n)
|
||||
{
|
||||
GUF_ASSERT_RELEASE(result);
|
||||
GUF_ASSERT_RELEASE(n >= 0);
|
||||
@ -399,9 +793,9 @@ GUF_RAND_KWRDS void guf_rand_normal_sample_f32(guf_randstate *state, float mean,
|
||||
while (i < n) {
|
||||
float u1, u2;
|
||||
do {
|
||||
u1 = guf_rand_f32(state);
|
||||
u1 = guf_rand32_f32(state);
|
||||
} while (u1 == 0);
|
||||
u2 = guf_rand_f32(state);
|
||||
u2 = guf_rand32_f32(state);
|
||||
|
||||
const float mag = std_dev * sqrtf(-2.f * logf(u1));
|
||||
result[i++] = mag * cosf(TAU * u2) + mean;
|
||||
@ -411,23 +805,91 @@ GUF_RAND_KWRDS void guf_rand_normal_sample_f32(guf_randstate *state, float mean,
|
||||
}
|
||||
}
|
||||
|
||||
GUF_RAND_KWRDS double guf_rand_normal_sample_one_f64(guf_randstate *state, double mean, double std_dev)
|
||||
#ifdef UINT64_MAX
|
||||
GUF_RAND_KWRDS void guf_rand64_normal_sample_f32(guf_rand64_state *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_rand64_f32(state);
|
||||
} while (u1 == 0);
|
||||
u2 = guf_rand64_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;
|
||||
}
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
GUF_RAND_KWRDS void guf_rand_normal_sample_f32(guf_randstate *state, float mean, float std_dev, float *result, ptrdiff_t n)
|
||||
{
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
guf_rand32_normal_sample_f32(state, mean, std_dev, result, n);
|
||||
#else
|
||||
guf_rand64_normal_sample_f32(state, mean, std_dev, result, n);
|
||||
#endif
|
||||
}
|
||||
|
||||
GUF_RAND_KWRDS double guf_rand32_normal_sample_one_f64(guf_rand32_state *state, double mean, double std_dev)
|
||||
{
|
||||
double result;
|
||||
guf_rand_normal_sample_f64(state, mean, std_dev, &result, 1);
|
||||
guf_rand32_normal_sample_f64(state, mean, std_dev, &result, 1);
|
||||
return result;
|
||||
}
|
||||
|
||||
GUF_RAND_KWRDS float guf_rand_normal_sample_one_f32(guf_randstate *state, float mean, float std_dev)
|
||||
#ifdef UINT64_MAX
|
||||
GUF_RAND_KWRDS double guf_rand64_normal_sample_one_f64(guf_rand64_state *state, double mean, double std_dev)
|
||||
{
|
||||
double result;
|
||||
guf_rand64_normal_sample_f64(state, mean, std_dev, &result, 1);
|
||||
return result;
|
||||
}
|
||||
#endif
|
||||
|
||||
GUF_RAND_KWRDS double guf_rand_normal_sample_one_f64(guf_randstate *state, double mean, double std_dev)
|
||||
{
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
return guf_rand32_normal_sample_one_f64(state, mean, std_dev);
|
||||
#else
|
||||
return guf_rand64_normal_sample_one_f64(state, mean, std_dev);
|
||||
#endif
|
||||
}
|
||||
|
||||
GUF_RAND_KWRDS float guf_rand32_normal_sample_one_f32(guf_rand32_state *state, float mean, float std_dev)
|
||||
{
|
||||
float result;
|
||||
guf_rand_normal_sample_f32(state, mean, std_dev, &result, 1);
|
||||
guf_rand32_normal_sample_f32(state, mean, std_dev, &result, 1);
|
||||
return result;
|
||||
}
|
||||
|
||||
#ifdef UINT64_MAX
|
||||
GUF_RAND_KWRDS float guf_rand64_normal_sample_one_f32(guf_rand64_state *state, float mean, float std_dev)
|
||||
{
|
||||
float result;
|
||||
guf_rand64_normal_sample_f32(state, mean, std_dev, &result, 1);
|
||||
return result;
|
||||
}
|
||||
#endif
|
||||
|
||||
GUF_RAND_KWRDS float guf_rand_normal_sample_one_f32(guf_randstate *state, float mean, float std_dev)
|
||||
{
|
||||
#ifdef GUF_RAND_32_BIT
|
||||
return guf_rand32_normal_sample_one_f32(state, mean, std_dev);
|
||||
#else
|
||||
return guf_rand64_normal_sample_one_f32(state, mean, std_dev);
|
||||
#endif
|
||||
}
|
||||
|
||||
#undef GUF_RAND_IMPL
|
||||
#undef GUF_RAND_IMPL_STATIC
|
||||
#endif /* end impl */
|
||||
|
||||
#undef GUF_RAND_KWRDS
|
||||
#undef GUF_RAND_32_BIT
|
||||
|
||||
@ -4,19 +4,69 @@
|
||||
|
||||
#ifndef GUF_UTILS_H
|
||||
#define GUF_UTILS_H
|
||||
#include "guf_common.h"
|
||||
#include "guf_assert.h"
|
||||
|
||||
static inline bool guf_platform_is_big_endian(void)
|
||||
static inline void guf_platform_assert_endianness(void)
|
||||
{
|
||||
unsigned i = 1;
|
||||
const unsigned i = 1;
|
||||
const char *bytes = (const char*)&i;
|
||||
return bytes[0] != 1;
|
||||
const bool is_big_endian = bytes[0] != 1;
|
||||
#if defined(GUF_PLATFORM_LITTLE_ENDIAN)
|
||||
GUF_ASSERT_RELEASE(!is_big_endian)
|
||||
#elif defined(GUF_PLATFORM_BIG_ENDIAN)
|
||||
GUF_ASSERT_RELEASE(is_big_endian)
|
||||
#endif
|
||||
}
|
||||
|
||||
static inline int guf_platform_native_word_bits(void)
|
||||
static inline void guf_platform_assert_native_word_bits(void)
|
||||
{
|
||||
const int bits = sizeof(void*) * CHAR_BIT;
|
||||
return bits;
|
||||
GUF_ASSERT_RELEASE(GUF_PLATFORM_BITS == bits);
|
||||
}
|
||||
|
||||
#ifdef NDEBUG
|
||||
#define GUF_DBG_STR "release"
|
||||
#else
|
||||
#define GUF_DBG_STR "debug"
|
||||
#endif
|
||||
|
||||
#if defined(GUF_STDC_AT_LEAST_C23)
|
||||
#ifdef GUF_PLATFORM_LITTLE_ENDIAN
|
||||
#define GUF_PLATFORM_STRING "C23 (or above) " GUF_STRINGIFY(GUF_PLATFORM_BITS) "-bit " "little-endian " GUF_DBG_STR
|
||||
#else
|
||||
#define GUF_PLATFORM_STRING "C23 (or above) " GUF_STRINGIFY(GUF_PLATFORM_BITS) "-bit " "big-endian " GUF_DBG_STR
|
||||
#endif
|
||||
#elif defined(GUF_STDC_AT_LEAST_C17)
|
||||
#ifdef GUF_PLATFORM_LITTLE_ENDIAN
|
||||
#define GUF_PLATFORM_STRING "C17 " GUF_STRINGIFY(GUF_PLATFORM_BITS) "-bit " "little-endian " GUF_DBG_STR
|
||||
#else
|
||||
#define GUF_PLATFORM_STRING "C17 " GUF_STRINGIFY(GUF_PLATFORM_BITS) "-bit " "big-endian " GUF_DBG_STR
|
||||
#endif
|
||||
#elif defined(GUF_STDC_AT_LEAST_C11)
|
||||
#ifdef GUF_PLATFORM_LITTLE_ENDIAN
|
||||
#define GUF_PLATFORM_STRING "C11 " GUF_STRINGIFY(GUF_PLATFORM_BITS) "-bit " "little-endian " GUF_DBG_STR
|
||||
#else
|
||||
#define GUF_PLATFORM_STRING "C11 " GUF_STRINGIFY(GUF_PLATFORM_BITS) "-bit " "big-endian " GUF_DBG_STR
|
||||
#endif
|
||||
#elif defined(GUF_STDC_AT_LEAST_C99)
|
||||
#ifdef GUF_PLATFORM_LITTLE_ENDIAN
|
||||
#define GUF_PLATFORM_STRING "C99 " GUF_STRINGIFY(GUF_PLATFORM_BITS) "-bit " "little-endian " GUF_DBG_STR
|
||||
#else
|
||||
#define GUF_PLATFORM_STRING "C99 " GUF_STRINGIFY(GUF_PLATFORM_BITS) "-bit " "big-endian " GUF_DBG_STR
|
||||
#endif
|
||||
#elif defined(GUF_STDCPP_AT_LEAST_CPP11)
|
||||
#ifdef GUF_PLATFORM_LITTLE_ENDIAN
|
||||
#define GUF_PLATFORM_STRING "C++11 (or above) " GUF_STRINGIFY(GUF_PLATFORM_BITS) "-bit " "little-endian " GUF_DBG_STR
|
||||
#else
|
||||
#define GUF_PLATFORM_STRING "C++11 (or above) " GUF_STRINGIFY(GUF_PLATFORM_BITS) "-bit " "big-endian " GUF_DBG_STR
|
||||
#endif
|
||||
#else
|
||||
#ifdef GUF_PLATFORM_LITTLE_ENDIAN
|
||||
#define GUF_PLATFORM_STRING "C?? " GUF_STRINGIFY(GUF_PLATFORM_BITS) "-bit " "little-endian " GUF_DBG_STR
|
||||
#else
|
||||
#define GUF_PLATFORM_STRING "C?? " GUF_STRINGIFY(GUF_PLATFORM_BITS) "-bit " "big-endian " GUF_DBG_STR
|
||||
#endif
|
||||
#endif
|
||||
|
||||
#endif
|
||||
|
||||
@ -8,6 +8,7 @@
|
||||
#include "guf_alloc_libc.h"
|
||||
#include "guf_cstr.h"
|
||||
#include "guf_linalg.h"
|
||||
#include "guf_utils.h"
|
||||
|
||||
#define GUF_T float
|
||||
#define GUF_SORT_IMPL_STATIC
|
||||
@ -46,13 +47,16 @@
|
||||
#include "guf_dbuf.h"
|
||||
|
||||
#define GUF_RAND_IMPL_STATIC
|
||||
// #define GUF_RAND_32_BIT
|
||||
#include "guf_rand.h"
|
||||
|
||||
#include "impls/dict_impl.h"
|
||||
|
||||
int main(void)
|
||||
{
|
||||
printf("libguf test: compiled with C %ld\n", __STDC_VERSION__);
|
||||
printf("libguf example: " GUF_PLATFORM_STRING "\n");
|
||||
guf_platform_assert_endianness();
|
||||
guf_platform_assert_native_word_bits();
|
||||
|
||||
guf_allocator test_allocator = guf_allocator_libc;
|
||||
guf_libc_alloc_ctx test_allocator_ctx = {.alloc_type_id = 0, .thread_id = 0, .zero_init = true};
|
||||
@ -206,7 +210,7 @@ int main(void)
|
||||
dbuf_int_free(&integers, NULL);
|
||||
printf("\n");
|
||||
guf_randstate rng;
|
||||
guf_randstate_init(&rng, time(NULL));
|
||||
guf_randstate_init(&rng, (guf_rand_seed_t)time(NULL));
|
||||
int heads = 0, tails = 0;
|
||||
int throws = 10;
|
||||
for (i = 0; i < throws; ++i) {
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user