commit 6f6cfadd6fee0b9a0c0fb0ad3cc4543e73d1e84f
parent f8156070b2c34631b1912763ad4f4498be653d2f
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 3 Oct 2016 15:54:26 +0200
Rearrange the random variate code
Split the distribution in several files
Diffstat:
9 files changed, 609 insertions(+), 506 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -71,9 +71,12 @@ set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(SSP_FILES_SRC
+ ssp_ran.c
+ ssp_ranst_discrete.c
+ ssp_ranst_gaussian.c
+ ssp_ranst_piecewise_linear.c
ssp_rng.c
- ssp_rng_proxy.c
- ssp_distributions.c)
+ ssp_rng_proxy.c)
set(SSP_FILES_INC
ssp_rng_c.h)
set(SSP_FILES_INC_API
diff --git a/src/ssp.h b/src/ssp.h
@@ -506,7 +506,6 @@ ssp_ran_hemisphere_cos_pdf(const float up[3], const float sample[3])
/*******************************************************************************
* Henyey Greenstein distribution
******************************************************************************/
-
/* Return the probability distribution for a Henyey-Greenstein random
* variate with an implicit incident direction in Z */
static INLINE float
@@ -583,7 +582,6 @@ ssp_ran_sphere_hg
/*******************************************************************************
* Gaussian distribution with state
******************************************************************************/
-
SSP_API res_T
ssp_ranst_gaussian_create
(struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */
@@ -592,8 +590,8 @@ ssp_ranst_gaussian_create
SSP_API res_T
ssp_ranst_gaussian_setup
(struct ssp_ranst_gaussian* ran,
- double mu,
- double sigma);
+ const double mu,
+ const double sigma);
SSP_API res_T
ssp_ranst_gaussian_ref_get
@@ -611,12 +609,11 @@ ssp_ranst_gaussian_get
SSP_API double
ssp_ranst_gaussian_pdf
(const struct ssp_ranst_gaussian* ran,
- double x);
+ const double x);
/*******************************************************************************
* Piecewise linear distribution with state
******************************************************************************/
-
SSP_API res_T
ssp_ranst_piecewise_linear_create
(struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */
@@ -627,7 +624,7 @@ ssp_ranst_piecewise_linear_setup
(struct ssp_ranst_piecewise_linear *ran,
const double* intervals,
const double* weights,
- size_t size);
+ const size_t size);
SSP_API res_T
ssp_ranst_piecewise_linear_ref_get
@@ -643,13 +640,16 @@ ssp_ranst_piecewise_linear_get
struct ssp_rng* rng);
/*******************************************************************************
- * Uniform disk distribution
+ * Uniform disk distribution.
+ *
+ * TODO: provide both local and world space version. Return the pdf in addition
+ * of the position. Add the pdf functions that return the pdf of a provided
+ * sample.
******************************************************************************/
-
static FINLINE double*
ssp_ran_uniform_disk
(struct ssp_rng* rng,
- double radius,
+ const double radius,
double pt[2])
{
double theta, r;
diff --git a/src/ssp_distributions.c b/src/ssp_distributions.c
@@ -1,436 +0,0 @@
-/* Copyright (C) |Meso|Star> 2015-2016 (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 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
- * 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 license and that you accept its terms. */
-
-#include <rsys/mem_allocator.h>
-#include <rsys/ref_count.h>
-#include <rsys/dynamic_array_double.h>
-#include <rsys/math.h>
-
-#include "ssp_rng_c.h"
-
-/*******************************************************************************
- * Gaussian distribution utilities
- ******************************************************************************/
-
-using normal_dist = RAN_NAMESPACE::normal_distribution<double>;
-
-struct ssp_ranst_gaussian
-{
- ref_T ref;
- struct mem_allocator* allocator;
- normal_dist* distrib;
- double mu;
- double K1; // 1.0 / sigma;
- double K2; // 1.0 / (sigma * SQRT_2_PI)
-};
-
-static void
-gaussian_release(ref_T* ref)
-{
- ssp_ranst_gaussian* ran;
- ASSERT(ref);
- ran = CONTAINER_OF(ref, ssp_ranst_gaussian, ref);
- if(ran->distrib) {
- ran->distrib->~normal_dist();
- MEM_RM(ran->allocator, ran->distrib);
- }
- MEM_RM(ran->allocator, ran);
-}
-
-/*******************************************************************************
- * Piecewise linear distribution utilities
- ******************************************************************************/
-
-using piecewise_dist = RAN_NAMESPACE::piecewise_linear_distribution<double>;
-
-struct ssp_ranst_piecewise_linear
-{
- ref_T ref;
- struct mem_allocator* allocator;
- piecewise_dist* distrib;
-};
-
-static void
-piecewise_release(ref_T* ref)
-{
- ssp_ranst_piecewise_linear* ran;
- ASSERT(ref);
- ran = CONTAINER_OF(ref, ssp_ranst_piecewise_linear, ref);
- if(ran->distrib) {
- ran->distrib->~piecewise_dist();
- MEM_RM(ran->allocator, ran->distrib);
- }
- MEM_RM(ran->allocator, ran);
-}
-
-/*******************************************************************************
- * Discrete distribution utilities
- ******************************************************************************/
-
-using discrete_dist=RAN_NAMESPACE::discrete_distribution<size_t>;
-
-struct ssp_ranst_discrete {
- struct darray_double pdf;
- ref_T ref;
- struct mem_allocator* allocator;
- discrete_dist* distrib;
-};
-
-static void
-discrete_release(ref_T* ref)
-{
- struct ssp_ranst_discrete* ran;
- ASSERT(ref);
- ran = CONTAINER_OF(ref, struct ssp_ranst_discrete, ref);
- if(ran->distrib) {
- ran->distrib->~discrete_dist();
- MEM_RM(ran->allocator, ran->distrib);
- }
- darray_double_release(&ran->pdf);
- MEM_RM(ran->allocator, ran);
-}
-
-/*******************************************************************************
- * Discrete distribution exported functions
- ******************************************************************************/
-
-res_T
-ssp_ranst_discrete_create
- (struct mem_allocator* allocator,
- struct ssp_ranst_discrete** out_ran)
-{
- ssp_ranst_discrete* ran = nullptr;
- void* mem = nullptr;
- res_T res = RES_OK;
-
- if(!out_ran) {
- res = RES_BAD_ARG;
- goto error;
- }
- allocator = allocator ? allocator : &mem_default_allocator;
-
- ran = static_cast<ssp_ranst_discrete*>(
- MEM_CALLOC(allocator, 1, sizeof(ssp_ranst_discrete)));
- if(!ran) {
- res = RES_MEM_ERR;
- goto error;
- }
- ref_init(&ran->ref);
-
- mem = MEM_ALLOC(allocator, sizeof(discrete_dist));
- if (!mem) {
- res = RES_MEM_ERR;
- goto error;
- }
- ran->allocator = allocator;
- ran->distrib = static_cast<discrete_dist*>(new(mem) discrete_dist);
- darray_double_init(allocator, &ran->pdf);
-
-exit:
- if(out_ran) *out_ran = ran;
- return res;
-error:
- if(ran) {
- SSP(ranst_discrete_ref_put(ran));
- ran = NULL;
- }
- goto exit;
-}
-
-res_T
-ssp_ranst_discrete_setup
- (struct ssp_ranst_discrete* ran,
- const double* weights,
- const size_t nweights)
-{
- double* pdf;
- double sum = 0;
- size_t i;
- res_T res = RES_OK;
-
- if(!ran || !weights || !nweights) {
- res = RES_BAD_ARG;
- goto error;
- }
-
- res = darray_double_resize(&ran->pdf, nweights);
- if(res != RES_OK) goto error;
-
- pdf = darray_double_data_get(&ran->pdf);
-
- for(i = 0; i < nweights; i++) {
- if(weights[i] < 0.f) {
- res = RES_BAD_ARG;
- goto error;
- }
- sum += weights[i];
- pdf[i] = weights[i];
- }
- if (sum == 0) {
- res = RES_BAD_ARG;
- goto error;
- }
-
- sum = 1 / sum;
- for(i = 0; i < nweights; i++) {
- pdf[i] *= sum;
- }
-
- {
- discrete_dist::param_type p{weights, weights + nweights};
- ran->distrib->param(p);
- }
-
-exit:
- return res;
-error:
- if(ran) {
- darray_double_clear(&ran->pdf);
- }
- goto exit;
-}
-
-res_T
-ssp_ranst_discrete_ref_get(struct ssp_ranst_discrete* ran)
-{
- if(!ran) return RES_BAD_ARG;
- ref_get(&ran->ref);
- return RES_OK;
-}
-
-res_T
-ssp_ranst_discrete_ref_put(struct ssp_ranst_discrete* ran)
-{
- if(!ran) return RES_BAD_ARG;
- ref_put(&ran->ref, discrete_release);
- return RES_OK;
-}
-
-size_t
-ssp_ranst_discrete_get(struct ssp_rng* rng, const struct ssp_ranst_discrete* ran)
-{
- class rng_cxx r(*rng);
- return ran->distrib->operator()(r);
-}
-
-double
-ssp_ranst_discrete_pdf(const size_t i, const struct ssp_ranst_discrete* ran)
-{
- ASSERT(ran && i < darray_double_size_get(&ran->pdf));
- return darray_double_cdata_get(&ran->pdf)[i];
-}
-
-/*******************************************************************************
- * Gaussian distribution exported functions
- ******************************************************************************/
-
-res_T
-ssp_ranst_gaussian_create
- (struct mem_allocator* allocator,
- struct ssp_ranst_gaussian** out_ran)
-{
- ssp_ranst_gaussian* ran = nullptr;
- void* mem = nullptr;
- res_T res = RES_OK;
-
- if (!out_ran)
- return RES_BAD_ARG;
-
- allocator = allocator ? allocator : &mem_default_allocator;
-
- ran = static_cast<ssp_ranst_gaussian*>(
- MEM_CALLOC(allocator, 1, sizeof(ssp_ranst_gaussian)));
- if (!ran) {
- res = RES_MEM_ERR;
- goto err;
- }
- ref_init(&ran->ref);
-
- mem = MEM_ALLOC(allocator, sizeof(normal_dist));
- if (!mem) {
- res = RES_MEM_ERR;
- goto err;
- }
- ran->allocator = allocator;
- ran->distrib = static_cast<normal_dist*>(new(mem) normal_dist);
- ran->K1 = -1; /* invalid */
- if (out_ran) *out_ran = ran;
- end:
- return res;
- err:
- if (ran) {
- SSP(ranst_gaussian_ref_put(ran));
- ran = nullptr;
- }
- goto end;
-}
-
-res_T
-ssp_ranst_gaussian_setup
- (struct ssp_ranst_gaussian* ran,
- double mu,
- double sigma)
-{
- if (!ran || sigma < 0)
- return RES_BAD_ARG;
-
- normal_dist::param_type p{mu, sigma};
- ran->distrib->param(p);
- ran->mu = mu;
- ran->K1 = 1.0 / sigma;
- ran->K2 = 1.0 / (sigma * SQRT_2_PI);
-
- return RES_OK;
-}
-
-res_T
-ssp_ranst_gaussian_ref_get
- (struct ssp_ranst_gaussian* ran)
-{
- if(!ran) return RES_BAD_ARG;
- ref_get(&ran->ref);
- return RES_OK;
-}
-
-res_T
-ssp_ranst_gaussian_ref_put
- (struct ssp_ranst_gaussian* ran)
-{
- if(!ran) return RES_BAD_ARG;
- ref_put(&ran->ref, gaussian_release);
- return RES_OK;
-}
-
-double
-ssp_ranst_gaussian_get
- (const struct ssp_ranst_gaussian* ran, struct ssp_rng* rng)
-{
- class rng_cxx r(*rng);
- ASSERT(ran->K1 > 0);
- return ran->distrib->operator()(r);
-}
-
-double
-ssp_ranst_gaussian_pdf
- (const struct ssp_ranst_gaussian* ran, double x)
-{
- const double tmp = (x - ran->mu) * ran->K1;
- ASSERT(ran->K1 > 0);
- return ran->K2 * exp(-0.5 * tmp * tmp);
-}
-
-/*******************************************************************************
- * Piecewise linear distribution exported functions
- ******************************************************************************/
-
-res_T
-ssp_ranst_piecewise_linear_create
- (struct mem_allocator* allocator,
- struct ssp_ranst_piecewise_linear** out_ran)
-{
- ssp_ranst_piecewise_linear* ran = nullptr;
- void* mem = nullptr;
- res_T res = RES_OK;
-
- if (!out_ran)
- return RES_BAD_ARG;
-
- allocator = allocator ? allocator : &mem_default_allocator;
-
- ran = static_cast<ssp_ranst_piecewise_linear*>(
- MEM_CALLOC(allocator, 1, sizeof(ssp_ranst_piecewise_linear)));
- if (!ran) {
- res = RES_MEM_ERR;
- goto err;
- }
- ref_init(&ran->ref);
-
- mem = MEM_ALLOC(allocator, sizeof(piecewise_dist));
- if (!mem) {
- res = RES_MEM_ERR;
- goto err;
- }
- ran->allocator = allocator;
- ran->distrib = static_cast<piecewise_dist*>(new(mem) piecewise_dist);
- if (out_ran) *out_ran = ran;
- end:
- return res;
- err:
- if (ran) {
- SSP(ranst_piecewise_linear_ref_put(ran));
- ran = nullptr;
- }
- goto end;
-}
-
-res_T
-ssp_ranst_piecewise_linear_setup
- (struct ssp_ranst_piecewise_linear* ran,
- const double* intervals,
- const double* weights,
- size_t size)
-{
- if (!ran || !intervals || !weights || size < 2)
- return RES_BAD_ARG;
-
- piecewise_dist::param_type p{intervals, intervals + size, weights};
- ran->distrib->param(p);
-
- return RES_OK;
-}
-
-res_T
-ssp_ranst_piecewise_linear_ref_get
- (struct ssp_ranst_piecewise_linear* ran)
-{
- if(!ran) return RES_BAD_ARG;
- ref_get(&ran->ref);
- return RES_OK;
-}
-
-res_T
-ssp_ranst_piecewise_linear_ref_put
- (struct ssp_ranst_piecewise_linear* ran)
-{
- if(!ran) return RES_BAD_ARG;
- ref_put(&ran->ref, piecewise_release);
- return RES_OK;
-}
-
-double
-ssp_ranst_piecewise_linear_get
- (const struct ssp_ranst_piecewise_linear* ran,
- struct ssp_rng* rng)
-{
- class rng_cxx r(*rng);
- return ran->distrib->operator()(r);
-}
-
diff --git a/src/ssp_ran.c b/src/ssp_ran.c
@@ -0,0 +1,90 @@
+/* Copyright (C) |Meso|Star> 2015-2016 (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 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
+ * 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 license and that you accept its terms. */
+
+#include "ssp_rng_c.h"
+
+/*******************************************************************************
+ * Exported state free random variates
+ ******************************************************************************/
+double
+ssp_ran_exp(struct ssp_rng* rng, const double mu)
+{
+ ASSERT(rng);
+ rng_cxx rng_cxx(*rng);
+ RAN_NAMESPACE::exponential_distribution<double> distribution(mu);
+ return distribution(rng_cxx);
+}
+
+double
+ssp_ran_exp_pdf(const double x, const double mu)
+{
+ ASSERT(x >= 0);
+ return mu * exp(-x * mu);
+}
+
+double
+ssp_ran_gaussian
+ (struct ssp_rng* rng,
+ const double mu,
+ const double sigma)
+{
+ ASSERT(rng);
+ rng_cxx rng_cxx(*rng);
+ RAN_NAMESPACE::normal_distribution<double> distribution(mu, sigma);
+ return distribution(rng_cxx);
+}
+
+double
+ssp_ran_gaussian_pdf(const double x, const double mu, const double sigma)
+{
+ const double tmp = (x - mu) / sigma;
+ return 1.0 / (sigma * SQRT_2_PI) * exp(-0.5 * tmp * tmp);
+}
+
+double
+ssp_ran_lognormal
+ (struct ssp_rng* rng,
+ const double zeta,
+ const double sigma)
+{
+ ASSERT(rng);
+ rng_cxx rng_cxx(*rng);
+ RAN_NAMESPACE::lognormal_distribution<double> distribution(zeta, sigma);
+ return distribution(rng_cxx);
+}
+
+double
+ssp_ran_lognormal_pdf(const double x, const double zeta, const double sigma)
+{
+ const double tmp = log(x) - zeta;
+ return 1.0 / (sigma * x * SQRT_2_PI) * exp(-(tmp*tmp) / (2*sigma*sigma));
+}
+
diff --git a/src/ssp_ranst_discrete.c b/src/ssp_ranst_discrete.c
@@ -0,0 +1,192 @@
+/* Copyright (C) |Meso|Star> 2015-2016 (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 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
+ * 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 license and that you accept its terms. */
+
+#include "ssp_rng_c.h"
+
+#include <rsys/dynamic_array_double.h>
+#include <rsys/mem_allocator.h>
+#include <rsys/ref_count.h>
+
+using discrete_dist=RAN_NAMESPACE::discrete_distribution<size_t>;
+
+struct ssp_ranst_discrete {
+ struct darray_double pdf;
+ ref_T ref;
+ struct mem_allocator* allocator;
+ discrete_dist* distrib;
+};
+
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+static void
+discrete_release(ref_T* ref)
+{
+ struct ssp_ranst_discrete* ran;
+ ASSERT(ref);
+ ran = CONTAINER_OF(ref, struct ssp_ranst_discrete, ref);
+ if(ran->distrib) {
+ ran->distrib->~discrete_dist();
+ MEM_RM(ran->allocator, ran->distrib);
+ }
+ darray_double_release(&ran->pdf);
+ MEM_RM(ran->allocator, ran);
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+ssp_ranst_discrete_create
+ (struct mem_allocator* allocator,
+ struct ssp_ranst_discrete** out_ran)
+{
+ struct ssp_ranst_discrete* ran = nullptr;
+ void* mem = nullptr;
+ res_T res = RES_OK;
+
+ if(!out_ran) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ allocator = allocator ? allocator : &mem_default_allocator;
+
+ ran = static_cast<ssp_ranst_discrete*>
+ (MEM_CALLOC(allocator, 1, sizeof(ssp_ranst_discrete)));
+ if(!ran) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ref_init(&ran->ref);
+
+ mem = MEM_ALLOC(allocator, sizeof(discrete_dist));
+ if(!mem) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ran->allocator = allocator;
+ ran->distrib = static_cast<discrete_dist*>(new(mem) discrete_dist);
+ darray_double_init(allocator, &ran->pdf);
+
+exit:
+ if(out_ran) *out_ran = ran;
+ return res;
+error:
+ if(ran) {
+ SSP(ranst_discrete_ref_put(ran));
+ ran = NULL;
+ }
+ goto exit;
+}
+
+res_T
+ssp_ranst_discrete_setup
+ (struct ssp_ranst_discrete* ran,
+ const double* weights,
+ const size_t nweights)
+{
+ double* pdf;
+ double sum = 0;
+ size_t i;
+ res_T res = RES_OK;
+
+ if(!ran || !weights || !nweights) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ res = darray_double_resize(&ran->pdf, nweights);
+ if(res != RES_OK) goto error;
+
+ pdf = darray_double_data_get(&ran->pdf);
+
+ FOR_EACH(i, 0, nweights) {
+ if(weights[i] < 0.f) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ sum += weights[i];
+ pdf[i] = weights[i];
+ }
+ if(sum == 0) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ sum = 1 / sum;
+ FOR_EACH(i, 0, nweights) {
+ pdf[i] *= sum;
+ }
+
+ {
+ discrete_dist::param_type p{weights, weights + nweights};
+ ran->distrib->param(p);
+ }
+
+exit:
+ return res;
+error:
+ if(ran) {
+ darray_double_clear(&ran->pdf);
+ }
+ goto exit;
+}
+
+res_T
+ssp_ranst_discrete_ref_get(struct ssp_ranst_discrete* ran)
+{
+ if(!ran) return RES_BAD_ARG;
+ ref_get(&ran->ref);
+ return RES_OK;
+}
+
+res_T
+ssp_ranst_discrete_ref_put(struct ssp_ranst_discrete* ran)
+{
+ if(!ran) return RES_BAD_ARG;
+ ref_put(&ran->ref, discrete_release);
+ return RES_OK;
+}
+
+size_t
+ssp_ranst_discrete_get(struct ssp_rng* rng, const struct ssp_ranst_discrete* ran)
+{
+ class rng_cxx r(*rng);
+ return ran->distrib->operator()(r);
+}
+
+double
+ssp_ranst_discrete_pdf(const size_t i, const struct ssp_ranst_discrete* ran)
+{
+ ASSERT(ran && i < darray_double_size_get(&ran->pdf));
+ return darray_double_cdata_get(&ran->pdf)[i];
+}
+
diff --git a/src/ssp_ranst_gaussian.c b/src/ssp_ranst_gaussian.c
@@ -0,0 +1,163 @@
+/* Copyright (C) |Meso|Star> 2015-2016 (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 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
+ * 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 license and that you accept its terms. */
+
+#include "ssp_rng_c.h"
+
+#include <rsys/mem_allocator.h>
+#include <rsys/ref_count.h>
+
+using normal_dist = RAN_NAMESPACE::normal_distribution<double>;
+
+struct ssp_ranst_gaussian {
+ ref_T ref;
+ struct mem_allocator* allocator;
+ normal_dist* distrib;
+ double mu;
+ double K1; /* 1.0 / sigma */
+ double K2; /* 1.0 / (sigma * SQRT_2_PI) */
+};
+
+/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+static void
+gaussian_release(ref_T* ref)
+{
+ ssp_ranst_gaussian* ran;
+ ASSERT(ref);
+ ran = CONTAINER_OF(ref, ssp_ranst_gaussian, ref);
+ if(ran->distrib) {
+ ran->distrib->~normal_dist();
+ MEM_RM(ran->allocator, ran->distrib);
+ }
+ MEM_RM(ran->allocator, ran);
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+ssp_ranst_gaussian_create
+ (struct mem_allocator* allocator,
+ struct ssp_ranst_gaussian** out_ran)
+{
+ ssp_ranst_gaussian* ran = nullptr;
+ void* mem = nullptr;
+ res_T res = RES_OK;
+
+ if(!out_ran)
+ return RES_BAD_ARG;
+
+ allocator = allocator ? allocator : &mem_default_allocator;
+
+ ran = static_cast<ssp_ranst_gaussian*>
+ (MEM_CALLOC(allocator, 1, sizeof(struct ssp_ranst_gaussian)));
+ if(!ran) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ref_init(&ran->ref);
+
+ mem = MEM_ALLOC(allocator, sizeof(normal_dist));
+ if(!mem) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ran->allocator = allocator;
+ ran->distrib = static_cast<normal_dist*>(new(mem) normal_dist);
+ ran->K1 = -1; /* invalid */
+ if(out_ran) *out_ran = ran;
+
+exit:
+ return res;
+error:
+ if(ran) {
+ SSP(ranst_gaussian_ref_put(ran));
+ ran = nullptr;
+ }
+ goto exit;
+}
+
+res_T
+ssp_ranst_gaussian_setup
+ (struct ssp_ranst_gaussian* ran,
+ const double mu,
+ const double sigma)
+{
+ if(!ran || sigma < 0)
+ return RES_BAD_ARG;
+
+ normal_dist::param_type p{mu, sigma};
+ ran->distrib->param(p);
+ ran->mu = mu;
+ ran->K1 = 1.0 / sigma;
+ ran->K2 = 1.0 / (sigma * SQRT_2_PI);
+
+ return RES_OK;
+}
+
+res_T
+ssp_ranst_gaussian_ref_get
+ (struct ssp_ranst_gaussian* ran)
+{
+ if(!ran) return RES_BAD_ARG;
+ ref_get(&ran->ref);
+ return RES_OK;
+}
+
+res_T
+ssp_ranst_gaussian_ref_put
+ (struct ssp_ranst_gaussian* ran)
+{
+ if(!ran) return RES_BAD_ARG;
+ ref_put(&ran->ref, gaussian_release);
+ return RES_OK;
+}
+
+double
+ssp_ranst_gaussian_get
+ (const struct ssp_ranst_gaussian* ran, struct ssp_rng* rng)
+{
+ class rng_cxx r(*rng);
+ ASSERT(ran->K1 > 0);
+ return ran->distrib->operator()(r);
+}
+
+double
+ssp_ranst_gaussian_pdf
+ (const struct ssp_ranst_gaussian* ran, const double x)
+{
+ const double tmp = (x - ran->mu) * ran->K1;
+ ASSERT(ran->K1 > 0);
+ return ran->K2 * exp(-0.5 * tmp * tmp);
+}
+
+
diff --git a/src/ssp_ranst_piecewise_linear.c b/src/ssp_ranst_piecewise_linear.c
@@ -0,0 +1,148 @@
+/* Copyright (C) |Meso|Star> 2015-2016 (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 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
+ * 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 license and that you accept its terms. */
+
+#include "ssp_rng_c.h"
+
+#include <rsys/mem_allocator.h>
+#include <rsys/ref_count.h>
+
+using piecewise_dist = RAN_NAMESPACE::piecewise_linear_distribution<double>;
+
+struct ssp_ranst_piecewise_linear {
+ ref_T ref;
+ struct mem_allocator* allocator;
+ piecewise_dist* distrib;
+};
+
+/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+static void
+piecewise_release(ref_T* ref)
+{
+ ssp_ranst_piecewise_linear* ran;
+ ASSERT(ref);
+ ran = CONTAINER_OF(ref, ssp_ranst_piecewise_linear, ref);
+ if(ran->distrib) {
+ ran->distrib->~piecewise_dist();
+ MEM_RM(ran->allocator, ran->distrib);
+ }
+ MEM_RM(ran->allocator, ran);
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+ssp_ranst_piecewise_linear_create
+ (struct mem_allocator* allocator,
+ struct ssp_ranst_piecewise_linear** out_ran)
+{
+ struct ssp_ranst_piecewise_linear* ran = nullptr;
+ void* mem = nullptr;
+ res_T res = RES_OK;
+
+ if(!out_ran) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ allocator = allocator ? allocator : &mem_default_allocator;
+
+ ran = static_cast<ssp_ranst_piecewise_linear*>
+ (MEM_CALLOC(allocator, 1, sizeof(ssp_ranst_piecewise_linear)));
+ if(!ran) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ref_init(&ran->ref);
+
+ mem = MEM_ALLOC(allocator, sizeof(piecewise_dist));
+ if(!mem) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ran->allocator = allocator;
+ ran->distrib = static_cast<piecewise_dist*>(new(mem) piecewise_dist);
+ if(out_ran) *out_ran = ran;
+
+exit:
+ return res;
+error:
+ if(ran) {
+ SSP(ranst_piecewise_linear_ref_put(ran));
+ ran = nullptr;
+ }
+ goto exit;
+}
+
+res_T
+ssp_ranst_piecewise_linear_setup
+ (struct ssp_ranst_piecewise_linear* ran,
+ const double* intervals,
+ const double* weights,
+ const size_t size)
+{
+ if(!ran || !intervals || !weights || size < 2)
+ return RES_BAD_ARG;
+
+ piecewise_dist::param_type p{intervals, intervals + size, weights};
+ ran->distrib->param(p);
+ return RES_OK;
+}
+
+res_T
+ssp_ranst_piecewise_linear_ref_get
+ (struct ssp_ranst_piecewise_linear* ran)
+{
+ if(!ran) return RES_BAD_ARG;
+ ref_get(&ran->ref);
+ return RES_OK;
+}
+
+res_T
+ssp_ranst_piecewise_linear_ref_put
+ (struct ssp_ranst_piecewise_linear* ran)
+{
+ if(!ran) return RES_BAD_ARG;
+ ref_put(&ran->ref, piecewise_release);
+ return RES_OK;
+}
+
+double
+ssp_ranst_piecewise_linear_get
+ (const struct ssp_ranst_piecewise_linear* ran,
+ struct ssp_rng* rng)
+{
+ class rng_cxx r(*rng);
+ return ran->distrib->operator()(r);
+}
+
diff --git a/src/ssp_rng.c b/src/ssp_rng.c
@@ -639,63 +639,6 @@ ssp_rng_discard(struct ssp_rng* rng, uint64_t n)
return rng->type.discard(rng->state, n);
}
-/*******************************************************************************
- * Exported distributions
- ******************************************************************************/
-double
-ssp_ran_exp(struct ssp_rng* rng, const double mu)
-{
- ASSERT(rng);
- rng_cxx rng_cxx(*rng);
- RAN_NAMESPACE::exponential_distribution<double> distribution(mu);
- return distribution(rng_cxx);
-}
-
-double
-ssp_ran_exp_pdf(const double x, const double mu)
-{
- ASSERT(x >= 0);
- return mu * exp(-x * mu);
-}
-
-double
-ssp_ran_gaussian
- (struct ssp_rng* rng,
- const double mu,
- const double sigma)
-{
- ASSERT(rng);
- rng_cxx rng_cxx(*rng);
- RAN_NAMESPACE::normal_distribution<double> distribution(mu, sigma);
- return distribution(rng_cxx);
-}
-
-double
-ssp_ran_gaussian_pdf(const double x, const double mu, const double sigma)
-{
- const double tmp = (x - mu) / sigma;
- return 1.0 / (sigma * SQRT_2_PI) * exp(-0.5 * tmp * tmp);
-}
-
-double
-ssp_ran_lognormal
- (struct ssp_rng* rng,
- const double zeta,
- const double sigma)
-{
- ASSERT(rng);
- rng_cxx rng_cxx(*rng);
- RAN_NAMESPACE::lognormal_distribution<double> distribution(zeta, sigma);
- return distribution(rng_cxx);
-}
-
-double
-ssp_ran_lognormal_pdf(const double x, const double zeta, const double sigma)
-{
- const double tmp = log(x) - zeta;
- return 1.0 / (sigma * x * SQRT_2_PI) * exp(-(tmp*tmp) / (2*sigma*sigma));
-}
-
#ifdef COMPILER_CL
#pragma warning(pop)
#endif
diff --git a/src/ssp_rng_c.h b/src/ssp_rng_c.h
@@ -74,7 +74,7 @@ public:
typedef uint64_t result_type;
private:
- ssp_rng& rng;
+ struct ssp_rng& rng;
};
#endif /* SSP_RNG_C_H */
