415 lines
14 KiB
C
415 lines
14 KiB
C
#ifndef GUF_RAND_H
|
|
#define GUF_RAND_H
|
|
|
|
#include "guf_common.h"
|
|
|
|
#if defined(GUF_IMPL_STATIC) || defined(GUF_STATIC)
|
|
#define GUF_FN_KEYWORDS static
|
|
#else
|
|
#define GUF_FN_KEYWORDS
|
|
#endif
|
|
|
|
#ifdef GUF_RAND_32_BIT
|
|
#define GUF_RAND_MAX UINT32_MAX
|
|
typedef struct guf_randstate { // State for xoshiro128** 1.1
|
|
uint32_t s[4];
|
|
} guf_randstate;
|
|
#else
|
|
#define GUF_RAND_MAX UINT64_MAX
|
|
typedef struct guf_randstate { // State for xoshiro256** 1.0
|
|
uint64_t s[4];
|
|
} guf_randstate;
|
|
#endif
|
|
|
|
GUF_FN_KEYWORDS uint64_t guf_rand_splitmix64(uint64_t *state);
|
|
|
|
GUF_FN_KEYWORDS void guf_randstate_init(guf_randstate *state, uint64_t seed);
|
|
void guf_randstate_jump(guf_randstate *state); // Advance the state; equivalent to 2^128 calls to guf_rand_u64(state)
|
|
|
|
// uniform distributions
|
|
GUF_FN_KEYWORDS uint32_t guf_rand_u32(guf_randstate *state); // [0, UINT32_MAX]
|
|
GUF_FN_KEYWORDS uint64_t guf_rand_u64(guf_randstate *state); // [0, UINT64_MAX]
|
|
GUF_FN_KEYWORDS double guf_rand_f64(guf_randstate *state); // [0.0, 1.0)
|
|
GUF_FN_KEYWORDS float guf_rand_f32(guf_randstate *state); // [0.f, 1.f)
|
|
|
|
// return true with a probability of p, false with a probability of (1 - p)
|
|
GUF_FN_KEYWORDS bool guf_rand_bernoulli_trial_f32(guf_randstate *state, float p);
|
|
GUF_FN_KEYWORDS bool guf_rand_bernoulli_trial_f64(guf_randstate *state, double p);
|
|
GUF_FN_KEYWORDS bool guf_rand_flip(guf_randstate *state); // Fair coin flip (bernoulli trial with p == 0.5)
|
|
|
|
GUF_FN_KEYWORDS double guf_randrange_f64(guf_randstate *state, double min, double end); // [min, end)
|
|
GUF_FN_KEYWORDS float guf_randrange_f32(guf_randstate *state, float min, float end); // [min, end)
|
|
|
|
GUF_FN_KEYWORDS int32_t guf_randrange_i32(guf_randstate *state, int32_t min, int32_t max); // [min, max]
|
|
GUF_FN_KEYWORDS uint32_t guf_randrange_u32(guf_randstate *state, uint32_t min, uint32_t max); // [min, max]
|
|
GUF_FN_KEYWORDS int64_t guf_randrange_i64(guf_randstate *state, int64_t min, int64_t max); // [min, max]
|
|
|
|
// 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);
|
|
|
|
#endif
|
|
|
|
#if defined(GUF_IMPL) || defined(GUF_IMPL_STATIC)
|
|
#include <math.h>
|
|
#include <float.h>
|
|
#include "guf_common.h"
|
|
#include "guf_assert.h"
|
|
#include "guf_math.h"
|
|
|
|
/*
|
|
splitmix64 (public domain) written in 2015 by Sebastiano Vigna (vigna@acm.org)
|
|
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);
|
|
#ifdef GUF_RAND_32_BIT
|
|
for (size_t i = 0; i < GUF_STATIC_BUF_SIZE(state->s); ++i) {
|
|
state->s[i] = (uint32_t)(guf_rand_splitmix64(&seed) >> 32);
|
|
}
|
|
|
|
if (!state->s[0] && !state->s[1] && !state->s[2] && !state->s[3]) { // State must not be only zeroes:
|
|
state->s[0] = 0x9e3779b9; // arbitrary constant != 0
|
|
seed = 0x9e3779b97f4a7c15;
|
|
for (size_t i = 1; i < GUF_STATIC_BUF_SIZE(state->s); ++i) {
|
|
state->s[i] = (uint32_t)(guf_rand_splitmix64(&seed) >> 32);
|
|
}
|
|
}
|
|
#else
|
|
for (size_t i = 0; i < GUF_STATIC_BUF_SIZE(state->s); ++i) {
|
|
state->s[i] = guf_rand_splitmix64(&seed);
|
|
}
|
|
if (!state->s[0] && !state->s[1] && !state->s[2] && !state->s[3]) { // State must not be only zeroes:
|
|
state->s[0] = 0x9e3779b97f4a7c15; // arbitrary constant != 0
|
|
seed = state->s[0];
|
|
for (size_t i = 1; i < GUF_STATIC_BUF_SIZE(state->s); ++i) {
|
|
state->s[i] = guf_rand_splitmix64(&seed);
|
|
}
|
|
}
|
|
#endif
|
|
}
|
|
|
|
GUF_FN_KEYWORDS uint32_t guf_rand_u32(guf_randstate *state)
|
|
{
|
|
GUF_ASSERT(state);
|
|
GUF_ASSERT(state->s[0] || state->s[1] || state->s[2] || state->s[3]);
|
|
|
|
#ifdef GUF_RAND_32_BIT
|
|
/*
|
|
xoshiro128** 1.1 (public domain) written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org)
|
|
cf. https://prng.di.unimi.it/xoshiro128starstar.c (last-retrieved 2025-02-11)
|
|
*/
|
|
const uint32_t result = guf_rotl_u32(state->s[1] * 5, 7) * 9;
|
|
const uint32_t t = state->s[1] << 9;
|
|
state->s[2] ^= state->s[0];
|
|
state->s[3] ^= state->s[1];
|
|
state->s[1] ^= state->s[2];
|
|
state->s[0] ^= state->s[3];
|
|
state->s[2] ^= t;
|
|
state->s[3] = guf_rotl_u32(state->s[3], 11);
|
|
return result;
|
|
#else
|
|
return (uint32_t)(guf_rand_u64(state) >> 32);
|
|
#endif
|
|
}
|
|
|
|
GUF_FN_KEYWORDS uint64_t guf_rand_u64(guf_randstate *state)
|
|
{
|
|
GUF_ASSERT(state);
|
|
GUF_ASSERT(state->s[0] || state->s[1] || state->s[2] || state->s[3]);
|
|
|
|
#ifdef GUF_RAND_32_BIT
|
|
const uint32_t lower_bits = guf_rand_u32(state);
|
|
const uint32_t upper_bits = guf_rand_u32(state);
|
|
return ((uint64_t)upper_bits << 32) | (uint64_t)lower_bits; // TODO: not sure if that's a good idea...
|
|
#else
|
|
/*
|
|
xoshiro256** 1.0 (public domain) written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org)
|
|
cf. https://prng.di.unimi.it/xoshiro256starstar.c (last-retrieved 2025-02-11)
|
|
*/
|
|
const uint64_t result = guf_rotl_u64(state->s[1] * 5, 7) * 9;
|
|
const uint64_t t = state->s[1] << 17;
|
|
state->s[2] ^= state->s[0];
|
|
state->s[3] ^= state->s[1];
|
|
state->s[1] ^= state->s[2];
|
|
state->s[0] ^= state->s[3];
|
|
state->s[2] ^= t;
|
|
state->s[3] = guf_rotl_u64(state->s[3], 45);
|
|
return result;
|
|
#endif
|
|
}
|
|
|
|
/*
|
|
Equivalent to 2^128 calls to guf_rand() (or 2^64 calls if GUF_RAND_32_BIT); it can be used to generate 2^128 (or 2^64)
|
|
non-overlapping subsequences for parallel computations.
|
|
*/
|
|
void guf_randstate_jump(guf_randstate *state)
|
|
{
|
|
GUF_ASSERT(state);
|
|
#ifdef GUF_RAND_32_BIT
|
|
static const uint32_t JUMP[] = { 0x8764000b, 0xf542d2d3, 0x6fa035c3, 0x77f2db5b };
|
|
uint32_t s0 = 0;
|
|
uint32_t s1 = 0;
|
|
uint32_t s2 = 0;
|
|
uint32_t s3 = 0;
|
|
for(size_t i = 0; i < sizeof JUMP / sizeof *JUMP; ++i) {
|
|
for(int b = 0; b < 32; ++b) {
|
|
if (JUMP[i] & UINT32_C(1) << b) {
|
|
s0 ^= state->s[0];
|
|
s1 ^= state->s[1];
|
|
s2 ^= state->s[2];
|
|
s3 ^= state->s[3];
|
|
}
|
|
guf_rand_u32(state);
|
|
}
|
|
}
|
|
state->s[0] = s0;
|
|
state->s[1] = s1;
|
|
state->s[2] = s2;
|
|
state->s[3] = s3;
|
|
#else
|
|
static const uint64_t JUMP[] = { 0x180ec6d33cfd0aba, 0xd5a61266f0c9392c, 0xa9582618e03fc9aa, 0x39abdc4529b1661c };
|
|
uint64_t s0 = 0;
|
|
uint64_t s1 = 0;
|
|
uint64_t s2 = 0;
|
|
uint64_t s3 = 0;
|
|
for (size_t i = 0; i < sizeof JUMP / sizeof *JUMP; ++i) {
|
|
for (int b = 0; b < 64; ++b) {
|
|
if (JUMP[i] & UINT64_C(1) << b) {
|
|
s0 ^= state->s[0];
|
|
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;
|
|
#endif
|
|
}
|
|
|
|
// 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)
|
|
{
|
|
#ifdef GUF_RAND_32_BIT
|
|
return (guf_rand_u32(state) >> 8) * 0x1.0p-24f; // 8 == 32 - 24; (float has a 24-bit mantissa/significand)
|
|
#else
|
|
return (guf_rand_u64(state) >> 40) * 0x1.0p-24f; // 40 == 64 - 24; (float has a 24-bit mantissa/significand)
|
|
#endif
|
|
}
|
|
|
|
GUF_FN_KEYWORDS bool guf_rand_bernoulli_trial_f32(guf_randstate *state, float p)
|
|
{
|
|
p = guf_clamp_f32(p, 0, 1);
|
|
return guf_rand_f32(state) < p; // never true for p = 0, always true for p = 1 since guf_rand_f64 is in range [0, 1)
|
|
}
|
|
|
|
GUF_FN_KEYWORDS bool guf_rand_bernoulli_trial_f64(guf_randstate *state, double p)
|
|
{
|
|
p = guf_clamp_f64(p, 0, 1);
|
|
return guf_rand_f64(state) < p; // never true for p = 0, always true for p = 1 since guf_rand_f64 is in range [0, 1)
|
|
}
|
|
|
|
GUF_FN_KEYWORDS bool guf_rand_flip(guf_randstate *state)
|
|
{
|
|
#ifdef GUF_RAND_32_BIT
|
|
return guf_rand_bernoulli_trial_f32(state, 0.5f);
|
|
#else
|
|
return guf_rand_bernoulli_trial_f64(state, 0.5);
|
|
#endif
|
|
}
|
|
|
|
// returns uniformly-distributed random double in range [min, end) (or min if min == end)
|
|
GUF_FN_KEYWORDS double guf_randrange_f64(guf_randstate *state, double min, double end)
|
|
{
|
|
if (min == (double)INFINITY) {
|
|
min = DBL_MAX;
|
|
} else if (min == (double)-INFINITY) {
|
|
min = -DBL_MAX;
|
|
}
|
|
if (end == (double)INFINITY) {
|
|
end = DBL_MAX;
|
|
} else if (end == (double)-INFINITY) {
|
|
end = -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 (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_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 = (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;
|
|
}
|
|
|
|
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 = (double)max - (double)min;
|
|
const double result = floor(guf_rand_f64(state) * (delta + 1.0) + min);
|
|
GUF_ASSERT(result >= min && result <= max);
|
|
return (uint32_t)result;
|
|
}
|
|
|
|
// returns uniformly-distributed random int64_t in range [min, max] (max is inclusive as opposed to the f32/f64 versions)
|
|
GUF_FN_KEYWORDS int64_t guf_randrange_i64(guf_randstate *state, int64_t min, int64_t max)
|
|
{
|
|
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;
|
|
}
|
|
|
|
|
|
// 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 /* endif GUF_IMPL/GUF_IMPL_STATIC */
|
|
|
|
#undef GUF_STATIC
|
|
#undef GUF_FN_KEYWORDS
|
|
#undef GUF_RAND_32_BIT
|