From d9905d3300eb5c9b6706688f74a42bac6920de1c Mon Sep 17 00:00:00 2001 From: Tyge Løvset Date: Wed, 16 Dec 2020 12:49:02 +0100 Subject: Minor fixing of unsafe macros. --- docs/crandom_api.md | 12 ++++++------ stc/crandom.h | 35 +++++++++++++++++++---------------- 2 files changed, 25 insertions(+), 22 deletions(-) diff --git a/docs/crandom_api.md b/docs/crandom_api.md index 53e8a0d1..c4b6fa98 100644 --- a/docs/crandom_api.md +++ b/docs/crandom_api.md @@ -5,12 +5,12 @@ This describes the API of module **crand**. It contains *pcg32* created by Melis ## Types -| Name | Type definition | Used to represent... | -|:-------------------|:-------------------------------------------|:-----------------------------| -| `crand_t` | `struct {uint64_t state[4];}` | The PRNG engine type | -| `crand_uniform_t` | `struct {int64_t offset; uint64_t range;}` | Integer uniform distribution | -| `crand_uniformf_t` | `struct {double offset, range;}` | Real number uniform distr. | -| `crand_normalf_t` | `struct {double mean, stddev;}` | Normal distribution type | +| Name | Type definition | Used to represent... | +|:-------------------|:------------------------------------------|:-----------------------------| +| `crand_t` | `struct {uint64_t state[4];}` | The PRNG engine type | +| `crand_uniform_t` | `struct {int64_t lower; uint64_t range;}` | Integer uniform distribution | +| `crand_uniformf_t` | `struct {double lower, range;}` | Real number uniform distr. | +| `crand_normalf_t` | `struct {double mean, stddev;}` | Normal distribution type | ## Header file diff --git a/stc/crandom.h b/stc/crandom.h index 5db3c160..8a67e3ce 100644 --- a/stc/crandom.h +++ b/stc/crandom.h @@ -43,10 +43,10 @@ int main() { #include #include -typedef struct {uint64_t state[4];} crand_t; -typedef struct {int64_t offset; uint64_t range, threshold;} crand_uniform_t; -typedef struct {double offset, range;} crand_uniformf_t; -typedef struct {double mean, stddev, next; bool has_next;} crand_normalf_t; +typedef struct {uint64_t state[4];} crand_t; +typedef struct {int64_t lower; uint64_t range, threshold;} crand_uniform_t; +typedef struct {double lower, range;} crand_uniformf_t; +typedef struct {double mean, stddev, next; bool has_next;} crand_normalf_t; /* int random number generator, range [0, 2^64). PRNG copyright Tyge Løvset, NORCE Research, 2020 */ @@ -69,7 +69,7 @@ STC_INLINE crand_uniformf_t crand_uniformf_init(double low, double high) { crand_uniformf_t dist = {low, high - low}; return dist; } STC_INLINE double crand_uniformf(crand_t* rng, crand_uniformf_t* dist) { - return dist->offset + crand_nextf(rng) * dist->range; + return crand_nextf(rng)*dist->range + dist->lower; } /* double normal distributed RNG. */ @@ -108,8 +108,7 @@ STC_DEF uint64_t crand_next(crand_t* rng) { return result; } -/* unbiased integer random number generator in range [low, high] */ - +/* Very fast unbiased uniform int RNG with bounds [low, high] */ STC_DEF crand_uniform_t crand_uniform_init(int64_t low, int64_t high) { crand_uniform_t dist = {low, (uint64_t) (high - low + 1)}; dist.threshold = (uint64_t)(-dist.range) % dist.range; @@ -117,23 +116,27 @@ STC_DEF crand_uniform_t crand_uniform_init(int64_t low, int64_t high) { } #if defined(__SIZEOF_INT128__) - #define cmul128(a, b, m) *(__uint128_t *)m = (__uint128_t)a * b + #define cmul128(a, b, lo, hi) \ + do { __uint128_t _z = (__uint128_t)(a) * (b); \ + *(lo) = (uint64_t)_z, *(hi) = _z >> 64; } while(0) #elif defined(_MSC_VER) && defined(_WIN64) #include - #define cmul128(a, b, m) m[0] = _umul128(a, b, &m[1]) + #define cmul128(a, b, lo, hi) (*(lo) = _umul128(a, b, hi), (void)0) #elif defined(__x86_64__) - #define cmul128(a, b, m) \ - asm("mulq %[rhs]" : "=a" (m[0]), "=d" (m[1]) \ + #define cmul128(a, b, lo, hi) \ + asm("mulq %[rhs]" : "=a" (*(lo)), "=d" (*(hi)) \ : [lhs] "0" (a), [rhs] "rm" (b)) #endif -STC_DEF int64_t crand_uniform(crand_t* rng, crand_uniform_t* dist) { - uint64_t m[2]; - do { cmul128(crand_next(rng), dist->range, m); } while (m[0] < dist->threshold); - return dist->offset + m[1]; +STC_DEF int64_t crand_uniform(crand_t* rng, crand_uniform_t* d) { + uint64_t lo, hi; + do { + cmul128(crand_next(rng), d->range, &lo, &hi); + } while (lo < d->threshold); + return d->lower + hi; } -/* Marsaglia polar method for gaussian distribution. */ +/* Marsaglia polar method for gaussian/normal distribution. */ STC_DEF double crand_normalf(crand_t* rng, crand_normalf_t* dist) { double u1, u2, s, m; if (dist->has_next) { -- cgit v1.2.3