star-sp

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

commit 32bd4ecec1dd9f70dd671dac8cb8d839890f9ae7
parent c1a0faaf2474f8e688f6f6035adfe8db4f9fd5c6
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed,  1 Apr 2015 17:41:16 +0200

Update the RNG API

Diffstat:
Mcmake/CMakeLists.txt | 2+-
Msrc/ssp.h | 42++++++++++++++++++++++++++++++------------
Msrc/ssp_rng.c | 187+++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------
Dsrc/test_ssp_rng.c | 97-------------------------------------------------------------------------------
Asrc/test_ssp_rng_kiss.c | 102+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
5 files changed, 264 insertions(+), 166 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -94,7 +94,7 @@ if(NOT NO_TEST) add_test(${_name} ${_name}) endfunction(new_test) - new_test(test_ssp_rng) + new_test(test_ssp_rng_kiss) endif() ################################################################################ diff --git a/src/ssp.h b/src/ssp.h @@ -52,21 +52,39 @@ #define SSP(Func) ssp_ ## Func #endif -enum ssp_rng_type { - SSP_RNG_KISS, - SSP_RNG_MT19937_64 -}; - +/* Forward declaration of opaque types */ struct ssp_rng; +/* Forward declaration of external types */ struct mem_allocator; +/* Generic Random Number Generator type descriptor */ +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 min; + unsigned long max; + size_t sizeof_state; +}; + BEGIN_DECLS +/* David Jones's Keep It Simple Stupid builtin PRNG type. Suitable for fast + * basic randomness */ +SSP_API struct ssp_rng_type ssp_rng_kiss; +/* 64-bits Mersenne Twister builtin PRNG type of Matsumoto and Nishimura, 2000 */ +SSP_API struct ssp_rng_type ssp_rng_mt19937_64; + +/******************************************************************************* + * API declaration + ******************************************************************************/ SSP_API res_T ssp_rng_create (struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */ - const enum ssp_rng_type type, + const struct ssp_rng_type* type, struct ssp_rng** rng); SSP_API res_T @@ -82,10 +100,16 @@ ssp_rng_set (struct ssp_rng* rng, const unsigned long seed); +/* Uniform random distribution in [ssp_rng_min(rng), ssp_rng_max(rng)] */ SSP_API unsigned long ssp_rng_get (struct ssp_rng* rng); +/* Uniform random distribution in [0, 1) */ +SSP_API double +ssp_rng_get_canonical + (struct ssp_rng* rng); + SSP_API unsigned long ssp_rng_min (struct ssp_rng* rng); @@ -94,12 +118,6 @@ SSP_API unsigned long ssp_rng_max (struct ssp_rng* rng); -static FINLINE double /* in [0, 1) */ -ssp_rng_get_canonical(struct ssp_rng* rng) -{ - return (double)ssp_rng_get(rng) / ((double)ssp_rng_max(rng) + 1.0); -} - END_DECLS #endif /* SSP_H */ diff --git a/src/ssp_rng.c b/src/ssp_rng.c @@ -36,27 +36,35 @@ #include <random> -struct rng_kiss { - uint32_t x, y, z, c; -}; - struct ssp_rng { - enum ssp_rng_type type; - union { - struct rng_kiss kiss; - std::mt19937 mt19937_64; - } algo; - + struct ssp_rng_type type; + void* state; struct mem_allocator* allocator; ref_T ref; }; /******************************************************************************* - * Helper functions + * KISS PRNG ******************************************************************************/ -static FINLINE uint32_t -rng_kiss_get(struct rng_kiss* kiss) +struct rng_kiss { uint32_t x, y, z, c; }; + +static void +rng_kiss_set(void* data, const unsigned long seed) { + struct rng_kiss* kiss = (struct rng_kiss*)data; + std::mt19937 rng_mt(seed % ((1L<<32)-1)); + ASSERT(kiss); + kiss->x = (uint32_t)rng_mt(); + while(!(kiss->y = (uint32_t)rng_mt())); /* y must be != 0 */ + kiss->z = (uint32_t)rng_mt(); + /* Offset c by 1 to avoid z=c=0; should be less than 698769069 */ + kiss->c = (uint32_t)rng_mt() % 698769068 + 1; +} + +static unsigned long +rng_kiss_get(void* data) +{ + struct rng_kiss* kiss = (struct rng_kiss*)data; uint64_t t; ASSERT(kiss); @@ -70,28 +78,106 @@ rng_kiss_get(struct rng_kiss* kiss) return (uint32_t)(kiss->x + kiss->y + kiss->z); } -static void -rng_kiss_set(struct rng_kiss* kiss, const unsigned long seed) +static double +rng_kiss_get_canonical(void* data) { - std::mt19937 rng_mt(seed % ((1L<<32)-1)); + struct rng_kiss* kiss = (struct rng_kiss*)data; ASSERT(kiss); - kiss->x = (uint32_t)rng_mt(); - while(!(kiss->y = (uint32_t)rng_mt())); /* y must be != 0 */ - kiss->z = (uint32_t)rng_mt(); - /* Offset c by 1 to avoid z=c=0; should be less than 698769069 */ - kiss->c = (uint32_t)rng_mt() % 698769068 + 1; + return (double)rng_kiss_get(data) / ((double)UINT32_MAX + 1); +} + +static void +rng_kiss_init(struct mem_allocator* allocator, void* data) +{ + (void)allocator; + rng_kiss_set(data, 0); +} + +static void +rng_kiss_release(void* data) +{ + (void)data; +} + +/* Exported type */ +struct ssp_rng_type ssp_rng_kiss = { + rng_kiss_init, + rng_kiss_release, + rng_kiss_set, + rng_kiss_get, + rng_kiss_get_canonical, + 0, + UINT32_MAX, + sizeof(struct rng_kiss) +}; + +/******************************************************************************* + * 64-bits Mersenne Twister PRNG + ******************************************************************************/ +static void +rng_mt19937_64_set(void* data, const unsigned long seed) +{ + std::mt19937_64* mt = (std::mt19937_64*)data; + ASSERT(mt); + mt->seed(seed); +} + +static unsigned long +rng_mt19937_64_get(void* data) +{ + std::mt19937_64* mt = (std::mt19937_64*)data; + ASSERT(mt); + return (*mt)(); +} + +static double +rng_mt19937_64_get_canonical(void* data) +{ + std::mt19937_64* mt = (std::mt19937_64*)data; + ASSERT(mt); + return std::generate_canonical<double, 48>(*mt); +} + +static void +rng_mt19937_64_init(struct mem_allocator* allocator, void* data) +{ + (void)allocator; + ASSERT(data); + new (data) std::mt19937_64; +} + +static void +rng_mt19937_64_release(void* data) +{ + std::mt19937_64* mt = (std::mt19937_64*)data; + ASSERT(mt); + mt->~mersenne_twister_engine(); } +/* Exported type */ +struct ssp_rng_type ssp_rng_mt19937_64 = { + rng_mt19937_64_init, + rng_mt19937_64_release, + rng_mt19937_64_set, + rng_mt19937_64_get, + rng_mt19937_64_get_canonical, + std::mt19937_64::min(), + std::mt19937_64::max(), + sizeof(std::mt19937_64) +}; + +/******************************************************************************* + * Helper functions + ******************************************************************************/ static void rng_release(ref_T* ref) { struct ssp_rng* rng; ASSERT(ref); rng = CONTAINER_OF(ref, struct ssp_rng, ref); - switch(rng->type) { - case SSP_RNG_KISS: /* Do nothing */ break; - case SSP_RNG_MT19937_64: rng->algo.mt19937_64.~mersenne_twister_engine(); break; - default: FATAL("Unreachable code\n"); break; + if(rng->state) { + rng->type.release(rng->state); + MEM_FREE(rng->allocator, rng->state); } MEM_FREE(rng->allocator, rng); } @@ -102,14 +188,14 @@ rng_release(ref_T* ref) res_T ssp_rng_create (struct mem_allocator* mem_allocator, - const enum ssp_rng_type type, + const struct ssp_rng_type* type, struct ssp_rng** out_rng) { struct mem_allocator* allocator; struct ssp_rng* rng = NULL; res_T res = RES_OK; - if(!out_rng) { + if(!type || !out_rng) { res = RES_BAD_ARG; goto error; } @@ -117,18 +203,19 @@ ssp_rng_create allocator = mem_allocator ? mem_allocator : &mem_default_allocator; rng = (struct ssp_rng*)MEM_CALLOC(allocator, 1, sizeof(struct ssp_rng)); if(!rng) { - res = RES_BAD_ARG; + res = RES_MEM_ERR; goto error; } rng->allocator = allocator; ref_init(&rng->ref); - switch(type) { - case SSP_RNG_KISS: rng_kiss_set(&rng->algo.kiss, 0); break; - case SSP_RNG_MT19937_64: new (&rng->algo.mt19937_64) std::mt19937_64; break; - default: res = RES_BAD_ARG; goto error; + rng->state = MEM_CALLOC(rng->allocator, 1, sizeof(type->sizeof_state)); + if(!rng->state) { + res = RES_MEM_ERR; + goto error; } - rng->type = type; + type->init(allocator, rng->state); + rng->type = *type; exit: if(out_rng) *out_rng = rng; @@ -161,47 +248,35 @@ unsigned long ssp_rng_get(struct ssp_rng* rng) { if(!rng) FATAL("The Random Number Generator is NULL\n"); - switch(rng->type) { - case SSP_RNG_KISS: return rng_kiss_get(&rng->algo.kiss); break; - case SSP_RNG_MT19937_64: return rng->algo.mt19937_64(); break; - default: FATAL("Unexpected Random Number Generator type\n"); break; - } - return 0; + return rng->type.get(rng->state); +} + +double +ssp_rng_get_canonical(struct ssp_rng* rng) +{ + if(!rng) FATAL("The Random Number Generator is NULL\n"); + return rng->type.get_canonical(rng->state); } unsigned long ssp_rng_min(struct ssp_rng* rng) { if(!rng) FATAL("The Random Number Generator is NULL\n"); - switch(rng->type) { - case SSP_RNG_KISS: return 0; - case SSP_RNG_MT19937_64: return std::mt19937_64::min(); - default: FATAL("Unexpected Random Number Generator type\n"); break; - } - return 0; + return rng->type.min; } unsigned long ssp_rng_max(struct ssp_rng* rng) { if(!rng) FATAL("The Random Number Generator is NULL\n"); - switch(rng->type) { - case SSP_RNG_KISS: return UINT32_MAX; - case SSP_RNG_MT19937_64: return std::mt19937_64::max(); - default: FATAL("Unexpected Random Number Generator type\n"); break; - } - return 0; + return rng->type.max; } res_T ssp_rng_set(struct ssp_rng* rng, const unsigned long seed) { if(!rng) return RES_BAD_ARG; - switch(rng->type) { - case SSP_RNG_KISS: rng_kiss_set(&rng->algo.kiss, (uint32_t)seed); break; - case SSP_RNG_MT19937_64: rng->algo.mt19937_64.seed(seed); break; - default: return RES_BAD_ARG; - } + rng->type.set(rng->state, seed); return RES_OK; } diff --git a/src/test_ssp_rng.c b/src/test_ssp_rng.c @@ -1,97 +0,0 @@ -/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com) - * - * This software is a program whose purpose is to test the spp library. - * - * This software is governed by the CeCILL-C license under French law and - * abiding by the rules of distribution of free software. You can use, - * modify and/or redistribute the software under the terms of the CeCILL-C - * license as circulated by CEA, CNRS and INRIA at the following URL - * "http://www.cecill.info". - * - * As a counterpart to the access to the source code and rights to copy, - * modify and redistribute granted by the license, users are provided only - * with a limited warranty and the software's author, the holder of the - * economic rights, and the successive licensors have only limited - * liability. - * - * In this respect, the user's attention is drawn to the risks associated - * with loading, using, modifying and/or developing or reproducing the - * software by the user in light of its specific status of free software, - * that may mean that it is complicated to manipulate, and that also - * therefore means that it is reserved for developers and experienced - * professionals having in-depth computer knowledge. Users are therefore - * encouraged to load and test the software's suitability as regards their - * requirements in conditions enabling the security of their systems and/or - * data to be ensured and, more generally, to use and operate it in the - * same conditions as regards security. - * - * The fact that you are presently reading this means that you have had - * knowledge of the CeCILL-C license and that you accept its terms. */ - -#include "ssp.h" -#include "test_ssp_utils.h" - -#define NRAND 1024 - -int -main(int argc, char** argv) -{ - struct ssp_rng* rng; - struct mem_allocator allocator; - unsigned long datai0[NRAND]; - unsigned long datai1[NRAND]; - double dataf[NRAND]; - int i, j; - (void)argc, (void)argv; - - mem_init_proxy_allocator(&allocator, &mem_default_allocator); - - CHECK(ssp_rng_create(NULL, SSP_RNG_KISS, NULL), RES_BAD_ARG); - CHECK(ssp_rng_create(&allocator, SSP_RNG_KISS, NULL), RES_BAD_ARG); - CHECK(ssp_rng_create(NULL, SSP_RNG_KISS, &rng), RES_OK); - - CHECK(ssp_rng_ref_get(NULL), RES_BAD_ARG); - CHECK(ssp_rng_ref_get(rng), RES_OK); - CHECK(ssp_rng_ref_put(NULL), RES_BAD_ARG); - CHECK(ssp_rng_ref_put(rng), RES_OK); - CHECK(ssp_rng_ref_put(rng), RES_OK); - - CHECK(ssp_rng_create(&allocator, SSP_RNG_KISS, &rng), RES_OK); - CHECK(ssp_rng_set(NULL, 0), RES_BAD_ARG); - CHECK(ssp_rng_set(rng, 0), RES_OK); - - FOR_EACH(i, 0, NRAND) datai0[i] = ssp_rng_get(rng); - FOR_EACH(i, 0, NRAND) { - FOR_EACH(j, 0, NRAND) { - if(i == j) continue; - NCHECK(datai0[i], datai0[j]); - } - } - CHECK(ssp_rng_set(rng, 0xDECAFBAD), RES_OK); - FOR_EACH(i, 0, NRAND) datai1[i] = ssp_rng_get(rng); - FOR_EACH(i, 0, NRAND) { - NCHECK(datai0[i], datai1[i]); - FOR_EACH(j, 0, NRAND) { - if(i == j) continue; - NCHECK(datai1[i], datai1[j]); - } - } - FOR_EACH(i, 0, NRAND) - dataf[i] = ssp_rng_get_canonical(rng); - FOR_EACH(i, 0, NRAND) { - CHECK(dataf[i] >= 0.0, 1); - CHECK(dataf[i] < 1.0, 1); - FOR_EACH(j, 0, NRAND) { - if(i == j) continue; - NCHECK(dataf[i], dataf[j]); - } - } - CHECK(ssp_rng_min(rng), 0); - CHECK(ssp_rng_max(rng), UINT32_MAX); - CHECK(ssp_rng_ref_put(rng), RES_OK); - - check_memory_allocator(&allocator); - mem_shutdown_proxy_allocator(&allocator); - CHECK(mem_allocated_size(), 0); - return 0; -} diff --git a/src/test_ssp_rng_kiss.c b/src/test_ssp_rng_kiss.c @@ -0,0 +1,102 @@ +/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com) + * + * This software is a program whose purpose is to test the spp library. + * + * This software is governed by the CeCILL-C license under French law and + * abiding by the rules of distribution of free software. You can use, + * modify and/or redistribute the software under the terms of the CeCILL-C + * license as circulated by CEA, CNRS and INRIA at the following URL + * "http://www.cecill.info". + * + * As a counterpart to the access to the source code and rights to copy, + * modify and redistribute granted by the license, users are provided only + * with a limited warranty and the software's author, the holder of the + * economic rights, and the successive licensors have only limited + * liability. + * + * In this respect, the user's attention is drawn to the risks associated + * with loading, using, modifying and/or developing or reproducing the + * software by the user in light of its specific status of free software, + * that may mean that it is complicated to manipulate, and that also + * therefore means that it is reserved for developers and experienced + * professionals having in-depth computer knowledge. Users are therefore + * encouraged to load and test the software's suitability as regards their + * requirements in conditions enabling the security of their systems and/or + * data to be ensured and, more generally, to use and operate it in the + * same conditions as regards security. + * + * The fact that you are presently reading this means that you have had + * knowledge of the CeCILL-C license and that you accept its terms. */ + +#include "ssp.h" +#include "test_ssp_utils.h" + +#define NRAND 1024 + +int +main(int argc, char** argv) +{ + struct ssp_rng* rng; + struct mem_allocator allocator; + unsigned long datai0[NRAND]; + unsigned long datai1[NRAND]; + double dataf[NRAND]; + int i, j; + (void)argc, (void)argv; + + mem_init_proxy_allocator(&allocator, &mem_default_allocator); + + CHECK(ssp_rng_create(NULL, NULL, NULL), RES_BAD_ARG); + CHECK(ssp_rng_create(&allocator, NULL, NULL), RES_BAD_ARG); + CHECK(ssp_rng_create(NULL, &ssp_rng_kiss, NULL), RES_BAD_ARG); + CHECK(ssp_rng_create(&allocator, &ssp_rng_kiss, NULL), RES_BAD_ARG); + CHECK(ssp_rng_create(NULL, NULL, &rng), RES_BAD_ARG); + CHECK(ssp_rng_create(&allocator, NULL, &rng), RES_BAD_ARG); + CHECK(ssp_rng_create(NULL, &ssp_rng_kiss, &rng), RES_OK); + + CHECK(ssp_rng_ref_get(NULL), RES_BAD_ARG); + CHECK(ssp_rng_ref_get(rng), RES_OK); + CHECK(ssp_rng_ref_put(NULL), RES_BAD_ARG); + CHECK(ssp_rng_ref_put(rng), RES_OK); + CHECK(ssp_rng_ref_put(rng), RES_OK); + + CHECK(ssp_rng_create(&allocator, &ssp_rng_kiss, &rng), RES_OK); + CHECK(ssp_rng_set(NULL, 0), RES_BAD_ARG); + CHECK(ssp_rng_set(rng, 0), RES_OK); + + FOR_EACH(i, 0, NRAND) datai0[i] = ssp_rng_get(rng); + FOR_EACH(i, 0, NRAND) { + FOR_EACH(j, 0, NRAND) { + if(i == j) continue; + NCHECK(datai0[i], datai0[j]); + } + } + CHECK(ssp_rng_set(rng, 0xDECAFBAD), RES_OK); + FOR_EACH(i, 0, NRAND) datai1[i] = ssp_rng_get(rng); + FOR_EACH(i, 0, NRAND) { + NCHECK(datai0[i], datai1[i]); + FOR_EACH(j, 0, NRAND) { + if(i == j) continue; + NCHECK(datai1[i], datai1[j]); + } + } + FOR_EACH(i, 0, NRAND) + dataf[i] = ssp_rng_get_canonical(rng); + FOR_EACH(i, 0, NRAND) { + CHECK(dataf[i] >= 0.0, 1); + CHECK(dataf[i] < 1.0, 1); + FOR_EACH(j, 0, NRAND) { + if(i == j) continue; + NCHECK(dataf[i], dataf[j]); + } + } + CHECK(ssp_rng_min(rng), 0); + CHECK(ssp_rng_max(rng), UINT32_MAX); + CHECK(ssp_rng_ref_put(rng), RES_OK); + + check_memory_allocator(&allocator); + mem_shutdown_proxy_allocator(&allocator); + CHECK(mem_allocated_size(), 0); + return 0; +} +