From 6b222eafc154c2e5d8745236c17f7fc75c5f16d3 Mon Sep 17 00:00:00 2001 From: jun <83899451+zeichensystem@users.noreply.github.com> Date: Sat, 15 Feb 2025 03:00:15 +0100 Subject: [PATCH] Add guf_rand --- src/guf_dbuf.h | 15 +-- src/guf_int.h | 37 +++--- src/guf_math.h | 61 ++++++++++ src/guf_rand.h | 321 +++++++++++++++++++++++++++++++++++++++++++++++++ src/guf_sort.h | 2 +- src/guf_test.c | 46 +++++++ 6 files changed, 454 insertions(+), 28 deletions(-) create mode 100644 src/guf_math.h create mode 100644 src/guf_rand.h diff --git a/src/guf_dbuf.h b/src/guf_dbuf.h index 7ac729b..d9fe2a8 100644 --- a/src/guf_dbuf.h +++ b/src/guf_dbuf.h @@ -59,7 +59,7 @@ typedef struct GUF_CNT_NAME { /* - Regular iterator: base is always NULL - - Reverse iterator: base points the element after ptr + - Reverse iterator: base points the element after ptr (unless dbuf->size == 0, then ptr and base are NULL for both iterator types) Examples: - rbegin(): base points to end().ptr, and ptr to end().ptr - 1 - rend(): base points to begin().ptr and ptr to NULL @@ -68,7 +68,7 @@ typedef struct GUF_CNT_NAME { typedef struct GUF_CAT(GUF_CNT_NAME, _iter) { GUF_T *ptr; - GUF_T *base; + GUF_T *base; // Not NULL For reverse iterators (unless dbuf->size == 0, then ptr and base are NULL for both iterator types) } GUF_CAT(GUF_CNT_NAME, _iter); GUF_FN_KEYWORDS bool GUF_CAT(GUF_CNT_NAME, _valid)(const GUF_CNT_NAME* dbuf); @@ -769,11 +769,8 @@ GUF_FN_KEYWORDS void GUF_CAT(GUF_CNT_NAME, _try_shrink_to_fit)(GUF_CNT_NAME *dbu guf_err_set_if_not_null(err, GUF_ERR_NONE); } - /* Iterator functions */ -#define ITER_PTR(it) !it.reverse ? (it.iter.ptr) : (it.rev_iter.ptr - 1) - GUF_FN_KEYWORDS GUF_CAT(GUF_CNT_NAME, _iter) GUF_CAT(GUF_CNT_NAME, _begin)(const GUF_CNT_NAME* dbuf) { GUF_ASSERT_RELEASE(GUF_CAT(GUF_CNT_NAME, _valid)(dbuf)); @@ -822,9 +819,9 @@ GUF_FN_KEYWORDS GUF_CAT(GUF_CNT_NAME, _iter) GUF_CAT(GUF_CNT_NAME, _iter_at_idx return it; } - if (idx < 0) { + if (idx <= 0) { // begin() it.ptr = dbuf->data; - } else if (idx > dbuf->size) { + } else if (idx >= dbuf->size) { // end() it.ptr = dbuf->data + dbuf->size; } else { it.ptr = dbuf->data + idx; @@ -843,10 +840,10 @@ GUF_FN_KEYWORDS GUF_CAT(GUF_CNT_NAME, _iter) GUF_CAT(GUF_CNT_NAME, _reverse_ite return it; } - if (idx < 0) { + if (idx < 0) { // rend() it.base = dbuf->data; it.ptr = NULL; - } else if (idx >= dbuf->size) { + } else if (idx >= dbuf->size) { // rbegin() it.base = dbuf->data + dbuf->size; it.ptr = it.base - 1; } else { diff --git a/src/guf_int.h b/src/guf_int.h index 85d3efd..946e6cc 100644 --- a/src/guf_int.h +++ b/src/guf_int.h @@ -3,42 +3,43 @@ #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) {return GUF_CAT(guf_max_, int_type_name)(GUF_CAT(guf_min_, int_type_name)(x, max), min);} + 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(short, short) -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) -GUF_DEFINE_MIN_MAX_CLAMP(unsigned char, unsigned_char) +GUF_DEFINE_MIN_MAX_CLAMP(unsigned char, uchar) GUF_DEFINE_MIN_MAX_CLAMP(unsigned, unsigned) -GUF_DEFINE_MIN_MAX_CLAMP(unsigned short, u_short) -GUF_DEFINE_MIN_MAX_CLAMP(unsigned long, unsigned_long) -GUF_DEFINE_MIN_MAX_CLAMP(unsigned long long, unsigned_long_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) +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 - -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;} - #endif diff --git a/src/guf_math.h b/src/guf_math.h new file mode 100644 index 0000000..49dbfb3 --- /dev/null +++ b/src/guf_math.h @@ -0,0 +1,61 @@ +#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 float guf_clamp_f32(float x, float min, float max) +{ + if (x < min) return min; + if (x > max) return max; + return x; +} + +static inline double guf_clamp_f64(double x, double min, double max) +{ + if (x < min) return min; + if (x > max) return max; + return x; +} + +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; +} + +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); +} + +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;} + +#endif diff --git a/src/guf_rand.h b/src/guf_rand.h new file mode 100644 index 0000000..ed2e54a --- /dev/null +++ b/src/guf_rand.h @@ -0,0 +1,321 @@ +#ifndef GUF_RAND_H +#define GUF_RAND_H +#include "guf_common.h" +#include "guf_assert.h" +#include "guf_math.h" +#include +#include + +#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 +#else + #define GUF_FN_KEYWORDS +#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] +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_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] + +// normal distributions +GUF_FN_KEYWORDS void guf_rand_normal_sample_f64(guf_randstate *state, double mean, double std_dev, double *result, ptrdiff_t n); +GUF_FN_KEYWORDS void guf_rand_normal_sample_f32(guf_randstate *state, float mean, float std_dev, float *result, ptrdiff_t n); +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); + +#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_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); +} + +GUF_FN_KEYWORDS void guf_randstate_init(guf_randstate *state, uint64_t seed) +{ + GUF_ASSERT_RELEASE(state); + + 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]; + for (size_t i = 1; i < GUF_STATIC_BUF_SIZE(state->s); ++i) { + state->s[i] = guf_rand_splitmix64(&seed); + } + } +} + +/* + 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 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]); + + 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; +} + +/* + 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_randstate_jump(guf_randstate *state) +{ + GUF_ASSERT(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_rand_u64(state); + } + } + state->s[0] = s0; + state->s[1] = s1; + state->s[2] = s2; + state->s[3] = s3; +} + +// Generate double in the unit interval [0, 1) +GUF_FN_KEYWORDS double guf_rand_f64(guf_randstate *state) +{ + // cf. https://prng.di.unimi.it/ and https://dotat.at/@/2023-06-23-random-double.html (last-retrieved 2025-02-11) + return (guf_rand_u64(state) >> 11) * 0x1.0p-53; // 11 == 64 - 53 (double has a 53-bit mantissa/significand) +} + +// Generate float in the unit interval [0, 1) +GUF_FN_KEYWORDS float guf_rand_f32(guf_randstate *state) +{ + return (guf_rand_u64(state) >> 40) * 0x1.0p-24f; // 40 == 64 - 24; (float has a 24-bit mantissa/significand) +} + +GUF_FN_KEYWORDS bool guf_rand_bernoulli_trial(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_FN_KEYWORDS bool guf_rand_flip(guf_randstate *state) +{ + return guf_rand_bernoulli_trial(state, 0.5); +} + +// 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 (end == (double)INFINITY) { + end = DBL_MAX; + } + if (min == (double)-INFINITY) { + min = -DBL_MAX; + } + GUF_ASSERT_RELEASE(end >= min); + return guf_rand_f64(state) * (end - min) + min; +} + +// 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 (end == INFINITY) { + end = FLT_MAX; + } + if (min == -INFINITY) { + min = -FLT_MAX; + } + GUF_ASSERT_RELEASE(end >= min); + return guf_rand_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_randrange_i32(guf_randstate *state, int32_t min, int32_t max) +{ + GUF_ASSERT_RELEASE(max >= min); + if (min == max) { + return min; + } + const double delta = (int64_t)max - (int64_t)min; // Cast to int64_t to avoid overflow. + + // 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; +} + +GUF_FN_KEYWORDS 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 double delta = max - min; // Cannot overflow here. + + 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) +{ + GUF_ASSERT_RELEASE(max >= min); + if (min == max) { + 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 delta = max - min; + if (delta > rand_max) { + 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 + 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) >> rand_max_shift; + } while (step >= limit); + step = step / bin_size; + + const int64_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_rand_normal_sample_f64(guf_randstate *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_rand_f64(state); + } while (u1 == 0); + u2 = guf_rand_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; + } + } +} + +GUF_FN_KEYWORDS void guf_rand_normal_sample_f32(guf_randstate *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_rand_f32(state); + } while (u1 == 0); + u2 = guf_rand_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 double guf_rand_normal_sample_one_f64(guf_randstate *state, double mean, double std_dev) +{ + double result; + guf_rand_normal_sample_f64(state, mean, std_dev, &result, 1); + return result; +} + +GUF_FN_KEYWORDS float guf_rand_normal_sample_one_f32(guf_randstate *state, float mean, float std_dev) +{ + float result; + guf_rand_normal_sample_f32(state, mean, std_dev, &result, 1); + return result; +} + +#undef GUF_IMPL +#undef GUF_IMPL_STATIC +#endif + +#undef GUF_FN_KEYWORDS +#endif diff --git a/src/guf_sort.h b/src/guf_sort.h index 29db5ed..435ef7c 100644 --- a/src/guf_sort.h +++ b/src/guf_sort.h @@ -180,9 +180,9 @@ GUF_FN_KEYWORDS bool GUF_CAT(GUF_FN_NAME_PREFIX, _is_sorted)(GUF_T *arr, ptrdiff #undef GUF_IMPL #undef GUF_IMPL_STATIC -#undef GUF_FN_KEYWORDS #endif /* end #ifdef GUF_IMPL */ +#undef GUF_FN_KEYWORDS #undef GUF_T #undef GUF_FN_NAME_PREFIX #endif /* end #ifdef GUF_T */ diff --git a/src/guf_test.c b/src/guf_test.c index aca8bad..09651dd 100644 --- a/src/guf_test.c +++ b/src/guf_test.c @@ -1,6 +1,7 @@ #include #include #include +#include #include "guf_init.h" /* Must be included once (sets up the global panic handler) */ #include "guf_alloc_libc.h" @@ -10,6 +11,10 @@ #define GUF_T float #include "guf_sort.h" +#define GUF_IMPL_STATIC +#define GUF_T int +#include "guf_sort.h" + #define GUF_CNT_NAME dbuf_int #define GUF_T int #define GUF_T_IS_INTEGRAL_TYPE @@ -38,6 +43,9 @@ #define GUF_IMPL_STATIC #include "guf_dbuf.h" +#define GUF_IMPL_STATIC +#include "guf_rand.h" + int main(void) { printf("libguf test: compiled with C %ld\n", __STDC_VERSION__); @@ -154,5 +162,43 @@ int main(void) } dbuf_int_free(&integers, NULL); + + guf_randstate rng; + guf_randstate_init(&rng, time(NULL)); + guf_randstate_jump(&rng); + + printf("\n"); + int heads = 0, tails = 0; + int throws = 10; + for (i = 0; i < throws; ++i) { + bool is_head = guf_rand_flip(&rng); + if (is_head) { + puts("head"); + ++heads; + } else { + puts("tail"); + ++tails; + } + } + printf("n: %d\nheads: %d\ntails: %d\n", throws, heads, tails); + + 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 (size_t n = 50; n <= 150; ++n) { + printf("%zu:\t", n); + for (int j = 0; j < result[n] / 8; ++j) { + putc('#', stdout); + } + puts(""); + } + return EXIT_SUCCESS; }