fruit/libguf
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add guf_rand

Browse Source
This commit is contained in:
jun 2025-02-15 03:00:15 +01:00
parent fd58daa8b5
commit 6b222eafc1
6 changed files with 454 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
src/guf_dbuf.h
View File

@ -59,7 +59,7 @@ typedef struct GUF_CNT_NAME {
/* /*
- Regular iterator: base is always NULL - 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: Examples:
- rbegin(): base points to end().ptr, and ptr to end().ptr - 1 - rbegin(): base points to end().ptr, and ptr to end().ptr - 1
- rend(): base points to begin().ptr and ptr to NULL - 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) { typedef struct GUF_CAT(GUF_CNT_NAME, _iter) {
GUF_T *ptr; 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_CAT(GUF_CNT_NAME, _iter);
GUF_FN_KEYWORDS bool GUF_CAT(GUF_CNT_NAME, _valid)(const GUF_CNT_NAME* dbuf); 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); guf_err_set_if_not_null(err, GUF_ERR_NONE);
} }
/* Iterator functions */ /* 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_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)); 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; return it;
} }
if (idx < 0) { if (idx <= 0) { // begin()
it.ptr = dbuf->data; it.ptr = dbuf->data;
} else if (idx > dbuf->size) { } else if (idx >= dbuf->size) { // end()
it.ptr = dbuf->data + dbuf->size; it.ptr = dbuf->data + dbuf->size;
} else { } else {
it.ptr = dbuf->data + idx; 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; return it;
} }
if (idx < 0) { if (idx < 0) { // rend()
it.base = dbuf->data; it.base = dbuf->data;
it.ptr = NULL; it.ptr = NULL;
} else if (idx >= dbuf->size) { } else if (idx >= dbuf->size) { // rbegin()
it.base = dbuf->data + dbuf->size; it.base = dbuf->data + dbuf->size;
it.ptr = it.base - 1; it.ptr = it.base - 1;
} else { } else {

37
src/guf_int.h
View File

@ -3,42 +3,43 @@
#include "guf_common.h" #include "guf_common.h"
#include "guf_assert.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)\ #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_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_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(char, char)
GUF_DEFINE_MIN_MAX_CLAMP(int, int) 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(int8_t, i8)
GUF_DEFINE_MIN_MAX_CLAMP(int16_t, i16) GUF_DEFINE_MIN_MAX_CLAMP(int16_t, i16)
GUF_DEFINE_MIN_MAX_CLAMP(int32_t, i32) GUF_DEFINE_MIN_MAX_CLAMP(int32_t, i32)
GUF_DEFINE_MIN_MAX_CLAMP(int64_t, i64) GUF_DEFINE_MIN_MAX_CLAMP(int64_t, i64)
GUF_DEFINE_MIN_MAX_CLAMP(ptrdiff_t, ptrdiff_t) 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, 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(uint8_t, u8)
GUF_DEFINE_MIN_MAX_CLAMP(uint16_t, u16) GUF_DEFINE_MIN_MAX_CLAMP(uint16_t, u16)
GUF_DEFINE_MIN_MAX_CLAMP(uint32_t, u32) GUF_DEFINE_MIN_MAX_CLAMP(uint32_t, u32)
GUF_DEFINE_MIN_MAX_CLAMP(uint64_t, u64) GUF_DEFINE_MIN_MAX_CLAMP(uint64_t, u64)
GUF_DEFINE_MIN_MAX_CLAMP(size_t, size_t) 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 #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 #endif

61
src/guf_math.h Normal file
View File

@ -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

321
src/guf_rand.h Normal file
View File

@ -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 <math.h>
#include <float.h>
#define GUF_RAND_MAX UINT64_MAX
typedef struct guf_randstate { // State for xoshiro256** 1.0
uint64_t s[4];
} guf_randstate;
#ifdef GUF_IMPL_STATIC
#define GUF_FN_KEYWORDS static
#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

2
src/guf_sort.h
View File

@ -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
#undef GUF_IMPL_STATIC #undef GUF_IMPL_STATIC
#undef GUF_FN_KEYWORDS
#endif /* end #ifdef GUF_IMPL */ #endif /* end #ifdef GUF_IMPL */
#undef GUF_FN_KEYWORDS
#undef GUF_T #undef GUF_T
#undef GUF_FN_NAME_PREFIX #undef GUF_FN_NAME_PREFIX
#endif /* end #ifdef GUF_T */ #endif /* end #ifdef GUF_T */

46
src/guf_test.c
View File

@ -1,6 +1,7 @@
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include <stdio.h> #include <stdio.h>
#include <time.h>
#include "guf_init.h" /* Must be included once (sets up the global panic handler) */ #include "guf_init.h" /* Must be included once (sets up the global panic handler) */
#include "guf_alloc_libc.h" #include "guf_alloc_libc.h"
@ -10,6 +11,10 @@
#define GUF_T float #define GUF_T float
#include "guf_sort.h" #include "guf_sort.h"
#define GUF_IMPL_STATIC
#define GUF_T int
#include "guf_sort.h"
#define GUF_CNT_NAME dbuf_int #define GUF_CNT_NAME dbuf_int
#define GUF_T int #define GUF_T int
#define GUF_T_IS_INTEGRAL_TYPE #define GUF_T_IS_INTEGRAL_TYPE
@ -38,6 +43,9 @@
#define GUF_IMPL_STATIC #define GUF_IMPL_STATIC
#include "guf_dbuf.h" #include "guf_dbuf.h"
#define GUF_IMPL_STATIC
#include "guf_rand.h"
int main(void) int main(void)
{ {
printf("libguf test: compiled with C %ld\n", __STDC_VERSION__); printf("libguf test: compiled with C %ld\n", __STDC_VERSION__);
@ -154,5 +162,43 @@ int main(void)
} }
dbuf_int_free(&integers, NULL); 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; return EXIT_SUCCESS;
} }