commit e8e67b3fe8aa1d4b0e13f087b51d52e7b762203c
parent 4da2e6bed1087fa334d8d11e41f2057c07cff459
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 18 Jun 2015 23:05:58 +0200
Add and test the uniform sphere random variate
Diffstat:
5 files changed, 120 insertions(+), 30 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -141,6 +141,7 @@ if(NOT NO_TEST)
register_test(test_ssp_rng_aes test_ssp_rng aes)
endif()
new_test(test_ssp_ran_hemisphere ${MATH_LIB})
+ new_test(test_ssp_ran_sphere ${MATH_LIB})
new_test(test_ssp_rng_proxy)
endif()
diff --git a/src/ssp.h b/src/ssp.h
@@ -89,15 +89,12 @@ SSP_API const struct ssp_rng_type ssp_rng_kiss;
SSP_API const struct ssp_rng_type ssp_rng_mt19937_64;
/* 48-bits RANLUX builtin PRNG type of Lusher and James, 1994 */
SSP_API const struct ssp_rng_type ssp_rng_ranlux48;
-/*
-random_device generator from Boost / C++11
-std::random_device is a uniformly-distributed integer random number generator
-that produces non-deterministic random numbers.
-std::random_device may be implemented in terms of an implementation-defined
-pseudo-random number engine if a non-deterministic source (e.g. a hardware device)
-is not available to the implementation.
-In this case each std::random_device object may generate the same number sequence.
-*/
+/* A random_device RNG is a uniformly-distributed integer random number generator
+ * that produces non-deterministic random numbers. It may may be implemented in
+ * terms of an implementation-defined pseudo-random number engine if a
+ * non-deterministic source (e.g. a hardware device) is not available to the
+ * implementation. In this case each random_device object may generate the same
+ * number sequence. */
SSP_API const struct ssp_rng_type ssp_rng_random_device;
/* Counter Based RNG threefry of Salmon, Moraes, Dror & Shaw */
SSP_API const struct ssp_rng_type ssp_rng_threefry;
@@ -133,14 +130,14 @@ SSP_API uint64_t
ssp_rng_get
(struct ssp_rng* rng);
-/* Uniform random integer distribution in [lower, upper] */
+/* Uniform random integer distribution in [lower, upper) */
SSP_API uint64_t
ssp_rng_uniform_uint64
(struct ssp_rng* rng,
const uint64_t lower,
const uint64_t upper);
-/* Uniform random floating point distribution in [lower, upper] */
+/* Uniform random floating point distribution in [lower, upper) */
SSP_API double
ssp_rng_uniform_double
(struct ssp_rng* rng,
@@ -172,7 +169,7 @@ ssp_rng_read
SSP_API double
ssp_rng_entropy
-(const struct ssp_rng* rng);
+ (const struct ssp_rng* rng);
/*******************************************************************************
* Proxy of Random Number Generators - It manages a list of `nbuckets' RNGs
@@ -202,6 +199,34 @@ ssp_rng_proxy_create_rng
struct ssp_rng** rng);
/*******************************************************************************
+ * Sphere distribution
+ ******************************************************************************/
+/* Uniform sampling of an unit sphere. The PDF of the generated sample is
+ * stored in sample[3] */
+static INLINE float*
+ssp_ran_sphere_uniform(struct ssp_rng* rng, float sample[4])
+{
+ double phi, cos_theta, sin_theta, v;
+ ASSERT(rng && sample);
+ phi = ssp_rng_uniform_double(rng, 0.0, 2.0*PI);
+ v = ssp_rng_canonical(rng);
+ cos_theta = 1.0 - 2.0*v;
+ sin_theta = 2.0 * sqrt(v * (1.0 - v));
+ sample[0] = (float)(cos(phi) * sin_theta);
+ sample[1] = (float)(sin(phi) * sin_theta);
+ sample[2] = (float)cos_theta;
+ sample[3] = (float)(1.0/(4.0*PI));
+ return sample;
+}
+
+/* Return the probability distribution for a sphere uniform random variate */
+static INLINE float
+ssp_ran_sphere_pdf(void)
+{
+ return (float)(1.0/(4.0*PI));
+}
+
+/*******************************************************************************
* 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
@@ -48,26 +48,26 @@
#ifdef COMPILER_GCC
#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wconversion" /* conversion to ‘yyy’ from ‘xxx’ may alter its value */
+ #pragma GCC diagnostic ignored "-Wconversion" /* unsafe conversion */
#elif defined(COMPILER_CL)
-// disable some warnings in Random123 includes
-#pragma warning(push)
-#pragma warning(disable:4100) /* unreferenced formal parameter */
-#pragma warning(disable:4127) /* conditional expression is constant */
-#pragma warning(disable:4512) /* assignment operator could not be generated */
-#pragma warning(disable:4521) /* multiple copy constructors specified */
+ /* disable some warnings in Random123 includes */
+ #pragma warning(push)
+ #pragma warning(disable:4100) /* unreferenced formal parameter */
+ #pragma warning(disable:4127) /* conditional expression is constant */
+ #pragma warning(disable:4512) /* assignment operator could not be generated */
+ #pragma warning(disable:4521) /* multiple copy constructors specified */
#endif
-#include "Random123/conventional/Engine.hpp"
-#include "Random123/threefry.h"
+#include <Random123/conventional/Engine.hpp>
+#include <Random123/threefry.h>
#ifdef WITH_R123_AES
-#include "Random123/aes.h"
+ #include <Random123/aes.h>
#endif
#ifdef COMPILER_GCC
-#pragma GCC diagnostic pop
+ #pragma GCC diagnostic pop
#elif defined(COMPILER_CL)
-#pragma warning(pop)
+ #pragma warning(pop)
#endif
#include <sstream>
@@ -81,6 +81,7 @@ struct rng_kiss { uint32_t x, y, z, c; };
static res_T
rng_kiss_set(void* data, const size_t seed)
+
{
struct rng_kiss* kiss = (struct rng_kiss*)data;
RAN_NAMESPACE::mt19937 rng_mt(seed % UINT32_MAX);
@@ -332,7 +333,8 @@ template<>
double
rng_cxx_entropy<RAN_NAMESPACE::random_device>(const void* data)
{
- const RAN_NAMESPACE::random_device* rng = (const RAN_NAMESPACE::random_device*)data;
+ const RAN_NAMESPACE::random_device* rng =
+ (const RAN_NAMESPACE::random_device*)data;
ASSERT(rng);
return rng->entropy();
}
@@ -389,9 +391,8 @@ const struct ssp_rng_type ssp_rng_random_device = {
};
/*******************************************************************************
-* Random123 Counter Based RNG
-******************************************************************************/
-
+ * Random123 Counter Based RNG
+ ******************************************************************************/
typedef r123::Engine<r123::Threefry4x64> threefry_t;
#ifdef WITH_R123_AES
diff --git a/src/test_ssp_ran_hemisphere.c b/src/test_ssp_ran_hemisphere.c
@@ -1,7 +1,5 @@
/* 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 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
diff --git a/src/test_ssp_ran_sphere.c b/src/test_ssp_ran_sphere.c
@@ -0,0 +1,65 @@
+/* Copyright (C) |Meso|Star> 2015 (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. */
+
+#include "ssp.h"
+#include "test_ssp_utils.h"
+
+#define NSAMPS 128
+
+int
+main(int argc, char** argv)
+{
+ struct ssp_rng* rng;
+ struct mem_allocator allocator;
+ float samps[NSAMPS][4];
+ float* f = NULL;
+ int i = 0, j = 0;
+ (void)argc, (void)argv;
+
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+
+ CHECK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng), RES_OK);
+
+ FOR_EACH(i, 0, NSAMPS) {
+ f = ssp_ran_sphere_uniform(rng, samps[i]);
+ CHECK(f, samps[i]);
+ CHECK(f3_is_normalized(f), 1);
+ CHECK(eq_epsf(samps[i][3], (float)(1.0/(4.0*PI)), 1.e-6f), 1);
+ CHECK(eq_epsf(samps[i][3], ssp_ran_sphere_pdf(), 1.e-6f), 1);
+ FOR_EACH(j, 0, i) {
+ CHECK(f3_eq_eps(samps[j], samps[i], 1.e-6f), 0);
+ }
+ }
+
+ ssp_rng_ref_put(rng);
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+
+ CHECK(mem_allocated_size(), 0);
+ return 0;
+}
