star-sp

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

commit b6b1eb7b36edb65cd7e290680baf77ea737b0029
parent 6d71a3bac06adafda9f8949d53a50512cfac1439
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed, 17 Jan 2018 10:08:41 +0100

Add and test the uniform circle random variate

Diffstat:
Mcmake/CMakeLists.txt | 1+
Msrc/ssp.h | 29++++++++++++++++++++++++++++-
Msrc/ssp_ran.c | 30++++++++++++++++++++++++++++++
Asrc/test_ssp_ran_circle.c | 42++++++++++++++++++++++++++++++++++++++++++
Asrc/test_ssp_ran_circle.h | 101+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
5 files changed, 202 insertions(+), 1 deletion(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -164,6 +164,7 @@ if(NOT NO_TEST) new_test(test_ssp_ran_discrete) new_test(test_ssp_ran_hemisphere ${MATH_LIB}) new_test(test_ssp_ran_sphere ${MATH_LIB}) + new_test(test_ssp_ran_circle ${MATH_LIB}) new_test(test_ssp_ran_hg ${MATH_LIB}) new_test(test_ssp_ran_triangle ${MATH_LIB}) new_test(test_ssp_rng_proxy) diff --git a/src/ssp.h b/src/ssp.h @@ -387,6 +387,33 @@ ssp_ran_sphere_uniform_float_pdf(void) } /******************************************************************************* + * Circle distribution + ******************************************************************************/ +SSP_API double* +ssp_ran_circle_uniform + (struct ssp_rng* rng, + double sample[2], + double* pdf); /* Can be NULL => pdf not computed */ + +SSP_API float* +ssp_ran_circle_uniform_float + (struct ssp_rng* rng, + float sample[2], + float* pdf); /* Can be NULL => pdf not computed */ + +static INLINE double +ssp_ran_circle_uniform_pdf(void) +{ + return 1/(2*PI); +} + +static INLINE float +ssp_ran_circle_uniform_float_pdf(void) +{ + return 1/(2*(float)PI); +} + +/******************************************************************************* * Triangle distribution ******************************************************************************/ /* Uniform sampling of a triangle */ @@ -842,7 +869,7 @@ ssp_ranst_piecewise_linear_float_pdf float x); /******************************************************************************* - * Uniform disk distribution. + * Uniform distribution of a point into a disk. ******************************************************************************/ SSP_API double* ssp_ran_uniform_disk_local diff --git a/src/ssp_ran.c b/src/ssp_ran.c @@ -205,6 +205,36 @@ ssp_ran_sphere_uniform_float } double* +ssp_ran_circle_uniform + (struct ssp_rng* rng, + double sample[2], + double* pdf) +{ + double theta; + ASSERT(rng && sample); + theta = ssp_rng_uniform_double(rng, 0, 2*PI); + sample[0] = cos(theta); + sample[1] = sin(theta); + if(pdf) *pdf = 1/(2*PI); + return sample; +} + +float* +ssp_ran_circle_uniform_float + (struct ssp_rng* rng, + float sample[2], + float* pdf) +{ + float theta; + ASSERT(rng && sample); + theta = ssp_rng_uniform_float(rng, 0, 2*(float)PI); + sample[0] = cosf(theta); + sample[1] = sinf(theta); + if(pdf) *pdf = 1/(2*(float)PI); + return sample; +} + +double* ssp_ran_triangle_uniform (struct ssp_rng* rng, const double v0[3], diff --git a/src/test_ssp_ran_circle.c b/src/test_ssp_ran_circle.c @@ -0,0 +1,42 @@ +/* Copyright (C) |Meso|Star> 2015-2017 (contact@meso-star.com) + * + * 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. */ + +#define TYPE_FLOAT 0 +#include "test_ssp_ran_circle.h" + +#define TYPE_FLOAT 1 +#include "test_ssp_ran_circle.h" + +int +main(int argc, char** argv) +{ + (void)argc, (void)argv; + test_float(); + test_double(); + return 0; +} diff --git a/src/test_ssp_ran_circle.h b/src/test_ssp_ran_circle.h @@ -0,0 +1,101 @@ +/* Copyright (C) |Meso|Star> 2015-2017 (contact@meso-star.com) + * + * 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. */ + +#ifndef TEST_SSP_RAN_CIRCLE_H +#define TEST_SSP_RAN_CIRCLE_H + +#include "ssp.h" +#include "test_ssp_utils.h" + +#define NSAMPS 128 + +#endif /* TEST_SSP_RAN_CIRCLE_H */ + +#if TYPE_FLOAT == 0 + #include <rsys/double2.h> + + #define REAL double + #define TEST test_double + #define RAN_CIRCLE_UNIFORM ssp_ran_circle_uniform + #define RAN_CIRCLE_UNIFORM_PDF ssp_ran_circle_uniform_pdf + #define EQ_EPS eq_eps + #define R2_EQ_EPS d2_eq_eps + #define R2_IS_NORMALIZED d2_is_normalized + +#elif TYPE_FLOAT == 1 + #include <rsys/float2.h> + + #define REAL float + #define TEST test_float + #define RAN_CIRCLE_UNIFORM ssp_ran_circle_uniform_float + #define RAN_CIRCLE_UNIFORM_PDF ssp_ran_circle_uniform_float_pdf + #define EQ_EPS eq_epsf + #define R2_EQ_EPS f2_eq_eps + #define R2_IS_NORMALIZED f2_is_normalized +#else + #error "TYPE_FLOAT must be defined either 0 or 1" +#endif + +static void +TEST(void) +{ + struct ssp_rng* rng; + struct mem_allocator allocator; + REAL samps[NSAMPS][4]; + int i = 0, j = 9; + + CHK(mem_init_proxy_allocator(&allocator, &mem_default_allocator) == RES_OK); + CHK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng) == RES_OK); + CHK(EQ_EPS(RAN_CIRCLE_UNIFORM_PDF(), 1/(2*(REAL)PI), (REAL)1.e-6)); + + FOR_EACH(i, 0, NSAMPS) { + CHK(RAN_CIRCLE_UNIFORM(rng, samps[i], &samps[i][3]) == samps[i]); + CHK(R2_IS_NORMALIZED(samps[i])); + CHK(EQ_EPS(samps[i][3], 1/(2*(REAL)PI), (REAL)1.e-6)); + FOR_EACH(j, 0, i) { + CHK(!R2_EQ_EPS(samps[j], samps[i], (REAL)1.e-6)); + } + printf("%g %g\n", SPLIT2(samps[i])); + } + + ssp_rng_ref_put(rng); + check_memory_allocator(&allocator); + mem_shutdown_proxy_allocator(&allocator); + + CHK(mem_allocated_size() == 0); +} + +#undef REAL +#undef TEST +#undef RAN_CIRCLE_UNIFORM +#undef RAN_CIRCLE_UNIFORM_PDF +#undef EQ_EPS +#undef R2_EQ_EPS +#undef R2_IS_NORMALIZED +#undef TYPE_FLOAT +