star-sp

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

commit 1506e802be885d6662140854166579fee88d89e2
parent cdfd1506dbf8930664b6e924652f5860f60b8fa0
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue, 14 Apr 2015 19:33:03 +0200

Add and test the proxy RNG API

Diffstat:
Mcmake/CMakeLists.txt | 5+++--
Msrc/ssp.h | 28+++++++++++++++++++++++++++-
Msrc/ssp_rng.c | 10+---------
Asrc/ssp_rng_c.h | 46++++++++++++++++++++++++++++++++++++++++++++++
Asrc/ssp_rng_proxy.c | 445+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/test_ssp_rng_proxy.c | 92+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6 files changed, 614 insertions(+), 12 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -54,8 +54,8 @@ set(VERSION_MINOR 0) set(VERSION_PATCH 0) set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH}) -set(SSP_FILES_SRC ssp_rng.c) -set(SSP_FILES_INC ) +set(SSP_FILES_SRC ssp_rng.c ssp_rng_proxy.c) +set(SSP_FILES_INC ssp_rng_c.h) set(SSP_FILES_INC_API ssp.h) # Prepend each file in the `SSP_FILES_<SRC|INC>' list by `SSP_SOURCE_DIR' @@ -104,6 +104,7 @@ if(NOT NO_TEST) add_test(test_ssp_rng_ranlunx48 test_ssp_rng ranlux48) new_test(test_ssp_ran_hemisphere m) + new_test(test_ssp_rng_proxy) endif() ################################################################################ diff --git a/src/ssp.h b/src/ssp.h @@ -56,6 +56,7 @@ /* Forward declaration of opaque types */ struct ssp_rng; +struct ssp_rng_proxy; /* Forward declaration of external types */ struct mem_allocator; @@ -90,7 +91,7 @@ SSP_API const struct ssp_rng_type ssp_rng_mt19937_64; SSP_API const struct ssp_rng_type ssp_rng_ranlux48; /******************************************************************************* - * API declaration + * Random Number Generator API ******************************************************************************/ SSP_API res_T ssp_rng_create @@ -154,6 +155,31 @@ ssp_rng_read FILE* stream); /******************************************************************************* + * Proxy of Random Number Generators. A proxy ensures that the RNG that it + * manages, provide unconnected random numbers + ******************************************************************************/ +SSP_API res_T +ssp_rng_proxy_create + (struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */ + const struct ssp_rng_type* type, + const size_t nbuckets, + struct ssp_rng_proxy** proxy); + +SSP_API res_T +ssp_rng_proxy_ref_get + (struct ssp_rng_proxy* proxy); + +SSP_API res_T +ssp_rng_proxy_ref_put + (struct ssp_rng_proxy* proxy); + +SSP_API res_T +ssp_rng_proxy_create_rng + (struct ssp_rng_proxy* proxy, + const size_t ibucket, + struct ssp_rng** rng); + +/******************************************************************************* * Hemisphere distribution ******************************************************************************/ /* Uniform sampling of an unit hemisphere whose up direction is implicitly the diff --git a/src/ssp_rng.c b/src/ssp_rng.c @@ -29,21 +29,13 @@ * 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 "ssp_rng_c.h" #include <rsys/mem_allocator.h> -#include <rsys/ref_count.h> #include <random> #include <sstream> -struct ssp_rng { - struct ssp_rng_type type; - void* state; - struct mem_allocator* allocator; - ref_T ref; -}; - /******************************************************************************* * KISS PRNG ******************************************************************************/ diff --git a/src/ssp_rng_c.h b/src/ssp_rng_c.h @@ -0,0 +1,46 @@ +/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com) + * + * This software is a library whose purpose is to generate [pseudo] random + * numbers and random variates. + * + * 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. */ + +#ifndef SSP_RNG_C_H +#define SSP_RNG_C_H + +#include "ssp.h" +#include <rsys/ref_count.h> + +struct ssp_rng { + struct ssp_rng_type type; + void* state; + struct mem_allocator* allocator; + ref_T ref; +}; + +#endif /* SSP_RNG_C_H */ + diff --git a/src/ssp_rng_proxy.c b/src/ssp_rng_proxy.c @@ -0,0 +1,445 @@ +/* Copyright (C) |Meso|Star> 2015 (contact@meso-star.com) + * + * This software is a library whose purpose is to generate [pseudo] random + * numbers and random variates. + * + * 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_rng_c.h" + +#include <rsys/mutex.h> +#include <rsys/ref_count.h> +#include <rsys/stretchy_array.h> + +#include <limits.h> + +#define BUCKET_SIZE 1000000 + +struct rng_state_cache { + FILE* stream; + long head, tail; /* Position onto the head/tail `states' stream */ +}; + +struct ssp_rng_proxy { + struct ssp_rng_type type; /* Wrapped RNG type */ + + struct ssp_rng* rng; /* Main RNG */ + + int* buckets; /* Flag that defines which buckets are created */ + struct ssp_rng** pools; /* Cache of bucketized RNGs */ + struct rng_state_cache* states; + + struct mutex* mutex; + struct mem_allocator* allocator; + ref_T ref; +}; + +static struct ssp_rng* +rng_proxy_next_ran_pool + (struct ssp_rng_proxy* proxy, + const size_t bucket_id); + +/******************************************************************************* + * State cache + ******************************************************************************/ +static res_T +rng_state_cache_init(struct rng_state_cache* cache) +{ + ASSERT(cache); + cache->stream = tmpfile(); + if(!cache->stream) return RES_IO_ERR; + cache->head = cache->tail = ftell(cache->stream); + return RES_OK; +} + +static void +rng_state_cache_release(struct rng_state_cache* cache) +{ + ASSERT(cache); + if(cache->stream) fclose(cache->stream); +} + +static char +rng_state_cache_is_empty(struct rng_state_cache* cache) +{ + ASSERT(cache); + return cache->head == cache->tail; +} + +static res_T +rng_state_cache_read(struct rng_state_cache* cache, struct ssp_rng* rng) +{ + res_T res; + ASSERT(cache && rng && !rng_state_cache_is_empty(cache)); + + /* Read the rng state from the stream */ + fseek(cache->stream, cache->head, SEEK_SET); + res = ssp_rng_read(rng, cache->stream); + if(res != RES_OK) return res; + cache->head = ftell(cache->stream); + + /* Flush the state cache stream if ther is no more RNG state */ + if(rng_state_cache_is_empty(cache)) { + fclose(cache->stream); + cache->stream = tmpfile(); + if(!cache->stream) return RES_IO_ERR; + cache->head = cache->tail = ftell(cache->stream); + } + + return RES_OK; +} + +static res_T +rng_state_cache_write(struct rng_state_cache* cache, struct ssp_rng* rng) +{ + res_T res; + ASSERT(cache && rng); + fseek(cache->stream, cache->tail, SEEK_SET); + res = ssp_rng_write(rng, cache->stream); + if(res != RES_OK) return res; + cache->tail = ftell(cache->stream); + return RES_OK; +} + +/******************************************************************************* + * Bucketized RNG + ******************************************************************************/ +struct rng_bucket { + struct ssp_rng* pool; /* Wrapped RNG providing a pool of BUCKET_SIZE RNs */ + struct ssp_rng_proxy* proxy; /* The RNG proxy */ + size_t name; /* Unique bucket identifier */ + size_t count; /* Remaining unique random numbers in `pool' */ +}; + +static INLINE void +rng_bucket_next_ran_pool(struct rng_bucket* rng) +{ + ASSERT(rng); + rng->pool = rng_proxy_next_ran_pool(rng->proxy, rng->name); + rng->count = BUCKET_SIZE; +} + +static void +rng_bucket_set(void* data, const unsigned long seed) +{ + (void)data, (void)seed; + FATAL("The seed can't be set explicitly on proxy managed RNGs\n"); +} + +static unsigned long +rng_bucket_get(void* data) +{ + struct rng_bucket* rng = (struct rng_bucket*)data; + ASSERT(data); + if(!rng->count) rng_bucket_next_ran_pool(rng); + --rng->count; + return ssp_rng_get(rng->pool); +} + +static unsigned long +rng_bucket_uniform_int + (void* data, const unsigned long lower, const unsigned long upper) +{ + struct rng_bucket* rng = (struct rng_bucket*)data; + ASSERT(data); + if(!rng->count) rng_bucket_next_ran_pool(rng); + --rng->count; + return ssp_rng_uniform_int(rng->pool, lower, upper); +} + +static double +rng_bucket_uniform_double(void* data, const double lower, const double upper) +{ + struct rng_bucket* rng = (struct rng_bucket*)data; + ASSERT(data); + if(!rng->count) rng_bucket_next_ran_pool(rng); + --rng->count; + return ssp_rng_uniform_double(rng->pool, lower, upper); +} + +static double +rng_bucket_canonical(void* data) +{ + struct rng_bucket* rng = (struct rng_bucket*)data; + ASSERT(data); + if(!rng->count) rng_bucket_next_ran_pool(rng); + --rng->count; + return ssp_rng_canonical(rng->pool); +} + +static res_T +rng_bucket_read(void* data, FILE* file) +{ + (void)data, (void)file; + FATAL("A proxy managed RNG can't be de-serialized\n"); + return RES_BAD_OP; +} + +static res_T +rng_bucket_write(const void* data, FILE* file) +{ + (void)data, (void)file; + FATAL("A proxy managed RNG can't be serialized\n"); + return RES_BAD_OP; +} + +static void +rng_bucket_init(struct mem_allocator* allocator, void* data) +{ + struct rng_bucket* rng = (struct rng_bucket*)data; + ASSERT(data); + (void)allocator; + rng->proxy = NULL; + rng->pool = NULL; + rng->name = SIZE_MAX; + rng->count = 0; +} + +static void +rng_bucket_release(void* data) +{ + struct rng_bucket* rng = (struct rng_bucket*)data; + ASSERT(data); + ATOMIC_SET(&rng->proxy->buckets[rng->name], 0); + SSP(rng_proxy_ref_put(rng->proxy)); +} + +static const struct ssp_rng_type rng_bucket = { + rng_bucket_init, + rng_bucket_release, + rng_bucket_set, + rng_bucket_get, + rng_bucket_uniform_int, + rng_bucket_uniform_double, + rng_bucket_canonical, + rng_bucket_read, + rng_bucket_write, + 0, /* Dummy value */ + ULONG_MAX, /* Dummy value */ + sizeof(struct rng_bucket) +}; + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +struct ssp_rng* +rng_proxy_next_ran_pool + (struct ssp_rng_proxy* proxy, + const size_t bucket_name) +{ + res_T res = RES_OK; + ASSERT(proxy && bucket_name <= sa_size(proxy->buckets)); + + mutex_lock(proxy->mutex); + if(rng_state_cache_is_empty(proxy->states + bucket_name)) { + size_t ibucket; + /* Register a new state for *all* buckets */ + FOR_EACH(ibucket, 0, sa_size(proxy->states)) { + size_t irand; + res = rng_state_cache_write(proxy->states + ibucket, proxy->rng); + if(res != RES_OK) FATAL("RNG proxy: Unrecoverable IO error\n"); + FOR_EACH(irand, 0, BUCKET_SIZE) ssp_rng_get(proxy->rng); /* Discard */ + } + } + mutex_unlock(proxy->mutex); + + /* Set the new RNG pool state */ + res = rng_state_cache_read + (proxy->states + bucket_name, proxy->pools[bucket_name]); + if(res != RES_OK) FATAL("RNG proxy: Unrecoverable IO error\n"); + + return proxy->pools[bucket_name]; +} + +void +rng_proxy_clear(struct ssp_rng_proxy* proxy) +{ + size_t ibucket; + ASSERT(proxy); + ASSERT(sa_size(proxy->pools) == sa_size(proxy->buckets)); + ASSERT(sa_size(proxy->pools) == sa_size(proxy->states)); + + FOR_EACH(ibucket, 0, sa_size(proxy->pools)) { + ASSERT(proxy->buckets[ibucket] == 0); /* No bucket RNG should be created */ + SSP(rng_ref_put(proxy->pools[ibucket])); + rng_state_cache_release(proxy->states + ibucket); + } + sa_clear(proxy->buckets); + sa_clear(proxy->pools); + sa_clear(proxy->states); +} + +static res_T +rng_proxy_setup(struct ssp_rng_proxy* proxy, const size_t nbuckets) +{ + size_t ibucket; + res_T res = RES_OK; + + ASSERT(proxy && nbuckets); + rng_proxy_clear(proxy); + + sa_add(proxy->states, nbuckets); + sa_add(proxy->pools, nbuckets); + sa_add(proxy->buckets, nbuckets); + + FOR_EACH(ibucket, 0, nbuckets) { + res = rng_state_cache_init(proxy->states+ibucket); + if(res != RES_OK) goto error; + + res = ssp_rng_create(proxy->allocator, &proxy->type, proxy->pools+ibucket); + if(res != RES_OK) goto error; + + proxy->buckets[ibucket] = 0; + } + +exit: + return res; +error: + if(proxy) rng_proxy_clear(proxy); + goto exit; +} + +void +rng_proxy_release(ref_T* ref) +{ + struct ssp_rng_proxy* proxy; + ASSERT(ref); + proxy = CONTAINER_OF(ref, struct ssp_rng_proxy, ref); + rng_proxy_clear(proxy); + sa_release(proxy->states); + sa_release(proxy->pools); + sa_release(proxy->buckets); + SSP(rng_ref_put(proxy->rng)); + if(proxy->mutex) mutex_destroy(proxy->mutex); + MEM_FREE(proxy->allocator, proxy); +} + +/******************************************************************************* + * Exported functions + ******************************************************************************/ +res_T +ssp_rng_proxy_create + (struct mem_allocator* mem_allocator, + const struct ssp_rng_type* type, + const size_t nbuckets, + struct ssp_rng_proxy** out_proxy) +{ + struct mem_allocator* allocator = NULL; + struct ssp_rng_proxy* proxy = NULL; + res_T res = RES_OK; + + if(!type || !out_proxy || !nbuckets) { + res = RES_BAD_ARG; + goto error; + } + + allocator = mem_allocator ? mem_allocator : &mem_default_allocator; + proxy = (struct ssp_rng_proxy*) + MEM_CALLOC(allocator, 1, sizeof(struct ssp_rng_proxy)); + if(!proxy) { + res = RES_MEM_ERR; + goto error; + } + proxy->allocator = allocator; + ref_init(&proxy->ref); + + res = ssp_rng_create(allocator, type, &proxy->rng); + if(res != RES_OK) goto error; + proxy->type = *type; + + proxy->mutex = mutex_create(); + if(!proxy->mutex) { + res = RES_MEM_ERR; + goto error; + } + + res = rng_proxy_setup(proxy, nbuckets); + if(res != RES_OK) goto error; + +exit: + if(out_proxy) *out_proxy = proxy; + return res; +error: + if(proxy) { + SSP(rng_proxy_ref_put(proxy)); + proxy = NULL; + } + goto exit; +} + +res_T +ssp_rng_proxy_ref_get(struct ssp_rng_proxy* proxy) +{ + if(!proxy) return RES_BAD_ARG; + ref_get(&proxy->ref); + return RES_OK; +} + +res_T +ssp_rng_proxy_ref_put(struct ssp_rng_proxy* proxy) +{ + if(!proxy) return RES_BAD_ARG; + ref_put(&proxy->ref, rng_proxy_release); + return RES_OK; +} + +res_T +ssp_rng_proxy_create_rng + (struct ssp_rng_proxy* proxy, + const size_t ibucket, + struct ssp_rng** out_rng) +{ + struct ssp_rng* rng = NULL; + res_T res = RES_OK; + + if(!proxy || ibucket >= sa_size(proxy->buckets) || !out_rng) { + res = RES_BAD_ARG; + goto error; + } + + if(ATOMIC_CAS(&proxy->buckets[ibucket], 1, 0) == 1) { + res = RES_BAD_ARG; + goto error; + } + + res = ssp_rng_create(proxy->allocator, &rng_bucket, &rng); + if(res != RES_OK) goto error; + ((struct rng_bucket*)rng->state)->name = ibucket; + ((struct rng_bucket*)rng->state)->proxy = proxy; + SSP(rng_proxy_ref_get(proxy)); + +exit: + if(out_rng) *out_rng = rng; + return res; +error: + if(rng) { + SSP(rng_ref_put(rng)); + rng = NULL; + } + goto exit; +} + diff --git a/src/test_ssp_rng_proxy.c b/src/test_ssp_rng_proxy.c @@ -0,0 +1,92 @@ +/* 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" + +int +main(int argc, char** argv) +{ + struct ssp_rng_proxy* proxy; + struct ssp_rng* rng[4]; + struct mem_allocator allocator; + size_t i = 0; + (void)argc, (void)argv; + + mem_init_proxy_allocator(&allocator, &mem_default_allocator); + + CHECK(ssp_rng_proxy_create(NULL, NULL, 0, NULL), RES_BAD_ARG); + CHECK(ssp_rng_proxy_create(NULL, &ssp_rng_mt19937_64, 0, NULL), RES_BAD_ARG); + CHECK(ssp_rng_proxy_create(NULL, NULL, 4, NULL), RES_BAD_ARG); + CHECK(ssp_rng_proxy_create(NULL, &ssp_rng_mt19937_64, 4, NULL), RES_BAD_ARG); + CHECK(ssp_rng_proxy_create(NULL, NULL, 0, &proxy), RES_BAD_ARG); + CHECK(ssp_rng_proxy_create(NULL, &ssp_rng_mt19937_64, 0, &proxy), RES_BAD_ARG); + CHECK(ssp_rng_proxy_create(NULL, NULL, 4, &proxy), RES_BAD_ARG); + CHECK(ssp_rng_proxy_create(NULL, &ssp_rng_mt19937_64, 4, &proxy), RES_OK); + + + CHECK(ssp_rng_proxy_ref_get(NULL), RES_BAD_ARG); + CHECK(ssp_rng_proxy_ref_get(proxy), RES_OK); + CHECK(ssp_rng_proxy_ref_put(NULL), RES_BAD_ARG); + CHECK(ssp_rng_proxy_ref_put(proxy), RES_OK); + CHECK(ssp_rng_proxy_ref_put(proxy), RES_OK); + + CHECK(ssp_rng_proxy_create(&allocator, &ssp_rng_mt19937_64, 4, &proxy), RES_OK); + + CHECK(ssp_rng_proxy_create_rng(NULL, 0, NULL), RES_BAD_ARG); + CHECK(ssp_rng_proxy_create_rng(proxy, 0, NULL), RES_BAD_ARG); + CHECK(ssp_rng_proxy_create_rng(NULL, 0, &rng[0]), RES_BAD_ARG); + CHECK(ssp_rng_proxy_create_rng(proxy, 0, &rng[0]), RES_OK); + CHECK(ssp_rng_proxy_create_rng(proxy, 0, &rng[1]), RES_BAD_ARG); + CHECK(ssp_rng_ref_put(rng[0]), RES_OK); + + CHECK(ssp_rng_proxy_create_rng(proxy, 4, &rng[0]), RES_BAD_ARG); + FOR_EACH(i, 0, 4) { + CHECK(ssp_rng_proxy_create_rng(proxy, i, &rng[i]), RES_OK); + } + CHECK(ssp_rng_proxy_ref_put(proxy), RES_OK); + + FOR_EACH(i, 0, 2000000) { + unsigned long r[4]; + int j; + FOR_EACH(j, 0, 4) { + int k; + r[j] = ssp_rng_get(rng[j]); + FOR_EACH(k, 0, j) + NCHECK(r[k], r[j]); + } + } + + FOR_EACH(i, 0, 4) CHECK(ssp_rng_ref_put(rng[i]), RES_OK); + + check_memory_allocator(&allocator); + mem_shutdown_proxy_allocator(&allocator); + return 0; +}