diff --git a/src/guf_common.h b/src/guf_common.h index 86c5404..30c9c30 100644 --- a/src/guf_common.h +++ b/src/guf_common.h @@ -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) - #define GUF_STDC_AT_LEAST_C11 +#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) - #define GUF_STDCPP_AT_LEAST_CPP11 +#if defined(__cplusplus) + #if __cplusplus >= 201103L + #define GUF_STDCPP_AT_LEAST_CPP11 + #endif #endif /* diff --git a/src/guf_hash.h b/src/guf_hash.h index 777a936..1ab380a 100644 --- a/src/guf_hash.h +++ b/src/guf_hash.h @@ -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,18 +66,20 @@ 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) -{ - 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 = UINT64_C(1099511628211); - for (ptrdiff_t i = 0; i < num_bytes; ++i) { - hash ^= data_bytes[i]; - hash *= FNV_64_PRIME; +#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... + const uint64_t FNV_64_PRIME = UINT64_C(1099511628211); + for (ptrdiff_t i = 0; i < num_bytes; ++i) { + hash ^= data_bytes[i]; + hash *= FNV_64_PRIME; + } + return hash; } - return hash; -} +#endif #undef GUF_HASH_IMPL #undef GUF_HASH_IMPL_STATIC diff --git a/src/guf_math.h b/src/guf_math.h index 72bd996..66c1a9b 100644 --- a/src/guf_math.h +++ b/src/guf_math.h @@ -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) { diff --git a/src/guf_rand.h b/src/guf_rand.h index 15771b9..2988b06 100644 --- a/src/guf_rand.h +++ b/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 - uint32_t s[4]; - } guf_randstate; -#else - #define GUF_RAND_MAX UINT64_MAX - typedef struct guf_randstate { // State for xoshiro256** 1.0 +typedef struct guf_rand32_state { + uint32_t s[4]; +} 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) -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 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 float guf_randrange_f32(guf_randstate *state, float min, float 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" -/* - 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_ASSERT(state); - uint64_t z = ((*state) += 0x9e3779b97f4a7c15); - z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9; - z = (z ^ (z >> 27)) * 0x94d049bb133111eb; - return z ^ (z >> 31); -} +#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_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 ) @@ -91,42 +168,49 @@ 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; - for (size_t i = 1; i < GUF_ARR_SIZE(state->s); ++i) { - state->s[i] = (uint32_t)(guf_rand_splitmix64(&seed) >> 32); - } - } - #else - for (size_t i = 0; i < GUF_ARR_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]; for (size_t i = 1; i < GUF_ARR_SIZE(state->s); ++i) { - state->s[i] = guf_rand_splitmix64(&seed); + state->s[i] = guf_rand_splitmix32(&seed); } } - #endif } -GUF_RAND_KWRDS uint32_t guf_rand_u32(guf_randstate *state) +#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); + } + 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]; + for (size_t i = 1; i < GUF_ARR_SIZE(state->s); ++i) { + 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_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,37 +224,140 @@ 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 - /* - 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; + 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) + */ + 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 + +#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^64 calls to guf_rand32_u32; can be used to generate 2^64 + non-overlapping subsequences for parallel computations. +*/ +GUF_RAND_KWRDS void guf_rand32_state_jump(guf_rand32_state *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; +} + + + +#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; + 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]; + s1 ^= state->s[1]; + s2 ^= state->s[2]; + s3 ^= state->s[3]; + } + 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. @@ -179,50 +366,26 @@ GUF_RAND_KWRDS 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; + guf_rand32_state_jump(state); #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]; - s1 ^= state->s[1]; - s2 ^= state->s[2]; - s3 ^= state->s[3]; - } - guf_rand_u64(state); - } - } - state->s[0] = s0; - state->s[1] = s1; - state->s[2] = s2; - state->s[3] = s3; + 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,39 +393,90 @@ 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 - return guf_rand_bernoulli_trial_f32(state, 0.5f); + return guf_rand_bernoulli_trial_f32(state, 0.5f); #else - return guf_rand_bernoulli_trial_f64(state, 0.5); + 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_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,82 +534,198 @@ 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; } +#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_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_randrange_i32(guf_randstate *state, int32_t min, int32_t max) +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 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); - return (int32_t)result; + + 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) { - GUF_ASSERT_RELEASE(max >= min); - if (min == max) { - return min; - } - - 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_i64) { - guf_panic(GUF_ERR_INT_OVERFLOW, GUF_ERR_MSG("in function guf_randrange_i64: interval [min, max] larger than INT64_MAX")); - return -1; - } - /* - We should not use the same approach as in guf_randrange_i32 because (max - min) might be close to 2^63 - 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 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) - const uint64_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^62 + 1 (I think...), the minimum limit is 2^63 - (2^62 + 1), - which means in the worst case, the chance of having to iterate (i.e. step >= limit) - is 1 - (2^63 - (2^62 + 1)) / 2^63 == 0.5 - */ - uint64_t step; - do { - step = guf_rand_u64(state) >> 1; // [0, 2^63 - 1] - } while (step >= limit); - step = step / bin_size; - - const int64_t rnd = min + step; - GUF_ASSERT(rnd >= min && rnd <= max); - return rnd; + #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; + } + + 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_i64) { + guf_panic(GUF_ERR_INT_OVERFLOW, GUF_ERR_MSG("in function guf_randrange_i64: interval [min, max] larger than INT64_MAX")); + return -1; + } + /* + We should not use the same approach as in guf_randrange_i32 because (max - min) might be close to 2^63 - 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 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) + const uint64_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^62 + 1 (I think...), the minimum limit is 2^63 - (2^62 + 1), + which means in the worst case, the chance of having to iterate (i.e. step >= limit) + is 1 - (2^63 - (2^62 + 1)) / 2^63 == 0.5 + */ + uint64_t step; + do { + step = guf_rand64_u64(state) >> 1; // [0, 2^63 - 1] + } while (step >= limit); + step = step / bin_size; + + 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 diff --git a/src/guf_utils.h b/src/guf_utils.h index 8ec946c..9b85090 100644 --- a/src/guf_utils.h +++ b/src/guf_utils.h @@ -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 diff --git a/src/test/example.c b/src/test/example.c index 7e948f9..41aad68 100644 --- a/src/test/example.c +++ b/src/test/example.c @@ -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) {