star-sp

Random number generators and distributions
git clone git://git.meso-star.fr/star-sp.git
Log | Files | Refs | README | LICENSE

commit 1aa2de5d96259fa4324bf8f4fc99782a5aaa345d
parent 24b84a59cb545064f176f3fed05822658103c954
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu,  2 Apr 2015 14:14:16 +0200

Add and test the uniform_<int|double> RNG distributions

Diffstat:
Msrc/ssp.h | 25++++++++++++++++++++++---
Msrc/ssp_rng.c | 85++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------
Msrc/test_ssp_rng.c | 22+++++++++++++++++++++-
3 files changed, 121 insertions(+), 11 deletions(-)

diff --git a/src/ssp.h b/src/ssp.h @@ -63,8 +63,13 @@ struct ssp_rng_type { void (*init)(struct mem_allocator* allocator, void* state); void (*release)(void* state); void (*set)(void* state, const unsigned long seed); + unsigned long (*get)(void* state); - double (*get_canonical)(void* state); + unsigned long (*uniform_int) + (void* state, const unsigned long min, const unsigned long max); + double (*uniform_double)(void* state, const double min, const double max); + double (*canonical)(void* state); + unsigned long min; unsigned long max; size_t sizeof_state; @@ -105,9 +110,23 @@ SSP_API unsigned long ssp_rng_get (struct ssp_rng* rng); -/* Uniform random distribution in [0, 1) */ +/* Uniform random integer distribution in [lower, upper] */ +SSP_API unsigned long +ssp_rng_uniform_int + (struct ssp_rng* rng, + const unsigned long lower, + const unsigned long upper); + +/* Uniform random floating point distribution in [lower, upper] */ +SSP_API double +ssp_rng_uniform_double + (struct ssp_rng* rng, + const double lower, + const double upper); + +/* Uniform random floating point distribution in [0, 1) */ SSP_API double -ssp_rng_get_canonical +ssp_rng_canonical (struct ssp_rng* rng); SSP_API unsigned long diff --git a/src/ssp_rng.c b/src/ssp_rng.c @@ -78,8 +78,27 @@ rng_kiss_get(void* data) return (uint32_t)(kiss->x + kiss->y + kiss->z); } +static unsigned long +rng_kiss_uniform_int + (void* data, const unsigned long lower, const unsigned long upper) +{ + ASSERT(lower <= upper); + return (uint32_t) + ( (double)rng_kiss_get(data)/(double)UINT32_MAX + * (double)(upper - lower) + (double)lower); +} + +static double +rng_kiss_uniform_double(void* data, const double lower, const double upper) +{ + ASSERT(lower <= upper); + ASSERT(lower <= upper); + return (double)rng_kiss_get(data)/(double)UINT32_MAX * (upper-lower) + lower; +} + + static double -rng_kiss_get_canonical(void* data) +rng_kiss_canonical(void* data) { return (double)rng_kiss_get(data) / ((double)UINT32_MAX + 1); } @@ -103,7 +122,9 @@ struct ssp_rng_type ssp_rng_kiss = { rng_kiss_release, rng_kiss_set, rng_kiss_get, - rng_kiss_get_canonical, + rng_kiss_uniform_int, + rng_kiss_uniform_double, + rng_kiss_canonical, 0, UINT32_MAX, sizeof(struct rng_kiss) @@ -128,8 +149,26 @@ rng_mt19937_64_get(void* data) return (*mt)(); } +static unsigned long +rng_mt19937_64_uniform_int + (void* data, const unsigned long lower, const unsigned long upper) +{ + std::mt19937_64* mt = (std::mt19937_64*)data; + ASSERT(mt && lower <= upper); + return std::uniform_int_distribution<unsigned long>(lower, upper)(*mt); +} + +static double +rng_mt19937_64_uniform_double + (void* data, const double lower, const double upper) +{ + std::mt19937_64* mt = (std::mt19937_64*)data; + ASSERT(mt && lower <= upper); + return std::uniform_real_distribution<>(lower, upper)(*mt); +} + static double -rng_mt19937_64_get_canonical(void* data) +rng_mt19937_64_canonical(void* data) { std::mt19937_64* mt = (std::mt19937_64*)data; ASSERT(mt); @@ -158,7 +197,9 @@ struct ssp_rng_type ssp_rng_mt19937_64 = { rng_mt19937_64_release, rng_mt19937_64_set, rng_mt19937_64_get, - rng_mt19937_64_get_canonical, + rng_mt19937_64_uniform_int, + rng_mt19937_64_uniform_double, + rng_mt19937_64_canonical, std::mt19937_64::min(), std::mt19937_64::max(), sizeof(std::mt19937_64) @@ -167,6 +208,20 @@ struct ssp_rng_type ssp_rng_mt19937_64 = { /******************************************************************************* * Helper functions ******************************************************************************/ +static char +check_type(const struct ssp_rng_type* type) +{ + ASSERT(type); + return type->init != NULL + && type->release != NULL + && type->set != NULL + && type->get != NULL + && type->uniform_int != NULL + && type->uniform_double != NULL + && type->canonical != NULL + && type->min <= type->max; +} + static void rng_release(ref_T* ref) { @@ -193,7 +248,7 @@ ssp_rng_create struct ssp_rng* rng = NULL; res_T res = RES_OK; - if(!type || !out_rng) { + if(!type || !out_rng || !check_type(type)) { res = RES_BAD_ARG; goto error; } @@ -249,11 +304,27 @@ ssp_rng_get(struct ssp_rng* rng) return rng->type.get(rng->state); } +unsigned long +ssp_rng_uniform_int + (struct ssp_rng* rng, const unsigned long lower, const unsigned long upper) +{ + if(!rng) FATAL("The Random Number Generator is NULL\n"); + return rng->type.uniform_int(rng->state, lower, upper); +} + +double +ssp_rng_uniform_double + (struct ssp_rng* rng, double lower, double upper) +{ + if(!rng) FATAL("The Random Number Generator is NULL\n"); + return rng->type.uniform_double(rng->state, lower, upper); +} + double -ssp_rng_get_canonical(struct ssp_rng* rng) +ssp_rng_canonical(struct ssp_rng* rng) { if(!rng) FATAL("The Random Number Generator is NULL\n"); - return rng->type.get_canonical(rng->state); + return rng->type.canonical(rng->state); } unsigned long diff --git a/src/test_ssp_rng.c b/src/test_ssp_rng.c @@ -87,7 +87,27 @@ test_rng(struct ssp_rng_type* type) } FOR_EACH(i, 0, NRAND) { - dataf[i] = ssp_rng_get_canonical(rng); + datai0[i] = ssp_rng_uniform_int(rng, 1, 79); + CHECK(datai0[i] >= 1, 1); + CHECK(datai0[i] <= 79, 1); + } + FOR_EACH(i, 0, NRAND) { + FOR_EACH(j, 0, NRAND) if(datai0[i] != datai0[j]) break; + NCHECK(j, NRAND); + } + + FOR_EACH(i, 0, NRAND) { + dataf[i] = ssp_rng_uniform_double(rng, -5.0, 11.3); + CHECK(dataf[i] >= -5.0, 1); + CHECK(dataf[i] <= 11.3, 1); + } + FOR_EACH(i, 0, NRAND) { + FOR_EACH(j, 0, NRAND) if(dataf[i] != dataf[j]) break; + NCHECK(j, NRAND); + } + + FOR_EACH(i, 0, NRAND) { + dataf[i] = ssp_rng_canonical(rng); CHECK(dataf[i] >= 0.0, 1); CHECK(dataf[i] < 1.0, 1); FOR_EACH(j, 0, i) {