commit 82da6ddf896e0797f14c662482d13c6379e44cc5
parent d45585825566b201f473ddfcd3671026bf46f1a4
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Wed, 15 Nov 2017 15:07:34 +0100
Merge branch 'release_0.5.0'
Diffstat:
24 files changed, 2361 insertions(+), 673 deletions(-)
diff --git a/README.md b/README.md
@@ -1,36 +1,51 @@
# Star SamPling
The purpose of this library is to generate [pseudo] random numbers and random
-variates. While it is based on C++11 random generators, its API remains pure C.
+variates. While it is partly based on C++11 random generators, its API remains
+pure C.
## How to build
The Star-Sampling library relies on the [CMake](http://www.cmake.org) and the
-[RCMake](https://gitlab.com/vaplv/rcmake/) package to build. It also depends on
-the [RSys](https://gitlab.com/vaplv/rsys/) library and C++11 random number
-generators.
-
-First ensure that CMake and a C++11 compiler is installed on your system. Then
-install the RCMake package as well as the RSys library. Finally Generate the
-project from the `cmake/CMakeLists.txt` file by appending to the
-`CMAKE_PREFIX_PATH` variable the `<RCMAKE_INSTALL_DIR>` and
-`<RSYS_INSTALL_DIR>` directories, where `<RCMAKE_INSTALL_DIR>`, and
-`<RSYS_INSTALL_DIR>` are the install directories of the RCMake package and the
-RSys library, respectively.
+[RCMake](https://gitlab.com/vaplv/rcmake/) packages to build. It also depends
+on the [RSys](https://gitlab.com/vaplv/rsys/) library and on random number
+generators from C++11 and the
+[Random123](https://github.com/DEShawResearch/Random123-Boost/) library.
+
+First ensure that CMake and a C++11 compiler are installed on your system. Then
+install the RCMake package as well as the RSys and Random123 libraries. Finally
+generate the project from the `cmake/CMakeLists.txt` file by appending to the
+`CMAKE_PREFIX_PATH` variable the `<RCMAKE_INSTALL_DIR>`, `<RSYS_INSTALL_DIR>`
+and `<RANDOM123_INSTALL_DIR>` directories, where `<RCMAKE_INSTALL_DIR>`,
+`<RSYS_INSTALL_DIR>` and `<RANDOM123_INSTALL_DIR>` are the install directories
+of the RCMake package and the RSys and Random123 libraries, respectively.
### Microsoft Visual Studio
-When built with Visual Studio, the Star-Sampling library depends on the
-[Boost](http://www.boost.org) random library rather than the STL that is
-actually not fully compliant with the C++11 standard.
+When building Star-Sampling using Visual Studio, the Star-Sampling library
+also depends on the [Boost](http://www.boost.org) random library as a
+replacement for the STL libray that is not fully compliant with the C++11
+standard on current Microsoft compilers.
-For this platform, install the Boost random library. Then generate the CMake
-project as above with, in addition, the Boost install directory added to the
-`CMAKE_PREFIX_PATH` variable and the `BOOST_INCLUDEDIR` cmake variable set to
-the path that contains the Boost header directory.
+For this platform install the Boost random library. Then generate the CMake
+project as above, excepted that you need to add the Boost install directory to
+the `CMAKE_PREFIX_PATH` variable and to set the `BOOST_INCLUDEDIR` cmake
+variable to the directory that contains the `boost` include directory.
## Release notes
+### Version 0.5
+
+- Rename the ssp_ran_uniform_disk API call into ssp_ran_uniform_disk_local.
+- Add a more general version of the uniform disk random variate allowing
+ users to provide the disk's normal.
+- Add a float equivalent for all the already defined double random variates.
+- Add some missing pdf API calls.
+- Change the API of some random variates that returned the pdf as an additional
+ vector component along with the sampled vector (i.e. filling up a vector[4]
+ instead of a vector[3]). The pdf is now returned through an optional
+ dedicated argument.
+
### Version 0.4
- Update the API of the random variates to return double precision reals
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -70,7 +70,7 @@ rcmake_append_runtime_dirs(_runtime_dirs RSys ${Boost_LIBRARY_DIRS})
# Configure and define targets
################################################################################
set(VERSION_MAJOR 0)
-set(VERSION_MINOR 4)
+set(VERSION_MINOR 5)
set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
@@ -126,7 +126,9 @@ if(NOT NO_TEST)
function(build_test _name)
add_executable(${_name}
${SSP_SOURCE_DIR}/${_name}.c
+ ${SSP_SOURCE_DIR}/${_name}.h
${SSP_SOURCE_DIR}/test_ssp_utils.h)
+ set_source_files_properties(${SSP_SOURCE_DIR}/${_name}.c PROPERTIES LANGUAGE CXX)
target_link_libraries(${_name} ssp RSys ${MATH_LIB})
set(_libraries ${ARGN})
foreach(_lib ${_libraries})
diff --git a/src/ssp.h b/src/ssp.h
@@ -33,6 +33,7 @@
#define SSP_H
#include <rsys/double33.h>
+#include <rsys/float33.h>
#include <rsys/math.h>
#include <rsys/rsys.h>
@@ -59,7 +60,9 @@ struct ssp_rng;
struct ssp_rng_proxy;
struct ssp_ranst_discrete;
struct ssp_ranst_gaussian;
+struct ssp_ranst_gaussian_float;
struct ssp_ranst_piecewise_linear;
+struct ssp_ranst_piecewise_linear_float;
/* Forward declaration of external types */
struct mem_allocator;
@@ -273,11 +276,21 @@ ssp_ran_exp
(struct ssp_rng* rng,
const double mu);
+SSP_API float
+ssp_ran_exp_float
+ (struct ssp_rng* rng,
+ const float mu);
+
SSP_API double
ssp_ran_exp_pdf
(const double x,
const double mu);
+SSP_API float
+ssp_ran_exp_float_pdf
+ (const float x,
+ const float mu);
+
/* Stateless random variate from the gaussian (or normal) distribution
* with mean `mu':
* P(x) dx = 1 / (sigma*sqrt(2*PI)) * exp(1/2*((x-mu)/sigma)^2) dx */
@@ -287,12 +300,25 @@ ssp_ran_gaussian
const double mu,
const double sigma);
+SSP_API float
+ssp_ran_gaussian_float
+ (struct ssp_rng* rng,
+ const float mu,
+ const float sigma);
+
+/* Return the probability distribution for a gaussian random variate */
SSP_API double
ssp_ran_gaussian_pdf
(const double x,
const double mu,
const double sigma);
+SSP_API float
+ssp_ran_gaussian_float_pdf
+ (const float x,
+ const float mu,
+ const float sigma);
+
/* Random variate from the lognormal distribution:
* P(x) dx = 1/(sigma*x*sqrt(2*PI)) * exp(-(ln(x)-zeta)^2 / (2*sigma^2)) dx */
SSP_API double
@@ -301,69 +327,78 @@ ssp_ran_lognormal
const double zeta,
const double sigma);
+SSP_API float
+ssp_ran_lognormal_float
+ (struct ssp_rng* rng,
+ const float zeta,
+ const float sigma);
+
+/* Return the probability distribution for a lognormal random variate */
SSP_API double
ssp_ran_lognormal_pdf
(const double x,
const double zeta,
const double sigma);
+SSP_API float
+ssp_ran_lognormal_float_pdf
+ (const float x,
+ const float zeta,
+ const float sigma);
+
/*******************************************************************************
* Sphere distribution
******************************************************************************/
/* Uniform sampling of an unit sphere. The PDF of the generated sample is
* stored in sample[3] */
-static INLINE double*
-ssp_ran_sphere_uniform(struct ssp_rng* rng, double 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] = cos(phi) * sin_theta;
- sample[1] = sin(phi) * sin_theta;
- sample[2] = cos_theta;
- sample[3] = 1.0/(4.0*PI);
- return sample;
-}
+SSP_API double*
+ssp_ran_sphere_uniform
+ (struct ssp_rng* rng,
+ double sample[3],
+ double* pdf); /* Can be NULL => pdf not computed */
+
+SSP_API float*
+ssp_ran_sphere_uniform_float
+ (struct ssp_rng* rng,
+ float sample[3],
+ float* pdf); /* Can be NULL => pdf not computed */
/* Return the probability distribution for a sphere uniform random variate */
static INLINE double
ssp_ran_sphere_uniform_pdf(void)
{
- return 1.0/(4.0*PI);
+ return 1 / (4*PI);
+}
+
+static INLINE float
+ssp_ran_sphere_uniform_float_pdf(void)
+{
+ return 1 / (4*(float)PI);
}
/*******************************************************************************
* Triangle distribution
******************************************************************************/
/* Uniform sampling of a triangle */
-static INLINE double*
+SSP_API double*
ssp_ran_triangle_uniform
(struct ssp_rng* rng,
const double v0[3], /* Position of the 1st vertex */
const double v1[3], /* Position of the 2nd vertex */
const double v2[3], /* Position of the 3rd vertex */
- double sample[3]) /* Sampled position */
-{
- double sqrt_u, v, one_minus_u;
- double vec0[3];
- double vec1[3];
- ASSERT(rng && v0 && v1 && v2 && sample);
+ double sample[3], /* Sampled position */
+ double* pdf); /* Can be NULL => pdf not computed */
- sqrt_u = sqrt(ssp_rng_canonical(rng));
- v = ssp_rng_canonical(rng);
- one_minus_u = 1.0 - sqrt_u;
-
- d3_sub(vec0, v0, v2);
- d3_sub(vec1, v1, v2);
- sample[0] = v2[0] + one_minus_u * vec0[0] + v * sqrt_u * vec1[0];
- sample[1] = v2[1] + one_minus_u * vec0[1] + v * sqrt_u * vec1[1];
- sample[2] = v2[2] + one_minus_u * vec0[2] + v * sqrt_u * vec1[2];
- return sample;
-}
+SSP_API float*
+ssp_ran_triangle_uniform_float
+ (struct ssp_rng* rng,
+ const float v0[3], /* Position of the 1st vertex */
+ const float v1[3], /* Position of the 2nd vertex */
+ const float v2[3], /* Position of the 3rd vertex */
+ float sample[3], /* Sampled position */
+ float* pdf); /* Can be NULL => pdf not computed */
+/* Return the probability distribution for a uniform point sampling on a triangle */
static INLINE double
ssp_ran_triangle_uniform_pdf
(const double v0[3],
@@ -375,28 +410,39 @@ ssp_ran_triangle_uniform_pdf
d3_sub(vec0, v0, v2);
d3_sub(vec1, v1, v2);
- return 1.0 / (d3_len(d3_cross(tmp, vec0, vec1)) * 0.5);
+ return 2 / d3_len(d3_cross(tmp, vec0, vec1));
+}
+
+static INLINE float
+ssp_ran_triangle_uniform_float_pdf
+ (const float v0[3],
+ const float v1[3],
+ const float v2[3])
+{
+ float vec0[3], vec1[3], tmp[3];
+ ASSERT(v0 && v1 && v2);
+
+ f3_sub(vec0, v0, v2);
+ f3_sub(vec1, v1, v2);
+ return 2 / f3_len(f3_cross(tmp, vec0, vec1));
}
/*******************************************************************************
* Hemisphere distribution
******************************************************************************/
/* Uniform sampling of an unit hemisphere whose up direction is implicitly the
- * Z axis. The PDF of the generated sample is stored in sample[3] */
-static INLINE double*
-ssp_ran_hemisphere_uniform_local(struct ssp_rng* rng, double sample[4])
-{
- double phi, cos_theta, sin_theta;
- ASSERT(rng && sample);
- phi = 2.0 * PI * ssp_rng_canonical(rng);
- cos_theta = ssp_rng_canonical(rng);
- sin_theta = cos2sin(cos_theta);
- sample[0] = cos(phi) * sin_theta;
- sample[1] = sin(phi) * sin_theta;
- sample[2] = cos_theta;
- sample[3] = 1.0/(2.0*PI);
- return sample;
-}
+ * Z axis. */
+SSP_API double*
+ssp_ran_hemisphere_uniform_local
+ (struct ssp_rng* rng,
+ double sample[3],
+ double* pdf); /* Can be NULL => pdf not computed */
+
+SSP_API float*
+ssp_ran_hemisphere_uniform_float_local
+ (struct ssp_rng* rng,
+ float sample[3],
+ float* pdf); /* Can be NULL => pdf not computed */
/* Return the probability distribution for an hemispheric uniform random
* variate with an implicit up direction in Z */
@@ -404,50 +450,79 @@ static INLINE double
ssp_ran_hemisphere_uniform_local_pdf(const double sample[3])
{
ASSERT(sample);
- return sample[2] < 0.f ? 0.f : 1.0/(2.0*PI);
+ return sample[2] < 0 ? 0 : 1/(2*PI);
+}
+
+static INLINE float
+ssp_ran_hemisphere_uniform_float_local_pdf(const float sample[3])
+{
+ ASSERT(sample);
+ return sample[2] < 0 ? 0 : 1/(2*(float)PI);
}
-/* Uniform sampling of an unit hemisphere whose up direction is `up'. The PDF
- * of the generated sample is stored in sample[3] */
+/* Uniform sampling of an unit hemisphere whose up direction is `up'. */
static INLINE double*
ssp_ran_hemisphere_uniform
- (struct ssp_rng* rng, const double up[3], double sample[4])
+ (struct ssp_rng* rng,
+ const double up[3],
+ double sample[3],
+ double* pdf) /* Can be NULL => pdf not computed */
{
double basis[9];
- double sample_local[4];
+ double sample_local[3];
ASSERT(rng && up && sample && d3_is_normalized(up));
- ssp_ran_hemisphere_uniform_local(rng, sample_local);
- sample[3] = sample_local[3];
+ ssp_ran_hemisphere_uniform_local(rng, sample_local, pdf);
return d33_muld3(sample, d33_basis(basis, up), sample_local);
}
+static INLINE float*
+ssp_ran_hemisphere_uniform_float
+ (struct ssp_rng* rng,
+ const float up[3],
+ float sample[3],
+ float* pdf) /* Can be NULL => pdf not computed */
+{
+ float basis[9];
+ float sample_local[3];
+ ASSERT(rng && up && sample && f3_is_normalized(up));
+ ssp_ran_hemisphere_uniform_float_local(rng, sample_local, pdf);
+ return f33_mulf3(sample, f33_basis(basis, up), sample_local);
+}
+
/* Return the probability distribution for an hemispheric uniform random
* variate with an explicit `up' direction */
static INLINE double
-ssp_ran_hemisphere_uniform_pdf(const double up[3], const double sample[3])
+ssp_ran_hemisphere_uniform_pdf
+ (const double up[3],
+ const double sample[3])
{
ASSERT(up && sample && d3_is_normalized(up));
- return d3_dot(sample, up) < 0.f ? 0.f : 1.0/(2.0*PI);
+ return d3_dot(sample, up) < 0 ? 0 : 1/(2*PI);
+}
+
+static INLINE float
+ssp_ran_hemisphere_uniform_float_pdf
+ (const float up[3],
+ const float sample[3])
+{
+ ASSERT(up && sample && f3_is_normalized(up));
+ return f3_dot(sample, up) < 0 ? 0 : 1/(2*(float)PI);
}
/* Cosine weighted sampling of an unit hemisphere whose up direction is
* implicitly the Z axis. The PDF of the generated sample is stored in
* sample[3] */
-static INLINE double*
-ssp_ran_hemisphere_cos_local(struct ssp_rng* rng, double sample[4])
-{
- double phi, cos_theta, sin_theta, v;
- ASSERT(rng && sample);
- phi = 2.0 * PI * ssp_rng_canonical(rng);
- v = ssp_rng_canonical(rng);
- cos_theta = sqrt(v);
- sin_theta = sqrt(1.0 - v);
- sample[0] = cos(phi) * sin_theta;
- sample[1] = sin(phi) * sin_theta;
- sample[2] = cos_theta;
- sample[3] = cos_theta / PI;
- return sample;
-}
+SSP_API double*
+ssp_ran_hemisphere_cos_local
+ (struct ssp_rng* rng,
+ double sample[3],
+ double* pdf); /* Can be NULL => pdf not computed */
+
+SSP_API float*
+ssp_ran_hemisphere_cos_float_local
+ (struct ssp_rng* rng,
+ float sample[3],
+ float* pdf); /* Can be NULL => pdf not computed */
/* Return the probability distribution for an hemispheric cosine weighted
* random variate with an implicit up direction in Z */
@@ -455,111 +530,145 @@ static INLINE double
ssp_ran_hemisphere_cos_local_pdf(const double sample[3])
{
ASSERT(sample);
- return sample[2] < 0.f ? 0.f : sample[2] / PI;
+ return sample[2] < 0 ? 0 : sample[2]/PI;
}
-/* Cosine weighted sampling of an unit hemisphere whose up direction is `up'.
- * The PDF of the generated sample is stored in sample[3] */
+static INLINE float
+ssp_ran_hemisphere_cos_float_local_pdf(const float sample[3])
+{
+ ASSERT(sample);
+ return sample[2] < 0 ? 0 : sample[2]/(float)PI;
+}
+
+/* Cosine weighted sampling of an unit hemisphere whose up direction is `up'. */
static INLINE double*
-ssp_ran_hemisphere_cos(struct ssp_rng* rng, const double up[3], double sample[4])
+ssp_ran_hemisphere_cos
+ (struct ssp_rng* rng,
+ const double up[3],
+ double sample[3],
+ double* pdf) /* Can be NULL => pdf not computed */
{
- double sample_local[4];
+ double sample_local[3];
double basis[9];
ASSERT(rng && up && sample && d3_is_normalized(up));
- ssp_ran_hemisphere_cos_local(rng, sample_local);
- sample[3] = sample_local[3];
+ ssp_ran_hemisphere_cos_local(rng, sample_local, pdf);
return d33_muld3(sample, d33_basis(basis, up), sample_local);
}
+static INLINE float*
+ssp_ran_hemisphere_cos_float
+ (struct ssp_rng* rng,
+ const float up[3],
+ float sample[3],
+ float* pdf) /* Can be NULL => pdf not computed */
+{
+ float sample_local[3];
+ float basis[9];
+ ASSERT(rng && up && sample && f3_is_normalized(up));
+ ssp_ran_hemisphere_cos_float_local(rng, sample_local, pdf);
+ return f33_mulf3(sample, f33_basis(basis, up), sample_local);
+}
+
/* Return the probability distribution for an hemispheric cosine weighted
* random variate with an explicit `up' direction */
static INLINE double
-ssp_ran_hemisphere_cos_pdf(const double up[3], const double sample[3])
+ssp_ran_hemisphere_cos_pdf
+ (const double up[3],
+ const double sample[3])
{
double dot;
ASSERT(up && sample);
dot = d3_dot(sample, up);
- return dot < 0.f ? 0.f : dot/PI;
+ return dot < 0 ? 0 : dot/PI;
+}
+
+static INLINE float
+ssp_ran_hemisphere_cos_float_pdf
+ (const float up[3],
+ const float sample[3])
+{
+ float dot;
+ ASSERT(up && sample);
+ dot = f3_dot(sample, up);
+ return dot < 0 ? 0 : dot/(float)PI;
}
/*******************************************************************************
* Henyey Greenstein distribution
******************************************************************************/
+/* Henyey-Greenstein sampling of an unit sphere for an incident direction
+ * that is implicitly the Z axis. */
+SSP_API double*
+ssp_ran_sphere_hg_local
+ (struct ssp_rng* rng,
+ const double g,
+ double sample[3],
+ double* pdf); /* Can be NULL => pdf not computed */
+
+SSP_API float*
+ssp_ran_sphere_hg_float_local
+ (struct ssp_rng* rng,
+ const float g,
+ float sample[3],
+ float* pdf); /* Can be NULL => pdf not computed */
+
/* Return the probability distribution for a Henyey-Greenstein random
- * variate with an implicit incident direction in Z */
-static INLINE double
-ssp_ran_sphere_hg_local_pdf(const double g, const double sample[3])
-{
- double epsilon_g, epsilon_mu;
- ASSERT(sample && d3_is_normalized(sample));
- if (g>0) {
- epsilon_g = 1-g;
- epsilon_mu = 1-sample[2];
- } else {
- epsilon_g = 1+g;
- epsilon_mu = 1+sample[2];
- }
- return 1/(4*PI)
- *epsilon_g*(2-epsilon_g)
- *pow(epsilon_g*epsilon_g+2*epsilon_mu*(1-epsilon_g),-1.5);
-}
+* variate with an implicit incident direction in Z */
+SSP_API double
+ssp_ran_sphere_hg_local_pdf
+ (const double g,
+ const double sample[3]);
+
+SSP_API float
+ssp_ran_sphere_hg_float_local_pdf
+ (const float g,
+ const float sample[3]);
/* Henyey-Greenstein sampling of an unit sphere for an incident direction
- * that is implicitly the Z axis.
- * The PDF of the generated sample is stored in sample[3] */
+ * that is `up'. */
static INLINE double*
-ssp_ran_sphere_hg_local(struct ssp_rng* rng, const double g, double sample[4])
+ssp_ran_sphere_hg
+ (struct ssp_rng* rng,
+ const double up[3],
+ const double g,
+ double sample[3],
+ double* pdf) /* Can be NULL => pdf not computed */
{
- double phi, R, cos_theta, sin_theta;
- ASSERT(rng && sample);
- phi = 2 * PI * ssp_rng_canonical(rng);
- R = ssp_rng_canonical(rng);
- cos_theta = 2*R*(1+g)*(1+g)*(g*(R-1)+1)
- /((1-g*(1-2*R))*(1-g*(1-2*R)))-1;
- sin_theta = cos2sin(cos_theta);
- sample[0] = cos(phi) * sin_theta;
- sample[1] = sin(phi) * sin_theta;
- sample[2] = cos_theta;
- sample[3] = ssp_ran_sphere_hg_local_pdf(g,sample);
- return sample;
+ double sample_local[3];
+ double basis[9];
+ ASSERT(-1 <= g && g <= +1 && rng && up && sample && d3_is_normalized(up));
+ ssp_ran_sphere_hg_local(rng, g, sample_local, pdf);
+ return d33_muld3(sample, d33_basis(basis, up), sample_local);
+}
+
+static INLINE float*
+ssp_ran_sphere_hg_float
+ (struct ssp_rng* rng,
+ const float up[3],
+ const float g,
+ float sample[3],
+ float* pdf) /* Can be NULL => pdf not computed */
+{
+ float sample_local[3];
+ float basis[9];
+ ASSERT(-1 <= g && g <= +1 && rng && up && sample && f3_is_normalized(up));
+ ssp_ran_sphere_hg_float_local(rng, g, sample_local, pdf);
+ return f33_mulf3(sample, f33_basis(basis, up), sample_local);
}
/* Return the probability distribution for a Henyey-Greenstein random
- * variate with an explicit incident direction that is `up' */
-static INLINE double
+* variate with an explicit incident direction that is `up' */
+SSP_API double
ssp_ran_sphere_hg_pdf
(const double up[3],
const double g,
- const double sample[3])
-{
- double epsilon_g, epsilon_mu;
- ASSERT(up && sample && d3_is_normalized(up) && d3_is_normalized(sample));
- if (g>0) {
- epsilon_g = 1-g;
- epsilon_mu = 1-d3_dot(sample,up);
- } else {
- epsilon_g = 1+g;
- epsilon_mu = 1+d3_dot(sample,up);
- }
- return 1/(4*PI)
- *epsilon_g*(2-epsilon_g)
- *pow(epsilon_g*epsilon_g+2*epsilon_mu*(1-epsilon_g),-1.5);
-}
+ const double sample[3]);
-/* Henyey-Greenstein sampling of an unit sphere for an incident direction
- * that is `up'.
- * The PDF of the generated sample is stored in sample[3] */
-static INLINE double*
-ssp_ran_sphere_hg
- (struct ssp_rng* rng, const double up[3], const double g, double sample[4])
-{
- double sample_local[4];
- double basis[9];
- ASSERT(rng && up && sample && d3_is_normalized(up));
- ssp_ran_sphere_hg_local(rng, g, sample_local);
- sample[3] = sample_local[3];
- return d33_muld3(sample, d33_basis(basis, up), sample_local);
-}
+SSP_API float
+ssp_ran_sphere_hg_float_pdf
+ (const float up[3],
+ const float g,
+ const float sample[3]);
/*******************************************************************************
* General discrete distribution with state
@@ -607,12 +716,27 @@ ssp_ranst_gaussian_create
struct ssp_ranst_gaussian** ran);
SSP_API res_T
+ssp_ranst_gaussian_float_create
+ (struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */
+ struct ssp_ranst_gaussian_float** ran);
+
+SSP_API res_T
ssp_ranst_gaussian_setup
(struct ssp_ranst_gaussian* ran,
const double mu,
const double sigma);
SSP_API res_T
+ssp_ranst_gaussian_float_setup
+ (struct ssp_ranst_gaussian_float* ran,
+ const float mu,
+ const float sigma);
+
+SSP_API res_T
+ssp_ranst_gaussian_float_ref_get
+ (struct ssp_ranst_gaussian_float* ran);
+
+SSP_API res_T
ssp_ranst_gaussian_ref_get
(struct ssp_ranst_gaussian* ran);
@@ -620,16 +744,30 @@ SSP_API res_T
ssp_ranst_gaussian_ref_put
(struct ssp_ranst_gaussian* ran);
+SSP_API res_T
+ssp_ranst_gaussian_float_ref_put
+ (struct ssp_ranst_gaussian_float* ran);
+
SSP_API double
ssp_ranst_gaussian_get
(const struct ssp_ranst_gaussian* ran,
struct ssp_rng* rng);
+SSP_API float
+ssp_ranst_gaussian_float_get
+ (const struct ssp_ranst_gaussian_float* ran,
+ struct ssp_rng* rng);
+
SSP_API double
ssp_ranst_gaussian_pdf
(const struct ssp_ranst_gaussian* ran,
const double x);
+SSP_API float
+ssp_ranst_gaussian_float_pdf
+ (const struct ssp_ranst_gaussian_float* ran,
+ const float x);
+
/*******************************************************************************
* Piecewise linear distribution with state
******************************************************************************/
@@ -639,6 +777,11 @@ ssp_ranst_piecewise_linear_create
struct ssp_ranst_piecewise_linear **ran);
SSP_API res_T
+ssp_ranst_piecewise_linear_float_create
+ (struct mem_allocator* allocator, /* May be NULL <=> Use default allocator */
+ struct ssp_ranst_piecewise_linear_float **ran);
+
+SSP_API res_T
ssp_ranst_piecewise_linear_setup
(struct ssp_ranst_piecewise_linear *ran,
const double* intervals,
@@ -646,41 +789,107 @@ ssp_ranst_piecewise_linear_setup
const size_t size);
SSP_API res_T
+ssp_ranst_piecewise_linear_float_setup
+ (struct ssp_ranst_piecewise_linear_float *ran,
+ const float* intervals,
+ const float* weights,
+ const size_t size);
+
+SSP_API res_T
ssp_ranst_piecewise_linear_ref_get
(struct ssp_ranst_piecewise_linear* ran);
SSP_API res_T
+ssp_ranst_piecewise_linear_float_ref_get
+ (struct ssp_ranst_piecewise_linear_float* ran);
+
+SSP_API res_T
ssp_ranst_piecewise_linear_ref_put
(struct ssp_ranst_piecewise_linear* ran);
+SSP_API res_T
+ssp_ranst_piecewise_linear_float_ref_put
+ (struct ssp_ranst_piecewise_linear_float* ran);
+
SSP_API double
ssp_ranst_piecewise_linear_get
(const struct ssp_ranst_piecewise_linear *ran,
struct ssp_rng* rng);
+SSP_API float
+ssp_ranst_piecewise_linear_float_get
+ (const struct ssp_ranst_piecewise_linear_float *ran,
+ struct ssp_rng* rng);
+
+SSP_API double
+ssp_ranst_piecewise_linear_pdf
+ (const struct ssp_ranst_piecewise_linear *ran,
+ double x);
+
+SSP_API float
+ssp_ranst_piecewise_linear_float_pdf
+ (const struct ssp_ranst_piecewise_linear_float *ran,
+ float x);
+
/*******************************************************************************
* 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_API double*
+ssp_ran_uniform_disk_local
+ (struct ssp_rng* rng,
+ const double radius,
+ double pt[3], /* pt[2] <= 0 */
+ double* pdf); /* Can be NULL => pdf not computed */
+
+SSP_API float*
+ssp_ran_uniform_disk_float_local
+ (struct ssp_rng* rng,
+ const float radius,
+ float pt[3], /* pt[2] <= 0 */
+ float* pdf); /* Can be NULL => pdf not computed */
+
+static INLINE double
+ssp_ran_uniform_disk_local_pdf(const double radius)
+{
+ return 1 / (radius * radius);
+}
+
+static INLINE float
+ssp_ran_uniform_disk_float_local_pdf(const float radius)
+{
+ return 1 / (radius * radius);
+}
+
+static INLINE double*
ssp_ran_uniform_disk
(struct ssp_rng* rng,
const double radius,
- double pt[2])
+ const double up[3],
+ double pt[3],
+ double* pdf) /* Can be NULL => pdf not computed */
+{
+ double sample_local[3];
+ double basis[9];
+ ASSERT(rng && up && pt && radius > 0);
+ ssp_ran_uniform_disk_local(rng, radius, sample_local, pdf);
+ return d33_muld3(pt, d33_basis(basis, up), sample_local);
+}
+
+static INLINE float*
+ssp_ran_uniform_disk_float
+ (struct ssp_rng* rng,
+ const float radius,
+ const float up[3],
+ float pt[3],
+ float* pdf) /* Can be NULL => pdf not computed */
{
- double theta, r;
- ASSERT(rng && pt && radius > 0);
- theta = ssp_rng_uniform_double(rng, 0, 2 * PI);
- r = radius * sqrt(ssp_rng_canonical(rng));
- pt[0] = r * cos(theta);
- pt[1] = r * sin(theta);
- return pt;
+ float sample_local[3];
+ float basis[9];
+ ASSERT(rng && up && pt && radius > 0);
+ ssp_ran_uniform_disk_float_local(rng, radius, sample_local, pdf);
+ return f33_mulf3(pt, f33_basis(basis, up), sample_local);
}
END_DECLS
#endif /* SSP_H */
-
diff --git a/src/ssp_ran.c b/src/ssp_ran.c
@@ -31,25 +31,53 @@
#include "ssp_rng_c.h"
+static FINLINE float
+sinf2cosf(const float d)
+{
+ return sqrtf(MMAX(0, 1 - d*d));
+}
+
+static FINLINE float
+cosf2sinf(const float d)
+{
+ return sinf2cosf(d);
+}
+
/*******************************************************************************
* Exported state free random variates
******************************************************************************/
double
ssp_ran_exp(struct ssp_rng* rng, const double mu)
{
- ASSERT(rng);
+ ASSERT(rng && mu > 0);
rng_cxx rng_cxx(*rng);
RAN_NAMESPACE::exponential_distribution<double> distribution(mu);
return distribution(rng_cxx);
}
+float
+ssp_ran_exp_float(struct ssp_rng* rng, const float mu)
+{
+ ASSERT(rng && mu > 0);
+ rng_cxx rng_cxx(*rng);
+ RAN_NAMESPACE::exponential_distribution<float> distribution(mu);
+ return distribution(rng_cxx);
+}
+
double
ssp_ran_exp_pdf(const double x, const double mu)
{
- ASSERT(x >= 0);
+ ASSERT(x >= 0 && mu > 0);
return mu * exp(-x * mu);
}
+float
+ssp_ran_exp_float_pdf(const float x, const float mu)
+{
+ ASSERT(x >= 0 && mu > 0);
+ return mu * expf(-x * mu);
+}
+
double
ssp_ran_gaussian
(struct ssp_rng* rng,
@@ -62,11 +90,36 @@ ssp_ran_gaussian
return distribution(rng_cxx);
}
+float
+ssp_ran_gaussian_float
+ (struct ssp_rng* rng,
+ const float mu,
+ const float sigma)
+{
+ ASSERT(rng);
+ rng_cxx rng_cxx(*rng);
+ RAN_NAMESPACE::normal_distribution<float> distribution(mu, sigma);
+ return distribution(rng_cxx);
+}
+
double
-ssp_ran_gaussian_pdf(const double x, const double mu, const double sigma)
+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);
+ return 1 / (sigma * SQRT_2_PI) * exp(-0.5 * tmp * tmp);
+}
+
+float
+ssp_ran_gaussian_float_pdf
+ (const float x,
+ const float mu,
+ const float sigma)
+{
+ const float tmp = (x - mu) / sigma;
+ return 1 / (sigma * (float)SQRT_2_PI) * expf(-0.5f * tmp * tmp);
}
double
@@ -81,10 +134,360 @@ ssp_ran_lognormal
return distribution(rng_cxx);
}
+float
+ssp_ran_lognormal_float
+ (struct ssp_rng* rng,
+ const float zeta,
+ const float sigma)
+{
+ ASSERT(rng);
+ rng_cxx rng_cxx(*rng);
+ RAN_NAMESPACE::lognormal_distribution<float> distribution(zeta, sigma);
+ return distribution(rng_cxx);
+}
+
double
-ssp_ran_lognormal_pdf(const double x, const double zeta, const double sigma)
+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));
+ return 1 / (sigma * x * SQRT_2_PI) * exp(-(tmp*tmp) / (2*sigma*sigma));
+}
+
+float
+ssp_ran_lognormal_float_pdf
+ (const float x,
+ const float zeta,
+ const float sigma)
+{
+ const float tmp = logf(x) - zeta;
+ return 1 / (sigma * x * (float)SQRT_2_PI) * expf(-(tmp*tmp) / (2 * sigma*sigma));
}
+double*
+ssp_ran_sphere_uniform
+ (struct ssp_rng* rng,
+ double sample[3],
+ double* pdf)
+{
+ double phi, cos_theta, sin_theta, v;
+ ASSERT(rng && sample);
+ phi = ssp_rng_uniform_double(rng, 0, 2 * PI);
+ v = ssp_rng_canonical(rng);
+ cos_theta = 1 - 2 * v;
+ sin_theta = 2 * sqrt(v * (1 - v));
+ sample[0] = cos(phi) * sin_theta;
+ sample[1] = sin(phi) * sin_theta;
+ sample[2] = cos_theta;
+ if(pdf) *pdf = 1 / (4*PI);
+ return sample;
+}
+
+float*
+ssp_ran_sphere_uniform_float
+ (struct ssp_rng* rng,
+ float sample[3],
+ float* pdf)
+{
+ float phi, cos_theta, sin_theta, v;
+ ASSERT(rng && sample);
+ phi = ssp_rng_uniform_float(rng, 0, 2 * (float)PI);
+ v = ssp_rng_canonical_float(rng);
+ cos_theta = 1 - 2 * v;
+ sin_theta = 2 * sqrtf(v * (1 - v));
+ sample[0] = cosf(phi) * sin_theta;
+ sample[1] = sinf(phi) * sin_theta;
+ sample[2] = cos_theta;
+ if (pdf) *pdf = 1 / (4*(float)PI);
+ return sample;
+}
+
+double*
+ssp_ran_triangle_uniform
+ (struct ssp_rng* rng,
+ const double v0[3],
+ const double v1[3],
+ const double v2[3],
+ double sample[3],
+ double* pdf)
+{
+ double sqrt_u, v, one_minus_u;
+ double vec0[3];
+ double vec1[3];
+ double tmp[3];
+ ASSERT(rng && v0 && v1 && v2 && sample);
+
+ sqrt_u = sqrt(ssp_rng_canonical(rng));
+ v = ssp_rng_canonical(rng);
+ one_minus_u = 1.0 - sqrt_u;
+
+ d3_sub(vec0, v0, v2);
+ d3_sub(vec1, v1, v2);
+ sample[0] = v2[0] + one_minus_u * vec0[0] + v * sqrt_u * vec1[0];
+ sample[1] = v2[1] + one_minus_u * vec0[1] + v * sqrt_u * vec1[1];
+ sample[2] = v2[2] + one_minus_u * vec0[2] + v * sqrt_u * vec1[2];
+ if (pdf) *pdf = 2 / d3_len(d3_cross(tmp, vec0, vec1));
+ return sample;
+}
+
+float*
+ssp_ran_triangle_uniform_float
+ (struct ssp_rng* rng,
+ const float v0[3],
+ const float v1[3],
+ const float v2[3],
+ float sample[3],
+ float* pdf)
+{
+ float sqrt_u, v, one_minus_u;
+ float vec0[3];
+ float vec1[3];
+ float tmp[3];
+ ASSERT(rng && v0 && v1 && v2 && sample);
+
+ sqrt_u = sqrtf(ssp_rng_canonical_float(rng));
+ v = ssp_rng_canonical_float(rng);
+ one_minus_u = 1 - sqrt_u;
+
+ f3_sub(vec0, v0, v2);
+ f3_sub(vec1, v1, v2);
+ sample[0] = v2[0] + one_minus_u * vec0[0] + v * sqrt_u * vec1[0];
+ sample[1] = v2[1] + one_minus_u * vec0[1] + v * sqrt_u * vec1[1];
+ sample[2] = v2[2] + one_minus_u * vec0[2] + v * sqrt_u * vec1[2];
+ if (pdf) *pdf = 2 / f3_len(f3_cross(tmp, vec0, vec1));
+ return sample;
+}
+
+double*
+ssp_ran_hemisphere_uniform_local
+ (struct ssp_rng* rng,
+ double sample[3],
+ double* pdf)
+{
+ double phi, cos_theta, sin_theta;
+ ASSERT(rng && sample);
+ phi = ssp_rng_uniform_double(rng, 0, 2 * PI);
+ cos_theta = ssp_rng_canonical(rng);
+ sin_theta = cos2sin(cos_theta);
+ sample[0] = cos(phi) * sin_theta;
+ sample[1] = sin(phi) * sin_theta;
+ sample[2] = cos_theta;
+ if (pdf) *pdf = 1 / (2*PI);
+ return sample;
+}
+
+float*
+ssp_ran_hemisphere_uniform_float_local
+ (struct ssp_rng* rng,
+ float sample[3],
+ float* pdf)
+{
+ float phi, cos_theta, sin_theta;
+ ASSERT(rng && sample);
+ phi = ssp_rng_uniform_float(rng, 0, 2 * (float)PI);
+ cos_theta = ssp_rng_canonical_float(rng);
+ sin_theta = cosf2sinf(cos_theta);
+ sample[0] = cosf(phi) * sin_theta;
+ sample[1] = sinf(phi) * sin_theta;
+ sample[2] = cos_theta;
+ if (pdf) *pdf = 1 / (2*(float)PI);
+ return sample;
+}
+
+double*
+ssp_ran_hemisphere_cos_local
+ (struct ssp_rng* rng,
+ double sample[3],
+ double* pdf)
+{
+ double phi, cos_theta, sin_theta, v;
+ ASSERT(rng && sample);
+ phi = ssp_rng_uniform_double(rng, 0, 2 * PI);
+ v = ssp_rng_canonical(rng);
+ cos_theta = sqrt(v);
+ sin_theta = sqrt(1 - v);
+ sample[0] = cos(phi) * sin_theta;
+ sample[1] = sin(phi) * sin_theta;
+ sample[2] = cos_theta;
+ if (pdf) *pdf = cos_theta / PI;
+ return sample;
+}
+
+float*
+ssp_ran_hemisphere_cos_float_local
+ (struct ssp_rng* rng,
+ float sample[3],
+ float* pdf)
+{
+ float phi, cos_theta, sin_theta, v;
+ ASSERT(rng && sample);
+ phi = ssp_rng_uniform_float(rng, 0, 2 * (float)PI);
+ v = ssp_rng_canonical_float(rng);
+ cos_theta = sqrtf(v);
+ sin_theta = sqrtf(1 - v);
+ sample[0] = cosf(phi) * sin_theta;
+ sample[1] = sinf(phi) * sin_theta;
+ sample[2] = cos_theta;
+ if (pdf) *pdf = cos_theta / (float)PI;
+ return sample;
+}
+
+double
+ssp_ran_sphere_hg_local_pdf
+ (const double g,
+ const double sample[3])
+{
+ double epsilon_g, epsilon_mu;
+ ASSERT(fabs(g) <= 1 && sample && d3_is_normalized(sample));
+ if(g>0) {
+ epsilon_g = 1 - g;
+ epsilon_mu = 1 - sample[2];
+ } else {
+ epsilon_g = 1 + g;
+ epsilon_mu = 1 + sample[2];
+ }
+ return 1 / (4 * PI)
+ *epsilon_g*(2 - epsilon_g)
+ *pow(epsilon_g*epsilon_g + 2 * epsilon_mu*(1 - epsilon_g), -1.5);
+}
+
+float
+ssp_ran_sphere_hg_float_local_pdf
+ (const float g,
+ const float sample[3])
+{
+ float epsilon_g, epsilon_mu;
+ ASSERT(fabsf(g) <= 1 && sample && f3_is_normalized(sample));
+ if(g>0) {
+ epsilon_g = 1 - g;
+ epsilon_mu = 1 - sample[2];
+ } else {
+ epsilon_g = 1 + g;
+ epsilon_mu = 1 + sample[2];
+ }
+ return 1 / (4 * (float)PI)
+ *epsilon_g*(2 - epsilon_g)
+ *powf(epsilon_g*epsilon_g + 2 * epsilon_mu*(1 - epsilon_g), -1.5f);
+}
+
+double*
+ssp_ran_sphere_hg_local
+ (struct ssp_rng* rng,
+ const double g,
+ double sample[3],
+ double* pdf)
+{
+ double phi, R, cos_theta, sin_theta;
+ ASSERT(fabs(g) <= 1 && rng && sample);
+ phi = ssp_rng_uniform_double(rng, 0, 2 * PI);
+ R = ssp_rng_canonical(rng);
+ cos_theta = 2 * R*(1 + g)*(1 + g)*(g*(R - 1) + 1)
+ / ((1 - g*(1 - 2 * R))*(1 - g*(1 - 2 * R))) - 1;
+ sin_theta = cos2sin(cos_theta);
+ sample[0] = cos(phi) * sin_theta;
+ sample[1] = sin(phi) * sin_theta;
+ sample[2] = cos_theta;
+ if (pdf) *pdf = ssp_ran_sphere_hg_local_pdf(g, sample);
+ return sample;
+}
+
+float*
+ssp_ran_sphere_hg_float_local
+ (struct ssp_rng* rng,
+ const float g,
+ float sample[3],
+ float* pdf)
+{
+ float phi, R, cos_theta, sin_theta;
+ ASSERT(fabsf(g) <= 1 && rng && sample);
+ phi = ssp_rng_uniform_float(rng, 0, 2 * (float)PI);
+ R = ssp_rng_canonical_float(rng);
+ cos_theta = 2 * R*(1 + g)*(1 + g)*(g*(R - 1) + 1)
+ / ((1 - g*(1 - 2 * R))*(1 - g*(1 - 2 * R))) - 1;
+ sin_theta = cosf2sinf(cos_theta);
+ sample[0] = cosf(phi) * sin_theta;
+ sample[1] = sinf(phi) * sin_theta;
+ sample[2] = cos_theta;
+ if (pdf) *pdf = ssp_ran_sphere_hg_float_local_pdf(g, sample);
+ return sample;
+}
+
+double
+ssp_ran_sphere_hg_pdf
+ (const double up[3],
+ const double g,
+ const double sample[3])
+{
+ double epsilon_g, epsilon_mu;
+ ASSERT(-1 <= g && g <= +1 &&
+ up && sample && d3_is_normalized(up) && d3_is_normalized(sample));
+ if(g>0) {
+ epsilon_g = 1 - g;
+ epsilon_mu = 1 - d3_dot(sample, up);
+ } else {
+ epsilon_g = 1 + g;
+ epsilon_mu = 1 + d3_dot(sample, up);
+ }
+ return 1 / (4 * PI)
+ *epsilon_g*(2 - epsilon_g)
+ *pow(epsilon_g*epsilon_g + 2 * epsilon_mu*(1 - epsilon_g), -1.5);
+}
+
+float
+ssp_ran_sphere_hg_float_pdf
+ (const float up[3],
+ const float g,
+ const float sample[3])
+{
+ float epsilon_g, epsilon_mu;
+ ASSERT(-1 <= g && g <= +1 &&
+ up && sample && f3_is_normalized(up) && f3_is_normalized(sample));
+ if(g>0) {
+ epsilon_g = 1 - g;
+ epsilon_mu = 1 - f3_dot(sample, up);
+ } else {
+ epsilon_g = 1 + g;
+ epsilon_mu = 1 + f3_dot(sample, up);
+ }
+ return 1 / (4 * (float)PI)
+ *epsilon_g*(2 - epsilon_g)
+ *powf(epsilon_g*epsilon_g + 2 * epsilon_mu*(1 - epsilon_g), -1.5f);
+}
+
+double*
+ssp_ran_uniform_disk_local
+ (struct ssp_rng* rng,
+ const double radius,
+ double pt[3],
+ double* pdf)
+{
+ double theta, r;
+ ASSERT(rng && pt && radius > 0);
+ theta = ssp_rng_uniform_double(rng, 0, 2 * PI);
+ r = radius * sqrt(ssp_rng_canonical(rng));
+ pt[0] = r * cos(theta);
+ pt[1] = r * sin(theta);
+ pt[2] = 0;
+ if (pdf) *pdf = 1 / (radius * radius);
+ return pt;
+}
+
+float*
+ssp_ran_uniform_disk_float_local
+ (struct ssp_rng* rng,
+ const float radius,
+ float pt[3],
+ float* pdf)
+{
+ float theta, r;
+ ASSERT(rng && pt && radius > 0);
+ theta = ssp_rng_uniform_float(rng, 0, 2 * (float)PI);
+ r = radius * sqrtf(ssp_rng_canonical_float(rng));
+ pt[0] = r * cosf(theta);
+ pt[1] = r * sinf(theta);
+ pt[2] = 0;
+ if (pdf) *pdf = 1 / (radius * radius);
+ return pt;
+}
+\ No newline at end of file
diff --git a/src/ssp_ranst_gaussian.c b/src/ssp_ranst_gaussian.c
@@ -36,6 +36,8 @@
using normal_dist = RAN_NAMESPACE::normal_distribution<double>;
+using normal_dist_float = RAN_NAMESPACE::normal_distribution<float>;
+
struct ssp_ranst_gaussian {
ref_T ref;
struct mem_allocator* allocator;
@@ -45,6 +47,15 @@ struct ssp_ranst_gaussian {
double K2; /* 1.0 / (sigma * SQRT_2_PI) */
};
+struct ssp_ranst_gaussian_float {
+ ref_T ref;
+ struct mem_allocator* allocator;
+ normal_dist_float* distrib;
+ float mu;
+ float K1; /* 1.0 / sigma */
+ float K2; /* 1.0 / (sigma * SQRT_2_PI) */
+};
+
/*******************************************************************************
* Helper function
******************************************************************************/
@@ -61,6 +72,19 @@ gaussian_release(ref_T* ref)
MEM_RM(ran->allocator, ran);
}
+static void
+gaussian_float_release(ref_T* ref)
+{
+ ssp_ranst_gaussian_float* ran;
+ ASSERT(ref);
+ ran = CONTAINER_OF(ref, ssp_ranst_gaussian_float, ref);
+ if(ran->distrib) {
+ ran->distrib->~normal_dist_float();
+ MEM_RM(ran->allocator, ran->distrib);
+ }
+ MEM_RM(ran->allocator, ran);
+}
+
/*******************************************************************************
* Exported functions
******************************************************************************/
@@ -107,6 +131,48 @@ error:
}
res_T
+ssp_ranst_gaussian_float_create
+ (struct mem_allocator* allocator,
+ struct ssp_ranst_gaussian_float** out_ran)
+{
+ ssp_ranst_gaussian_float* 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_float*>
+ (MEM_CALLOC(allocator, 1, sizeof(struct ssp_ranst_gaussian_float)));
+ if(!ran) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ref_init(&ran->ref);
+
+ mem = MEM_ALLOC(allocator, sizeof(normal_dist_float));
+ if(!mem) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ran->allocator = allocator;
+ ran->distrib = static_cast<normal_dist_float*>(new(mem) normal_dist_float);
+ ran->K1 = -1; /* invalid */
+ if (out_ran) *out_ran = ran;
+
+exit:
+ return res;
+error:
+ if(ran) {
+ SSP(ranst_gaussian_float_ref_put(ran));
+ ran = nullptr;
+ }
+ goto exit;
+}
+
+res_T
ssp_ranst_gaussian_setup
(struct ssp_ranst_gaussian* ran,
const double mu,
@@ -118,8 +184,26 @@ ssp_ranst_gaussian_setup
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);
+ ran->K1 = 1 / sigma;
+ ran->K2 = 1 / (sigma * SQRT_2_PI);
+
+ return RES_OK;
+}
+
+res_T
+ssp_ranst_gaussian_float_setup
+ (struct ssp_ranst_gaussian_float* ran,
+ const float mu,
+ const float sigma)
+{
+ if(!ran || sigma < 0)
+ return RES_BAD_ARG;
+
+ normal_dist_float::param_type p{ mu, sigma };
+ ran->distrib->param(p);
+ ran->mu = mu;
+ ran->K1 = 1 / sigma;
+ ran->K2 = 1 / (sigma * (float)SQRT_2_PI);
return RES_OK;
}
@@ -134,6 +218,15 @@ ssp_ranst_gaussian_ref_get
}
res_T
+ssp_ranst_gaussian_float_ref_get
+ (struct ssp_ranst_gaussian_float* 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)
{
@@ -142,6 +235,15 @@ ssp_ranst_gaussian_ref_put
return RES_OK;
}
+res_T
+ssp_ranst_gaussian_float_ref_put
+ (struct ssp_ranst_gaussian_float* ran)
+{
+ if(!ran) return RES_BAD_ARG;
+ ref_put(&ran->ref, gaussian_float_release);
+ return RES_OK;
+}
+
double
ssp_ranst_gaussian_get
(const struct ssp_ranst_gaussian* ran, struct ssp_rng* rng)
@@ -151,6 +253,15 @@ ssp_ranst_gaussian_get
return ran->distrib->operator()(r);
}
+float
+ssp_ranst_gaussian_float_get
+ (const struct ssp_ranst_gaussian_float* 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)
@@ -160,4 +271,12 @@ ssp_ranst_gaussian_pdf
return ran->K2 * exp(-0.5 * tmp * tmp);
}
+float
+ssp_ranst_gaussian_float_pdf
+ (const struct ssp_ranst_gaussian_float* ran, const float x)
+{
+ const float tmp = (x - ran->mu) * ran->K1;
+ ASSERT(ran->K1 > 0);
+ return ran->K2 * expf(-0.5f * tmp * tmp);
+}
diff --git a/src/ssp_ranst_piecewise_linear.c b/src/ssp_ranst_piecewise_linear.c
@@ -34,14 +34,24 @@
#include <rsys/mem_allocator.h>
#include <rsys/ref_count.h>
+#include <algorithm>
+
using piecewise_dist = RAN_NAMESPACE::piecewise_linear_distribution<double>;
+using piecewise_dist_float = RAN_NAMESPACE::piecewise_linear_distribution<float>;
+
struct ssp_ranst_piecewise_linear {
ref_T ref;
struct mem_allocator* allocator;
piecewise_dist* distrib;
};
+struct ssp_ranst_piecewise_linear_float {
+ ref_T ref;
+ struct mem_allocator* allocator;
+ piecewise_dist_float* distrib;
+};
+
/*******************************************************************************
* Helper function
******************************************************************************/
@@ -58,6 +68,19 @@ piecewise_release(ref_T* ref)
MEM_RM(ran->allocator, ran);
}
+static void
+piecewise_float_release(ref_T* ref)
+{
+ ssp_ranst_piecewise_linear_float* ran;
+ ASSERT(ref);
+ ran = CONTAINER_OF(ref, ssp_ranst_piecewise_linear_float, ref);
+ if(ran->distrib) {
+ ran->distrib->~piecewise_dist_float();
+ MEM_RM(ran->allocator, ran->distrib);
+ }
+ MEM_RM(ran->allocator, ran);
+}
+
/*******************************************************************************
* Exported functions
******************************************************************************/
@@ -105,21 +128,94 @@ error:
}
res_T
+ssp_ranst_piecewise_linear_float_create
+ (struct mem_allocator* allocator,
+ struct ssp_ranst_piecewise_linear_float** out_ran)
+{
+ struct ssp_ranst_piecewise_linear_float* 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_float*>
+ (MEM_CALLOC(allocator, 1, sizeof(ssp_ranst_piecewise_linear_float)));
+ if(!ran) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ref_init(&ran->ref);
+
+ mem = MEM_ALLOC(allocator, sizeof(piecewise_dist_float));
+ if(!mem) {
+ res = RES_MEM_ERR;
+ goto error;
+ }
+ ran->allocator = allocator;
+ ran->distrib
+ = static_cast<piecewise_dist_float*>(new(mem) piecewise_dist_float);
+ if (out_ran) *out_ran = ran;
+
+exit:
+ return res;
+error:
+ if(ran) {
+ SSP(ranst_piecewise_linear_float_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)
{
+ size_t i;
if(!ran || !intervals || !weights || size < 2)
return RES_BAD_ARG;
+ /* Checking param validity to avoid an assert when using ran */
+ for(i=0; i < size-1; i++) {
+ if(weights[i] < 0) return RES_BAD_ARG;
+ if(intervals[i+1] <= intervals[i]) return RES_BAD_ARG;
+ }
+ if (intervals[size-1] - intervals[size-2] <= 0) 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_float_setup
+ (struct ssp_ranst_piecewise_linear_float* ran,
+ const float* intervals,
+ const float* weights,
+ const size_t size)
+{
+ size_t i;
+ if(!ran || !intervals || !weights || size < 2)
+ return RES_BAD_ARG;
+
+ /* Checking param validity to avoid an assert when using ran */
+ for(i=0; i < size-1; i++) {
+ if(weights[i] < 0) return RES_BAD_ARG;
+ if(intervals[i+1] <= intervals[i]) return RES_BAD_ARG;
+ }
+ if (weights[size-1] < 0) return RES_BAD_ARG;
+ piecewise_dist_float::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)
{
@@ -129,6 +225,15 @@ ssp_ranst_piecewise_linear_ref_get
}
res_T
+ssp_ranst_piecewise_linear_float_ref_get
+ (struct ssp_ranst_piecewise_linear_float* 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)
{
@@ -137,6 +242,15 @@ ssp_ranst_piecewise_linear_ref_put
return RES_OK;
}
+res_T
+ssp_ranst_piecewise_linear_float_ref_put
+ (struct ssp_ranst_piecewise_linear_float* ran)
+{
+ if(!ran) return RES_BAD_ARG;
+ ref_put(&ran->ref, piecewise_float_release);
+ return RES_OK;
+}
+
double
ssp_ranst_piecewise_linear_get
(const struct ssp_ranst_piecewise_linear* ran,
@@ -146,3 +260,55 @@ ssp_ranst_piecewise_linear_get
return ran->distrib->operator()(r);
}
+float
+ssp_ranst_piecewise_linear_float_get
+ (const struct ssp_ranst_piecewise_linear_float* ran,
+ struct ssp_rng* rng)
+{
+ class rng_cxx r(*rng);
+ return ran->distrib->operator()(r);
+}
+
+double
+ssp_ranst_piecewise_linear_pdf
+ (const struct ssp_ranst_piecewise_linear *ran,
+ double x)
+{
+ ASSERT(ran);
+ if (x<ran->distrib->min() || x>ran->distrib->max())
+ return 0;
+
+ const auto& inter = ran->distrib->intervals();
+ const auto& dens = ran->distrib->densities();
+ auto b = std::lower_bound(inter.begin(), inter.end(), x);
+ size_t idx = b - inter.begin();
+ if (x == *b) return dens[idx];
+ idx--;
+ ASSERT(idx < inter.size() - 1);
+ return (dens[idx+1] * (x - inter[idx]) + dens[idx] * (inter[idx+1] - x))
+ / (inter[idx+1]- inter[idx]);
+}
+
+float
+ssp_ranst_piecewise_linear_float_pdf
+ (const struct ssp_ranst_piecewise_linear_float *ran,
+ float x)
+{
+ ASSERT(ran);
+ if (x<ran->distrib->min() || x>ran->distrib->max())
+ return 0;
+
+ const auto& inter = ran->distrib->intervals();
+ /* std::piecewise_linear_distribution<float>::densities()
+ * should be a std::vector<float>, but is a std::vector<double> with gcc
+ * => use explicit casts */
+ const auto& dens = ran->distrib->densities();
+ auto b = std::lower_bound(inter.begin(), inter.end(), x);
+ size_t idx = b - inter.begin();
+ if (x == *b) return (float)dens[idx];
+ idx--;
+ ASSERT(idx < inter.size() - 1);
+ return
+ ((float)dens[idx+1] * (x-inter[idx]) + (float)dens[idx] * (inter[idx+1]-x))
+ / (inter[idx+1] - inter[idx]);
+}
+\ No newline at end of file
diff --git a/src/ssp_rng.c b/src/ssp_rng.c
@@ -70,7 +70,7 @@
struct rng_kiss { uint32_t x, y, z, c; };
static res_T
-rng_kiss_set(void* data, const size_t seed)
+rng_kiss_set(void* data, const uint64_t seed)
{
struct rng_kiss* kiss = (struct rng_kiss*)data;
RAN_NAMESPACE::mt19937 rng_mt(seed % UINT32_MAX);
diff --git a/src/test_ssp_ran_discrete.h b/src/test_ssp_ran_discrete.h
@@ -0,0 +1,32 @@
+/* Copyright (C) |Meso|Star> 2015-2017 (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
+* 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. */
+
+
+/* This file is intentionally empty */
diff --git a/src/test_ssp_ran_gaussian.c b/src/test_ssp_ran_gaussian.c
@@ -28,91 +28,17 @@
* 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"
-#include <rsys/clock_time.h>
+#define TYPE_FLOAT 1
+#include "test_ssp_ran_gaussian.h"
-#define NBS 1000000
+#define TYPE_FLOAT 0
+#include "test_ssp_ran_gaussian.h"
int
main(int argc, char** argv)
{
- struct ssp_rng* rng;
- struct mem_allocator allocator;
- struct ssp_ranst_gaussian *gaussian;
- int i;
- double x = 0, x2 = 0;
- double exp_mean = 10, mean;
- double exp_std = 3, std;
- struct time start, end, res;
- char dump[512];
(void)argc, (void)argv;
-
- mem_init_proxy_allocator(&allocator, &mem_default_allocator);
-
- CHECK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng), RES_OK);
-
- CHECK(ssp_ranst_gaussian_create(NULL, NULL), RES_BAD_ARG);
- CHECK(ssp_ranst_gaussian_create(NULL, &gaussian), RES_OK);
- CHECK(ssp_ranst_gaussian_ref_put(gaussian), RES_OK);
-
- CHECK(ssp_ranst_gaussian_create(&allocator, NULL), RES_BAD_ARG);
- CHECK(ssp_ranst_gaussian_create(&allocator, &gaussian), RES_OK);
-
- CHECK(ssp_ranst_gaussian_setup(NULL, exp_mean, exp_std), RES_BAD_ARG);
- CHECK(ssp_ranst_gaussian_setup(gaussian, exp_mean, -1), RES_BAD_ARG);
- CHECK(ssp_ranst_gaussian_setup(gaussian, exp_mean, exp_std), RES_OK);
-
- CHECK(ssp_ranst_gaussian_ref_get(NULL), RES_BAD_ARG);
- CHECK(ssp_ranst_gaussian_ref_get(gaussian), RES_OK);
-
- CHECK(ssp_ranst_gaussian_ref_put(NULL), RES_BAD_ARG);
- CHECK(ssp_ranst_gaussian_ref_put(gaussian), RES_OK);
-
- time_current(&start);
- FOR_EACH(i, 0, NBS) {
- const double r = ssp_ranst_gaussian_get(gaussian, rng);
- x += r;
- x2 += r * r;
- }
- time_current(&end);
- time_sub(&res, &end, &start);
- time_dump(&res, TIME_SEC|TIME_MSEC|TIME_USEC, NULL, dump, sizeof(dump));
- printf("%s--\n", dump);
-
- mean = x/NBS;
- std = sqrt(x2/NBS - mean*mean);
- printf("%g %g\n", mean, std);
- CHECK(eq_eps(mean, exp_mean, 1e-2), 1);
- CHECK(eq_eps(std, exp_std, 1e-2), 1);
-
- /* Same test with inline gaussian generation */
- x = 0;
- x2 = 0;
-
- time_current(&start);
- FOR_EACH(i, 0, NBS) {
- const double r = ssp_ran_gaussian(rng, exp_mean, exp_std);
- x += r;
- x2 += r * r;
- }
- time_current(&end);
- time_sub(&res, &end, &start);
- time_dump(&res, TIME_SEC|TIME_MSEC|TIME_USEC, NULL, dump, sizeof(dump));
- printf("%s--\n", dump);
-
- mean = x/NBS;
- std = sqrt(x2/NBS - mean*mean);
- printf("%g %g\n", mean, std);
- CHECK(eq_eps(mean, exp_mean, 1.e-2), 1);
- CHECK(eq_eps(std, exp_std, 1e-2), 1);
-
- CHECK(ssp_ranst_gaussian_ref_put(gaussian), RES_OK);
-
- CHECK(ssp_rng_ref_put(rng), RES_OK);
-
- check_memory_allocator(&allocator);
- mem_shutdown_proxy_allocator(&allocator);
- CHECK(mem_allocated_size(), 0);
+ test_float();
+ test_double();
return 0;
}
diff --git a/src/test_ssp_ran_gaussian.h b/src/test_ssp_ran_gaussian.h
@@ -0,0 +1,164 @@
+/* Copyright (C) |Meso|Star> 2015-2017 (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. */
+
+#ifndef TEST_SSP_RAN_GAUSSIAN_H
+#define TEST_SSP_RAN_GAUSSIAN_H
+
+#include "ssp.h"
+#include "test_ssp_utils.h"
+
+#include <string.h>
+#include <rsys/clock_time.h>
+
+#define NBS 1000000
+
+#endif /* TEST_SSP_RAN_GAUSSIAN_H */
+
+#if TYPE_FLOAT==0
+#define REAL double
+#define TEST test_double
+#define RAN_GAUSSIAN ssp_ran_gaussian
+#define RANST_GAUSSIAN ssp_ranst_gaussian
+#define RANST_GAUSSIAN_CREATE ssp_ranst_gaussian_create
+#define RANST_GAUSSIAN_REF_GET ssp_ranst_gaussian_ref_get
+#define RANST_GAUSSIAN_REF_PUT ssp_ranst_gaussian_ref_put
+#define RANST_GAUSSIAN_SETUP ssp_ranst_gaussian_setup
+#define RANST_GAUSSIAN_GET ssp_ranst_gaussian_get
+#define SQRT sqrt
+
+#elif TYPE_FLOAT==1
+#define REAL float
+#define TEST test_float
+#define RAN_GAUSSIAN ssp_ran_gaussian_float
+#define RANST_GAUSSIAN ssp_ranst_gaussian_float
+#define RANST_GAUSSIAN_CREATE ssp_ranst_gaussian_float_create
+#define RANST_GAUSSIAN_REF_GET ssp_ranst_gaussian_float_ref_get
+#define RANST_GAUSSIAN_REF_PUT ssp_ranst_gaussian_float_ref_put
+#define RANST_GAUSSIAN_SETUP ssp_ranst_gaussian_float_setup
+#define RANST_GAUSSIAN_GET ssp_ranst_gaussian_float_get
+#define SQRT sqrtf
+
+#else
+#error "TYPE_SUFFIX must be defined either 0 or 1"
+#endif
+
+static void
+TEST()
+{
+ struct ssp_rng* rng;
+ struct mem_allocator allocator;
+ struct RANST_GAUSSIAN *gaussian;
+ int i;
+ REAL x = 0, x2 = 0;
+ REAL exp_mean = 10, mean;
+ REAL exp_std = 3, std;
+ struct time start, end, res;
+ char dump[512];
+
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+
+ CHECK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng), RES_OK);
+
+ CHECK(RANST_GAUSSIAN_CREATE(NULL, NULL), RES_BAD_ARG);
+ CHECK(RANST_GAUSSIAN_CREATE(NULL, &gaussian), RES_OK);
+ CHECK(RANST_GAUSSIAN_REF_PUT(gaussian), RES_OK);
+
+ CHECK(RANST_GAUSSIAN_CREATE(&allocator, NULL), RES_BAD_ARG);
+ CHECK(RANST_GAUSSIAN_CREATE(&allocator, &gaussian), RES_OK);
+
+ CHECK(RANST_GAUSSIAN_SETUP(NULL, exp_mean, exp_std), RES_BAD_ARG);
+ CHECK(RANST_GAUSSIAN_SETUP(gaussian, exp_mean, -1), RES_BAD_ARG);
+ CHECK(RANST_GAUSSIAN_SETUP(gaussian, exp_mean, exp_std), RES_OK);
+
+ CHECK(RANST_GAUSSIAN_REF_GET(NULL), RES_BAD_ARG);
+ CHECK(RANST_GAUSSIAN_REF_GET(gaussian), RES_OK);
+
+ CHECK(RANST_GAUSSIAN_REF_PUT(NULL), RES_BAD_ARG);
+ CHECK(RANST_GAUSSIAN_REF_PUT(gaussian), RES_OK);
+
+ time_current(&start);
+ FOR_EACH(i, 0, NBS) {
+ const REAL r = RANST_GAUSSIAN_GET(gaussian, rng);
+ x += r;
+ x2 += r * r;
+ }
+ time_current(&end);
+ time_sub(&res, &end, &start);
+ time_dump(&res, TIME_SEC | TIME_MSEC | TIME_USEC, NULL, dump, sizeof(dump));
+ printf("%s--\n", dump);
+
+ mean = x / NBS;
+ std = SQRT(x2 / NBS - mean*mean);
+ printf("%g %g\n", mean, std);
+ CHECK(eq_eps(mean, exp_mean, 1e-2), 1);
+ CHECK(eq_eps(std, exp_std, 1e-2), 1);
+
+ /* Same test with inline gaussian generation */
+ x = 0;
+ x2 = 0;
+
+ time_current(&start);
+ FOR_EACH(i, 0, NBS) {
+ const REAL r = RAN_GAUSSIAN(rng, exp_mean, exp_std);
+ x += r;
+ x2 += r * r;
+ }
+ time_current(&end);
+ time_sub(&res, &end, &start);
+ time_dump(&res, TIME_SEC | TIME_MSEC | TIME_USEC, NULL, dump, sizeof(dump));
+ printf("%s--\n", dump);
+
+ mean = x / NBS;
+ std = SQRT(x2 / NBS - mean*mean);
+ printf("%g %g\n", mean, std);
+ CHECK(eq_eps(mean, exp_mean, 1.e-2), 1);
+ CHECK(eq_eps(std, exp_std, 1e-2), 1);
+
+ CHECK(RANST_GAUSSIAN_REF_PUT(gaussian), RES_OK);
+
+ CHECK(ssp_rng_ref_put(rng), RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHECK(mem_allocated_size(), 0);
+}
+
+#undef TEST
+#undef REAL
+#undef RAN_GAUSSIAN
+#undef RANST_GAUSSIAN
+#undef RANST_GAUSSIAN_CREATE
+#undef RANST_GAUSSIAN_REF_GET
+#undef RANST_GAUSSIAN_REF_PUT
+#undef RANST_GAUSSIAN_SETUP
+#undef RANST_GAUSSIAN_GET
+#undef SQRT
+#undef TYPE_FLOAT
diff --git a/src/test_ssp_ran_hemisphere.c b/src/test_ssp_ran_hemisphere.c
@@ -26,144 +26,18 @@
* 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 TYPE_FLOAT 1
+#include "test_ssp_ran_hemisphere.h"
-#include <rsys/double4.h>
-
-#define NSAMPS 128
+#define TYPE_FLOAT 0
+#include "test_ssp_ran_hemisphere.h"
int
main(int argc, char** argv)
{
- struct ssp_rng* rng0, *rng1;
- struct mem_allocator allocator;
- double samps0[NSAMPS][4];
- double samps1[NSAMPS][4];
- double samps2[NSAMPS][4];
- double samps3[NSAMPS][4];
- int i = 0, j = 0;
- double* f;
(void)argc, (void)argv;
-
- mem_init_proxy_allocator(&allocator, &mem_default_allocator);
-
- CHECK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng0), RES_OK);
- CHECK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng1), RES_OK);
-
- ssp_rng_set(rng0, 0);
- FOR_EACH(i, 0, NSAMPS) {
- double frame[9];
- double up[3] = {0.f, 0.f, 1.f};
- double xyz[3];
- uint64_t seed = ssp_rng_get(rng0);
-
- ssp_rng_set(rng1, seed);
- f = ssp_ran_hemisphere_uniform_local(rng1, samps0[i]);
- CHECK(f, samps0[i]);
- CHECK(d3_is_normalized(f), 1);
- CHECK(eq_eps(f[3], (float)(1.0/(2.0*PI)), 1.e-6), 1);
- CHECK(eq_eps(f[3], ssp_ran_hemisphere_uniform_local_pdf(f), 1.e-6), 1);
- CHECK(d3_dot(f, up) >= 0.f, 1);
-
- ssp_rng_set(rng1, seed);
- f = ssp_ran_hemisphere_uniform(rng1, up, samps1[i]);
- CHECK(f, samps1[i]);
- CHECK(d4_eq_eps(f, samps0[i], 1.e-6), 1);
-
- up[0] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
- up[1] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
- up[2] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
- d3_normalize(up, up);
-
- ssp_rng_set(rng1, seed);
- f = ssp_ran_hemisphere_uniform(rng1, up, samps1[i]);
- CHECK(d3_eq_eps(samps0[i], samps1[i], 1.e-6), 0);
- CHECK(d3_is_normalized(f), 1);
- CHECK(d3_dot(f, up) >= 0.f, 1);
- CHECK(eq_eps(f[3], 1.0/(2.0*PI), 1.e-6 ), 1);
- CHECK(eq_eps(f[3], ssp_ran_hemisphere_uniform_pdf(up, f), 1.e-6), 1);
-
- d33_basis(frame, up);
- d33_muld3(xyz, frame, samps0[i]);
- CHECK(d3_eq_eps(samps1[i], xyz, 1.e-6), 1);
- FOR_EACH(j, 0, i) {
- CHECK(d3_eq_eps(samps0[i], samps0[j], 1.e-6), 0);
- CHECK(d3_eq_eps(samps1[i], samps1[j], 1.e-6), 0);
- }
- }
-
- samps1[1][0] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
- samps1[1][1] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
- samps1[1][2] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
- d3_normalize(samps1[1], samps1[1]);
-
- ssp_rng_set(rng0, 0);
- ssp_ran_hemisphere_uniform(rng0, samps1[1], samps0[0]);
- ssp_rng_set(rng0, 0);
- ssp_ran_hemisphere_uniform(rng0, samps1[1], samps1[1]);
- CHECK(d4_eq(samps0[0], samps1[1]), 1);
-
- ssp_rng_set(rng0, 0);
- FOR_EACH(i, 0, NSAMPS) {
- double frame[9];
- double up[3] = { 0.0, 0.0, 1.0 };
- double xyz[3];
- uint64_t seed = ssp_rng_get(rng0);
-
- ssp_rng_set(rng1, seed);
- f = ssp_ran_hemisphere_cos_local(rng1, samps2[i]);
- CHECK(f, samps2[i]);
- CHECK(d3_eq_eps(samps0[i], samps2[i], 1.e-6), 0);
- CHECK(d3_is_normalized(f), 1);
- CHECK(eq_eps(f[3], f[2]/PI, 1.e-6), 1);
- CHECK(eq_eps(f[3], ssp_ran_hemisphere_cos_local_pdf(f), 1.e-6), 1);
- CHECK(d3_dot(f, up) >= 0.f, 1);
-
- ssp_rng_set(rng1, seed);
- f = ssp_ran_hemisphere_cos(rng1, up, samps3[i]);
- CHECK(f, samps3[i]);
- CHECK(d4_eq_eps(f, samps2[i], 1.e-6), 1);
-
- up[0] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
- up[1] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
- up[2] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
- d3_normalize(up, up);
-
- ssp_rng_set(rng1, seed);
- f = ssp_ran_hemisphere_cos(rng1, up, samps3[i]);
- CHECK(d3_eq_eps(samps2[i], samps3[i], 1.e-6), 0);
- CHECK(d3_eq_eps(samps1[i], samps3[i], 1.e-6), 0);
- CHECK(d3_is_normalized(f), 1);
- CHECK(d3_dot(f, up) >= 0.f, 1);
- CHECK(eq_eps(f[3], d3_dot(f, up)/PI, 1.e-6), 1);
- CHECK(eq_eps(f[3], ssp_ran_hemisphere_cos_pdf(f, up), 1.e-6), 1);
-
- d33_basis(frame, up);
- d33_muld3(xyz, frame, samps2[i]);
- CHECK(d3_eq_eps(samps3[i], xyz, 1.e-6), 1);
- FOR_EACH(j, 0, i) {
- CHECK(d3_eq_eps(samps2[i], samps2[j], 1.e-6), 0);
- CHECK(d3_eq_eps(samps3[i], samps3[j], 1.e-6), 0);
- }
- }
-
- samps1[1][0] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
- samps1[1][1] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
- samps1[1][2] = ssp_rng_uniform_double(rng1, -1.0, 1.0);
- d3_normalize(samps1[1], samps1[1]);
-
- ssp_rng_set(rng0, 0);
- ssp_ran_hemisphere_cos(rng0, samps1[1], samps0[0]);
- ssp_rng_set(rng0, 0);
- ssp_ran_hemisphere_cos(rng0, samps1[1], samps1[1]);
- CHECK(d4_eq(samps0[0], samps1[1]), 1);
-
- ssp_rng_ref_put(rng0);
- ssp_rng_ref_put(rng1);
- check_memory_allocator(&allocator);
- mem_shutdown_proxy_allocator(&allocator);
-
- CHECK(mem_allocated_size(), 0);
+ test_float();
+ test_double();
return 0;
}
+
diff --git a/src/test_ssp_ran_hemisphere.h b/src/test_ssp_ran_hemisphere.h
@@ -0,0 +1,247 @@
+/* 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_HEMISPHERE_H
+#define TEST_SSP_RAN_HEMISPHERE_H
+
+#include "ssp.h"
+#include "test_ssp_utils.h"
+
+#include <rsys/double4.h>
+#include <rsys/float4.h>
+
+#define NSAMPS 128
+
+#endif /* TEST_SSP_RAN_HEMISPHERE_H */
+
+#if TYPE_FLOAT==0
+#define REAL double
+#define TEST test_double
+#define RNG_UNIFORM_R ssp_rng_uniform_double
+#define RAN_HEMISPHERE_UNIFORM_LOCAL ssp_ran_hemisphere_uniform_local
+#define RAN_HEMISPHERE_UNIFORM_LOCAL_PDF ssp_ran_hemisphere_uniform_local_pdf
+#define RAN_HEMISPHERE_COS_LOCAL ssp_ran_hemisphere_cos_local
+#define RAN_HEMISPHERE_COS_LOCAL_PDF ssp_ran_hemisphere_cos_local_pdf
+#define RAN_HEMISPHERE_UNIFORM ssp_ran_hemisphere_uniform
+#define RAN_HEMISPHERE_COS ssp_ran_hemisphere_cos
+#define RAN_HEMISPHERE_COS_PDF ssp_ran_hemisphere_cos_pdf
+#define RAN_HEMISPHERE_UNIFORM_PDF ssp_ran_hemisphere_uniform_pdf
+#define R3_NORMALIZE d3_normalize
+#define R3_IS_NORMALIZED d3_is_normalized
+#define R3_DOT d3_dot
+#define R4_EQ_EPS d4_eq_eps
+#define R4_EQ d4_eq
+#define R3_EQ_EPS d3_eq_eps
+#define EQ_EPS_R eq_eps
+#define R33_BASIS d33_basis
+#define R33_MULR3 d33_muld3
+
+#elif TYPE_FLOAT==1
+#define REAL float
+#define TEST test_float
+#define RNG_UNIFORM_R ssp_rng_uniform_float
+#define RAN_HEMISPHERE_UNIFORM_LOCAL ssp_ran_hemisphere_uniform_float_local
+#define RAN_HEMISPHERE_UNIFORM_LOCAL_PDF ssp_ran_hemisphere_uniform_float_local_pdf
+#define RAN_HEMISPHERE_COS_LOCAL ssp_ran_hemisphere_cos_float_local
+#define RAN_HEMISPHERE_COS_LOCAL_PDF ssp_ran_hemisphere_cos_float_local_pdf
+#define RAN_HEMISPHERE_UNIFORM ssp_ran_hemisphere_uniform_float
+#define RAN_HEMISPHERE_COS ssp_ran_hemisphere_cos_float
+#define RAN_HEMISPHERE_COS_PDF ssp_ran_hemisphere_cos_float_pdf
+#define RAN_HEMISPHERE_UNIFORM_PDF ssp_ran_hemisphere_uniform_float_pdf
+#define R3_NORMALIZE f3_normalize
+#define R3_IS_NORMALIZED f3_is_normalized
+#define R3_DOT f3_dot
+#define R4_EQ_EPS f4_eq_eps
+#define R4_EQ f4_eq
+#define R3_EQ_EPS f3_eq_eps
+#define EQ_EPS_R eq_eps
+#define R33_BASIS f33_basis
+#define R33_MULR3 f33_mulf3
+
+#else
+#error "TYPE_SUFFIX must be defined either 0 or 1"
+#endif
+
+static void
+TEST()
+{
+ struct ssp_rng* rng0, *rng1;
+ struct mem_allocator allocator;
+ REAL samps0[NSAMPS][4];
+ REAL samps1[NSAMPS][4];
+ REAL samps2[NSAMPS][4];
+ REAL samps3[NSAMPS][4];
+ int i = 0, j = 0;
+ REAL* f;
+
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+
+ CHECK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng0), RES_OK);
+ CHECK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng1), RES_OK);
+
+ samps0[0][0] = 0; samps0[0][1] = 0; samps0[0][2] = 1;
+ f = RAN_HEMISPHERE_UNIFORM(rng1, samps0[0], samps1[0], NULL);
+ f = RAN_HEMISPHERE_UNIFORM(rng1, samps0[0], samps2[0], &samps2[0][3]);
+
+ ssp_rng_set(rng0, 0);
+ FOR_EACH(i, 0, NSAMPS) {
+ REAL frame[9];
+ REAL up[3] = {0, 0, 1};
+ REAL xyz[3];
+ uint64_t seed = ssp_rng_get(rng0);
+
+ ssp_rng_set(rng1, seed);
+ f = RAN_HEMISPHERE_UNIFORM_LOCAL(rng1, samps0[i], &samps0[i][3]);
+ CHECK(f, samps0[i]);
+ CHECK(R3_IS_NORMALIZED(f), 1);
+ CHECK(EQ_EPS_R(f[3], (1/(2*(REAL)PI)), (REAL)1.e-6), 1);
+ CHECK(EQ_EPS_R(f[3], RAN_HEMISPHERE_UNIFORM_LOCAL_PDF(f), (REAL)1.e-6), 1);
+ CHECK(R3_DOT(f, up) >= 0, 1);
+
+ ssp_rng_set(rng1, seed);
+ f = RAN_HEMISPHERE_UNIFORM(rng1, up, samps1[i], &samps1[i][3]);
+ CHECK(f, samps1[i]);
+ CHECK(R4_EQ_EPS(f, samps0[i], (REAL)1.e-6), 1);
+
+ up[0] = RNG_UNIFORM_R(rng1, -1, 1);
+ up[1] = RNG_UNIFORM_R(rng1, -1, 1);
+ up[2] = RNG_UNIFORM_R(rng1, -1, 1);
+ R3_NORMALIZE(up, up);
+
+ ssp_rng_set(rng1, seed);
+ f = RAN_HEMISPHERE_UNIFORM(rng1, up, samps1[i], &samps1[i][3]);
+ CHECK(R3_EQ_EPS(samps0[i], samps1[i], (REAL)1.e-6), 0);
+ CHECK(R3_IS_NORMALIZED(f), 1);
+ CHECK(R3_DOT(f, up) >= 0, 1);
+ CHECK(EQ_EPS_R(f[3], 1/(2*(REAL)PI), (REAL)1.e-6 ), 1);
+ CHECK(EQ_EPS_R(f[3], RAN_HEMISPHERE_UNIFORM_PDF(up, f), (REAL)1.e-6), 1);
+
+ R33_BASIS(frame, up);
+ R33_MULR3(xyz, frame, samps0[i]);
+ CHECK(R3_EQ_EPS(samps1[i], xyz, (REAL)1.e-6), 1);
+ FOR_EACH(j, 0, i) {
+ CHECK(R3_EQ_EPS(samps0[i], samps0[j], (REAL)1.e-6), 0);
+ CHECK(R3_EQ_EPS(samps1[i], samps1[j], (REAL)1.e-6), 0);
+ }
+ }
+
+ samps1[1][0] = RNG_UNIFORM_R(rng1, -1, 1);
+ samps1[1][1] = RNG_UNIFORM_R(rng1, -1, 1);
+ samps1[1][2] = RNG_UNIFORM_R(rng1, -1, 1);
+ R3_NORMALIZE(samps1[1], samps1[1]);
+
+ ssp_rng_set(rng0, 0);
+ RAN_HEMISPHERE_UNIFORM(rng0, samps1[1], samps0[0], &samps0[0][3]);
+ ssp_rng_set(rng0, 0);
+ RAN_HEMISPHERE_UNIFORM(rng0, samps1[1], samps1[1], &samps1[1][3]);
+ CHECK(R4_EQ(samps0[0], samps1[1]), 1);
+
+ ssp_rng_set(rng0, 0);
+ FOR_EACH(i, 0, NSAMPS) {
+ REAL frame[9];
+ REAL up[3] = { 0, 0, 1 };
+ REAL xyz[3];
+ uint64_t seed = ssp_rng_get(rng0);
+
+ ssp_rng_set(rng1, seed);
+ f = RAN_HEMISPHERE_COS_LOCAL(rng1, samps2[i], &samps2[i][3]);
+ CHECK(f, samps2[i]);
+ CHECK(R3_EQ_EPS(samps0[i], samps2[i], (REAL)1.e-6), 0);
+ CHECK(R3_IS_NORMALIZED(f), 1);
+ CHECK(EQ_EPS_R(f[3], f[2]/(REAL)PI, (REAL)1.e-6), 1);
+ CHECK(EQ_EPS_R(f[3], RAN_HEMISPHERE_COS_LOCAL_PDF(f), (REAL)1.e-6), 1);
+ CHECK(R3_DOT(f, up) >= 0, 1);
+
+ ssp_rng_set(rng1, seed);
+ f = RAN_HEMISPHERE_COS(rng1, up, samps3[i], &samps3[i][3]);
+ CHECK(f, samps3[i]);
+ CHECK(R4_EQ_EPS(f, samps2[i], (REAL)1.e-6), 1);
+
+ up[0] = RNG_UNIFORM_R(rng1, -1, 1);
+ up[1] = RNG_UNIFORM_R(rng1, -1, 1);
+ up[2] = RNG_UNIFORM_R(rng1, -1, 1);
+ R3_NORMALIZE(up, up);
+
+ ssp_rng_set(rng1, seed);
+ f = RAN_HEMISPHERE_COS(rng1, up, samps3[i], &samps3[i][3]);
+ CHECK(R3_EQ_EPS(samps2[i], samps3[i], (REAL)1.e-6), 0);
+ CHECK(R3_EQ_EPS(samps1[i], samps3[i], (REAL)1.e-6), 0);
+ CHECK(R3_IS_NORMALIZED(f), 1);
+ CHECK(R3_DOT(f, up) >= 0.f, 1);
+ CHECK(EQ_EPS_R(f[3], R3_DOT(f, up)/PI, (REAL)1.e-6), 1);
+ CHECK(EQ_EPS_R(f[3], RAN_HEMISPHERE_COS_PDF(f, up), (REAL)1.e-6), 1);
+
+ R33_BASIS(frame, up);
+ R33_MULR3(xyz, frame, samps2[i]);
+ CHECK(R3_EQ_EPS(samps3[i], xyz, (REAL)1.e-6), 1);
+ FOR_EACH(j, 0, i) {
+ CHECK(R3_EQ_EPS(samps2[i], samps2[j], (REAL)1.e-6), 0);
+ CHECK(R3_EQ_EPS(samps3[i], samps3[j], (REAL)1.e-6), 0);
+ }
+ }
+
+ samps1[1][0] = RNG_UNIFORM_R(rng1, -1, 1);
+ samps1[1][1] = RNG_UNIFORM_R(rng1, -1, 1);
+ samps1[1][2] = RNG_UNIFORM_R(rng1, -1, 1);
+ R3_NORMALIZE(samps1[1], samps1[1]);
+
+ ssp_rng_set(rng0, 0);
+ RAN_HEMISPHERE_COS(rng0, samps1[1], samps0[0], &samps0[0][3]);
+ ssp_rng_set(rng0, 0);
+ RAN_HEMISPHERE_COS(rng0, samps1[1], samps1[1], &samps1[1][3]);
+ CHECK(R4_EQ(samps0[0], samps1[1]), 1);
+
+ ssp_rng_ref_put(rng0);
+ ssp_rng_ref_put(rng1);
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+
+ CHECK(mem_allocated_size(), 0);
+}
+
+#undef REAL
+#undef TEST
+#undef RNG_UNIFORM_R
+#undef RAN_HEMISPHERE_UNIFORM_LOCAL
+#undef RAN_HEMISPHERE_UNIFORM_LOCAL_PDF
+#undef RAN_HEMISPHERE_COS_LOCAL
+#undef RAN_HEMISPHERE_COS_LOCAL_PDF
+#undef RAN_HEMISPHERE_UNIFORM
+#undef RAN_HEMISPHERE_COS
+#undef RAN_HEMISPHERE_COS_PDF
+#undef RAN_HEMISPHERE_UNIFORM_PDF
+#undef R3_NORMALIZE
+#undef R3_IS_NORMALIZED
+#undef R3_DOT
+#undef R4_EQ_EPS
+#undef R4_EQ
+#undef R3_EQ_EPS
+#undef EQ_EPS_R
+#undef R33_BASIS
+#undef R33_MULR3
+#undef TYPE_FLOAT
diff --git a/src/test_ssp_ran_hg.c b/src/test_ssp_ran_hg.c
@@ -26,50 +26,17 @@
* 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 TYPE_FLOAT 1
+#include "test_ssp_ran_hg.h"
-#define NG 100
-#define NSAMPS 10000
+#define TYPE_FLOAT 0
+#include "test_ssp_ran_hg.h"
int
main(int argc, char** argv)
{
- struct ssp_rng* rng;
- struct mem_allocator allocator;
- int i, j;
(void)argc, (void)argv;
-
- mem_init_proxy_allocator(&allocator, &mem_default_allocator);
-
- CHECK(ssp_rng_create(&allocator, &ssp_rng_threefry, &rng), RES_OK);
-
- FOR_EACH(i, 0, NG) {
- /* for any value of g... */
- double g = ssp_rng_uniform_double(rng, -1, +1);
- double sum_cos = 0;
- double sum_cos_local = 0;
- double dir[4], up[4] = {0, 0, 1};
- FOR_EACH(j, 0, NSAMPS) {
- /* HG relative to the Z axis */
- ssp_ran_sphere_hg_local(rng, g, dir);
- sum_cos_local += d3_dot(up, dir);
- }
- FOR_EACH(j, 0, NSAMPS) {
- /* HG relative to a up uniformaly sampled */
- ssp_ran_hemisphere_uniform_local(rng, up);
- ssp_ran_sphere_hg(rng, up, g, dir);
- sum_cos += d3_dot(up, dir);
- }
- /* ...on average cos(up, dir) should be g */
- CHECK(eq_eps(sum_cos_local / NSAMPS, g, 2.5e-2), 1);
- CHECK(eq_eps(sum_cos / NSAMPS, g, 2.5e-2), 1);
- }
-
- ssp_rng_ref_put(rng);
- check_memory_allocator(&allocator);
- mem_shutdown_proxy_allocator(&allocator);
-
- CHECK(mem_allocated_size(), 0);
+ test_float();
+ test_double();
return 0;
}
diff --git a/src/test_ssp_ran_hg.h b/src/test_ssp_ran_hg.h
@@ -0,0 +1,111 @@
+/* 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_HG_H
+#define TEST_SSP_RAN_HG_H
+
+#include "ssp.h"
+#include "test_ssp_utils.h"
+
+#define NG 100
+#define NSAMPS 10000
+
+#endif /* TEST_SSP_RAN_HG_H */
+
+#if TYPE_FLOAT==0
+#define REAL double
+#define TEST test_double
+#define RNG_UNIFORM_R ssp_rng_uniform_double
+#define RAN_SPHERE_HG ssp_ran_sphere_hg
+#define RAN_SPHERE_HG_LOCAL ssp_ran_sphere_hg_local
+#define RAN_HEMISPHERE_UNIFORM_LOCAL ssp_ran_hemisphere_uniform_local
+#define EQ_EPS_R eq_eps
+#define R3_DOT d3_dot
+
+#elif TYPE_FLOAT==1
+#define REAL float
+#define TEST test_float
+#define RNG_UNIFORM_R ssp_rng_uniform_float
+#define RAN_SPHERE_HG ssp_ran_sphere_hg_float
+#define RAN_SPHERE_HG_LOCAL ssp_ran_sphere_hg_float_local
+#define RAN_HEMISPHERE_UNIFORM_LOCAL ssp_ran_hemisphere_uniform_float_local
+#define EQ_EPS_R eq_epsf
+#define R3_DOT f3_dot
+#else
+#error "TYPE_SUFFIX must be defined either 0 or 1"
+#endif
+
+static void
+TEST()
+{
+ struct ssp_rng* rng;
+ struct mem_allocator allocator;
+ int i, j;
+
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+
+ CHECK(ssp_rng_create(&allocator, &ssp_rng_threefry, &rng), RES_OK);
+
+ FOR_EACH(i, 0, NG) {
+ /* for any value of g... */
+ REAL g = RNG_UNIFORM_R(rng, -1, +1);
+ REAL sum_cos = 0;
+ REAL sum_cos_local = 0;
+ REAL dir[3], pdf, up[4] = {0, 0, 1};
+ FOR_EACH(j, 0, NSAMPS) {
+ /* HG relative to the Z axis */
+ RAN_SPHERE_HG_LOCAL(rng, g, dir, &pdf);
+ sum_cos_local += R3_DOT(up, dir);
+ }
+ FOR_EACH(j, 0, NSAMPS) {
+ /* HG relative to a up uniformaly sampled */
+ RAN_HEMISPHERE_UNIFORM_LOCAL(rng, up, &pdf);
+ RAN_SPHERE_HG(rng, up, g, dir, &pdf);
+ sum_cos += R3_DOT(up, dir);
+ }
+ /* ...on average cos(up, dir) should be g */
+ CHECK(EQ_EPS_R(sum_cos_local / NSAMPS, g, (REAL)2.5e-2), 1);
+ CHECK(EQ_EPS_R(sum_cos / NSAMPS, g, (REAL)2.5e-2), 1);
+ }
+
+ ssp_rng_ref_put(rng);
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+
+ CHECK(mem_allocated_size(), 0);
+}
+
+#undef REAL
+#undef TEST
+#undef RNG_UNIFORM_R
+#undef RAN_SPHERE_HG
+#undef RAN_SPHERE_HG_LOCAL
+#undef RAN_HEMISPHERE_UNIFORM_LOCAL
+#undef EQ_EPS_R
+#undef R3_DOT
+#undef TYPE_FLOAT
diff --git a/src/test_ssp_ran_piecewise_linear.c b/src/test_ssp_ran_piecewise_linear.c
@@ -28,72 +28,17 @@
* 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 TYPE_FLOAT 1
+#include "test_ssp_ran_piecewise_linear.h"
-#define NBS 1000000
+#define TYPE_FLOAT 0
+#include "test_ssp_ran_piecewise_linear.h"
int
main(int argc, char** argv)
{
- struct ssp_rng* rng;
- struct mem_allocator allocator;
- struct ssp_ranst_piecewise_linear *pwl;
- int i;
- double exp_mean = 5, mean;
- double exp_std = 10 / sqrt(12) /*sqrt((b - a)² / 12) */, std;
- double x = 0, x2 = 0;
- const double intervals[] = { 0, 1, 3, 5, 7, 8, 10 };
- const double weights[] = { 1, 1, 1, 1, 1, 1, 1 };
(void)argc, (void)argv;
-
- mem_init_proxy_allocator(&allocator, &mem_default_allocator);
-
- CHECK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng), RES_OK);
-
- CHECK(ssp_ranst_piecewise_linear_create(NULL, NULL), RES_BAD_ARG);
- CHECK(ssp_ranst_piecewise_linear_create(NULL, &pwl), RES_OK);
- CHECK(ssp_ranst_piecewise_linear_ref_put(pwl), RES_OK);
-
- CHECK(ssp_ranst_piecewise_linear_create(&allocator, NULL), RES_BAD_ARG);
- CHECK(ssp_ranst_piecewise_linear_create(&allocator, &pwl), RES_OK);
-
- CHECK(ssp_ranst_piecewise_linear_setup
- (NULL, intervals, weights, sizeof(intervals)/sizeof(double)), RES_BAD_ARG);
- CHECK(ssp_ranst_piecewise_linear_setup
- (pwl, NULL, weights, sizeof(intervals)/sizeof(double)), RES_BAD_ARG);
- CHECK(ssp_ranst_piecewise_linear_setup
- (pwl, intervals, NULL, sizeof(intervals)/sizeof(double)), RES_BAD_ARG);
- CHECK(ssp_ranst_piecewise_linear_setup
- (pwl, intervals, weights, 1), RES_BAD_ARG);
- CHECK(ssp_ranst_piecewise_linear_setup
- (pwl, intervals, weights, sizeof(intervals)/sizeof(double)), RES_OK);
-
- CHECK(ssp_ranst_piecewise_linear_ref_get(NULL), RES_BAD_ARG);
- CHECK(ssp_ranst_piecewise_linear_ref_get(pwl), RES_OK);
-
- CHECK(ssp_ranst_piecewise_linear_ref_put(NULL), RES_BAD_ARG);
- CHECK(ssp_ranst_piecewise_linear_ref_put(pwl), RES_OK);
-
- FOR_EACH(i, 0, NBS) {
- double r = ssp_ranst_piecewise_linear_get(pwl, rng);
- CHECK(0 <= r && r <= 10, 1);
- x += r;
- x2 += r * r;
- }
-
- mean = x/NBS;
- std = sqrt(x2/NBS - mean*mean);
- printf("%g %g\n", mean, std);
- CHECK(eq_eps(mean, exp_mean, 1e-2), 1);
- CHECK(eq_eps(std, exp_std, 1.e-2), 1);
-
- CHECK(ssp_ranst_piecewise_linear_ref_put(pwl), RES_OK);
-
- CHECK(ssp_rng_ref_put(rng), RES_OK);
-
- check_memory_allocator(&allocator);
- mem_shutdown_proxy_allocator(&allocator);
- CHECK(mem_allocated_size(), 0);
+ test_float();
+ test_double();
return 0;
}
diff --git a/src/test_ssp_ran_piecewise_linear.h b/src/test_ssp_ran_piecewise_linear.h
@@ -0,0 +1,174 @@
+/* Copyright (C) |Meso|Star> 2015-2017 (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
+ * 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_PIECEWISE_LINEAR_H
+#define TEST_SSP_RAN_PIECEWISE_LINEAR_H
+
+#include "ssp.h"
+#include "test_ssp_utils.h"
+
+#define NBS 1000000
+
+#endif /* TEST_SSP_RAN_PIECEWISE_LINEAR_H */
+
+#if TYPE_FLOAT==0
+#define REAL double
+#define TEST test_double
+#define RANST_PIECEWISE_LINEAR ssp_ranst_piecewise_linear
+#define RANST_PIECEWISE_LINEAR_CREATE ssp_ranst_piecewise_linear_create
+#define RANST_PIECEWISE_LINEAR_SETUP ssp_ranst_piecewise_linear_setup
+#define RANST_PIECEWISE_LINEAR_GET ssp_ranst_piecewise_linear_get
+#define RANST_PIECEWISE_LINEAR_PDF ssp_ranst_piecewise_linear_pdf
+#define RANST_PIECEWISE_LINEAR_REF_GET ssp_ranst_piecewise_linear_ref_get
+#define RANST_PIECEWISE_LINEAR_REF_PUT ssp_ranst_piecewise_linear_ref_put
+#define EQ_EPS_R eq_eps
+#define EPS_R DBL_EPSILON
+#define SQRT sqrt
+
+#elif TYPE_FLOAT==1
+#define REAL float
+#define TEST test_float
+#define RANST_PIECEWISE_LINEAR ssp_ranst_piecewise_linear_float
+#define RANST_PIECEWISE_LINEAR_CREATE ssp_ranst_piecewise_linear_float_create
+#define RANST_PIECEWISE_LINEAR_SETUP ssp_ranst_piecewise_linear_float_setup
+#define RANST_PIECEWISE_LINEAR_GET ssp_ranst_piecewise_linear_float_get
+#define RANST_PIECEWISE_LINEAR_PDF ssp_ranst_piecewise_linear_float_pdf
+#define RANST_PIECEWISE_LINEAR_REF_GET ssp_ranst_piecewise_linear_float_ref_get
+#define RANST_PIECEWISE_LINEAR_REF_PUT ssp_ranst_piecewise_linear_float_ref_put
+#define EQ_EPS_R eq_epsf
+#define EPS_R FLT_EPSILON
+#define SQRT sqrtf
+#else
+#error "TYPE_SUFFIX must be defined either 0 or 1"
+#endif
+
+static void
+TEST()
+{
+ struct ssp_rng* rng;
+ struct mem_allocator allocator;
+ struct RANST_PIECEWISE_LINEAR *pwl;
+ int i;
+ REAL exp_mean = 5, mean;
+ REAL exp_std = 10 / SQRT(12) /*sqrt((b - a)² / 12) */, std;
+ REAL x = 0, x2 = 0;
+ REAL intervals[] = { 0, 1, 3, 5, 7, 8, 10 };
+ REAL weights[] = { 1, 1, 1, 1, 1, 1, 1 };
+
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+
+ CHECK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng), RES_OK);
+
+ CHECK(RANST_PIECEWISE_LINEAR_CREATE(NULL, NULL), RES_BAD_ARG);
+ CHECK(RANST_PIECEWISE_LINEAR_CREATE(NULL, &pwl), RES_OK);
+ CHECK(RANST_PIECEWISE_LINEAR_REF_PUT(pwl), RES_OK);
+
+ CHECK(RANST_PIECEWISE_LINEAR_CREATE(&allocator, NULL), RES_BAD_ARG);
+ CHECK(RANST_PIECEWISE_LINEAR_CREATE(&allocator, &pwl), RES_OK);
+
+ CHECK(RANST_PIECEWISE_LINEAR_SETUP
+ (NULL, intervals, weights, sizeof(intervals)/sizeof(REAL)), RES_BAD_ARG);
+ CHECK(RANST_PIECEWISE_LINEAR_SETUP
+ (pwl, NULL, weights, sizeof(intervals)/sizeof(REAL)), RES_BAD_ARG);
+ CHECK(RANST_PIECEWISE_LINEAR_SETUP
+ (pwl, intervals, NULL, sizeof(intervals)/sizeof(REAL)), RES_BAD_ARG);
+ CHECK(RANST_PIECEWISE_LINEAR_SETUP
+ (pwl, intervals, weights, 1), RES_BAD_ARG);
+ CHECK(RANST_PIECEWISE_LINEAR_SETUP
+ (pwl, intervals, weights, sizeof(intervals)/sizeof(REAL)), RES_OK);
+
+ CHECK(RANST_PIECEWISE_LINEAR_REF_PUT(pwl), RES_OK);
+ CHECK(RANST_PIECEWISE_LINEAR_CREATE(&allocator, &pwl), RES_OK);
+
+ weights[1] = -1;
+ CHECK(RANST_PIECEWISE_LINEAR_SETUP
+ (pwl, intervals, weights, sizeof(intervals) / sizeof(REAL)), RES_BAD_ARG);
+ weights[1] = 1;
+
+ intervals[1] = 4;
+ CHECK(RANST_PIECEWISE_LINEAR_SETUP
+ (pwl, intervals, weights, sizeof(intervals) / sizeof(REAL)), RES_BAD_ARG);
+ intervals[1] = 1;
+
+ intervals[1] = 3;
+ CHECK(RANST_PIECEWISE_LINEAR_SETUP
+ (pwl, intervals, weights, sizeof(intervals) / sizeof(REAL)), RES_BAD_ARG);
+ intervals[1] = 1;
+
+ CHECK(RANST_PIECEWISE_LINEAR_SETUP
+ (pwl, intervals, weights, sizeof(intervals) / sizeof(REAL)), RES_OK);
+
+ CHECK(RANST_PIECEWISE_LINEAR_REF_GET(NULL), RES_BAD_ARG);
+ CHECK(RANST_PIECEWISE_LINEAR_REF_GET(pwl), RES_OK);
+
+ CHECK(RANST_PIECEWISE_LINEAR_REF_PUT(NULL), RES_BAD_ARG);
+ CHECK(RANST_PIECEWISE_LINEAR_REF_PUT(pwl), RES_OK);
+
+ FOR_EACH(i, 0, NBS) {
+ REAL r = RANST_PIECEWISE_LINEAR_GET(pwl, rng);
+ CHECK(0 <= r && r <= 10, 1);
+ REAL pdf = RANST_PIECEWISE_LINEAR_PDF(pwl, r);
+ CHECK(EQ_EPS_R(pdf, (REAL)0.1, EPS_R), 1);
+ x += r;
+ x2 += r * r;
+ }
+ CHECK(EQ_EPS_R(RANST_PIECEWISE_LINEAR_PDF(pwl, 0), (REAL)0.1, EPS_R), 1);
+ CHECK(EQ_EPS_R(RANST_PIECEWISE_LINEAR_PDF(pwl, 10), (REAL)0.1, EPS_R), 1);
+ CHECK(RANST_PIECEWISE_LINEAR_PDF(pwl, -1), 0);
+ CHECK(RANST_PIECEWISE_LINEAR_PDF(pwl, 11), 0);
+
+ mean = x/NBS;
+ std = SQRT(x2/NBS - mean*mean);
+ printf("%g %g\n", mean, std);
+ CHECK(EQ_EPS_R(mean, exp_mean, (REAL)1e-2), 1);
+ CHECK(EQ_EPS_R(std, exp_std, (REAL)1.e-2), 1);
+
+ CHECK(RANST_PIECEWISE_LINEAR_REF_PUT(pwl), RES_OK);
+
+ CHECK(ssp_rng_ref_put(rng), RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHECK(mem_allocated_size(), 0);
+}
+
+#undef REAL
+#undef TEST
+#undef RANST_PIECEWISE_LINEAR
+#undef RANST_PIECEWISE_LINEAR_CREATE
+#undef RANST_PIECEWISE_LINEAR_SETUP
+#undef RANST_PIECEWISE_LINEAR_GET
+#undef RANST_PIECEWISE_LINEAR_PDF
+#undef RANST_PIECEWISE_LINEAR_REF_GET
+#undef RANST_PIECEWISE_LINEAR_REF_PUT
+#undef EQ_EPS_R
+#undef EPS_R
+#undef SQRT
+#undef TYPE_FLOAT
diff --git a/src/test_ssp_ran_sphere.c b/src/test_ssp_ran_sphere.c
@@ -26,40 +26,17 @@
* 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 TYPE_FLOAT 1
+#include "test_ssp_ran_sphere.h"
-#define NSAMPS 128
+#define TYPE_FLOAT 0
+#include "test_ssp_ran_sphere.h"
int
main(int argc, char** argv)
{
- struct ssp_rng* rng;
- struct mem_allocator allocator;
- double samps[NSAMPS][4];
- double* 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(d3_is_normalized(f), 1);
- CHECK(eq_eps(samps[i][3], 1.0/(4.0*PI), 1.e-6), 1);
- CHECK(eq_eps(samps[i][3], ssp_ran_sphere_uniform_pdf(), 1.e-6), 1);
- FOR_EACH(j, 0, i) {
- CHECK(d3_eq_eps(samps[j], samps[i], 1.e-6), 0);
- }
- }
-
- ssp_rng_ref_put(rng);
- check_memory_allocator(&allocator);
- mem_shutdown_proxy_allocator(&allocator);
-
- CHECK(mem_allocated_size(), 0);
+ test_float();
+ test_double();
return 0;
}
diff --git a/src/test_ssp_ran_sphere.h b/src/test_ssp_ran_sphere.h
@@ -0,0 +1,99 @@
+/* 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_SPHERE_H
+#define TEST_SSP_RAN_SPHERE_H
+
+#include "ssp.h"
+#include "test_ssp_utils.h"
+
+#define NSAMPS 128
+
+#endif /* TEST_SSP_RAN_SPHERE_H */
+
+#if TYPE_FLOAT==0
+#define REAL double
+#define TEST test_double
+#define RAN_SPHERE_UNIFORM ssp_ran_sphere_uniform
+#define RAN_SPHERE_UNIFORM_PDF ssp_ran_sphere_uniform_pdf
+#define EQ_EPS eq_eps
+#define R3_EQ_EPS d3_eq_eps
+#define R3_IS_NORMALIZED d3_is_normalized
+
+#elif TYPE_FLOAT==1
+#define REAL float
+#define TEST test_float
+#define RAN_SPHERE_UNIFORM ssp_ran_sphere_uniform_float
+#define RAN_SPHERE_UNIFORM_PDF ssp_ran_sphere_uniform_float_pdf
+#define EQ_EPS eq_epsf
+#define R3_EQ_EPS f3_eq_eps
+#define R3_IS_NORMALIZED f3_is_normalized
+
+#else
+#error "TYPE_SUFFIX must be defined either 0 or 1"
+#endif
+
+static void
+TEST()
+{
+ struct ssp_rng* rng;
+ struct mem_allocator allocator;
+ REAL samps[NSAMPS][4];
+ REAL* f = NULL;
+ int i = 0, j = 0;
+
+ 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 = RAN_SPHERE_UNIFORM(rng, samps[i], &samps[i][3]);
+ CHECK(f, samps[i]);
+ CHECK(R3_IS_NORMALIZED(f), 1);
+ CHECK(EQ_EPS(samps[i][3], 1/(4*(REAL)PI), (REAL)1.e-6), 1);
+ CHECK(EQ_EPS(samps[i][3], RAN_SPHERE_UNIFORM_PDF(), (REAL)1.e-6), 1);
+ FOR_EACH(j, 0, i) {
+ CHECK(R3_EQ_EPS(samps[j], samps[i], (REAL)1.e-6), 0);
+ }
+ }
+
+ ssp_rng_ref_put(rng);
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+
+ CHECK(mem_allocated_size(), 0);
+}
+
+#undef REAL
+#undef TEST
+#undef RAN_SPHERE_UNIFORM
+#undef RAN_SPHERE_UNIFORM_PDF
+#undef EQ_EPS
+#undef R3_EQ_EPS
+#undef R3_IS_NORMALIZED
+#undef TYPE_FLOAT
diff --git a/src/test_ssp_ran_triangle.c b/src/test_ssp_ran_triangle.c
@@ -26,97 +26,17 @@
* 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 TYPE_FLOAT 1
+#include "test_ssp_ran_triangle.h"
-#include <rsys/float4.h>
-
-#define NSAMPS 128
+#define TYPE_FLOAT 0
+#include "test_ssp_ran_triangle.h"
int
main(int argc, char** argv)
{
- struct ssp_rng* rng;
- struct mem_allocator allocator;
- double samps[NSAMPS][3];
- double A[3], B[3], C[3];
- double v0[3], v1[3], v2[3], v3[3], v4[3], v5[3];
- double plane[4];
- size_t counter[2];
- size_t nsteps;
- size_t i;
(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, 3) {
- A[i] = ssp_rng_uniform_double(rng, 0.0, 1.0);
- B[i] = ssp_rng_uniform_double(rng, 0.0, 1.0);
- C[i] = ssp_rng_uniform_double(rng, 0.0, 1.0);
- }
-
- d3_sub(v0, B, A);
- d3_sub(v1, C, A);
- d3_sub(v2, C, B);
- d3_minus(v3, v0);
- d3_minus(v4, v1);
- d3_minus(v5, v2);
- d3_cross(plane, v0, v1);
- plane[3] = -d3_dot(plane, C);
-
- FOR_EACH(i, 0, NSAMPS) {
- double tmp0[3], tmp1[3];
- double dot = 0.0;
- double area = 0.0;
-
- CHECK(ssp_ran_triangle_uniform(rng, A, B, C, samps[i]), samps[i]);
- CHECK(eq_eps(d3_dot(plane, samps[i]), -plane[3], 1.e-6), 1);
-
- d3_sub(tmp0, samps[i], A);
- dot = d3_dot(d3_cross(tmp0, tmp0, v0), d3_cross(tmp1, v1, v0));
- CHECK(sign(dot), 1.0);
- d3_sub(tmp0, samps[i], B);
- dot = d3_dot(d3_cross(tmp0, tmp0, v2), d3_cross(tmp1, v3, v2));
- CHECK(sign(dot), 1.0);
- d3_sub(tmp0, samps[i], C);
- dot = d3_dot(d3_cross(tmp0, tmp0, v4), d3_cross(tmp1, v5, v4));
- CHECK(sign(dot), 1.0);
-
- area = d3_len(tmp1) * 0.5f;
- CHECK(eq_eps(1.0 / area, ssp_ran_triangle_uniform_pdf(A, B, C), 1.e-6), 1);
- }
-
- nsteps = 10000;
- counter[0] = counter[1] = 0;
- d3(A, -1.0, 0.0, 0.0);
- d3(B, 1.0, 0.0, 0.0);
- d3(C, 0.0, 1.0, 0.0);
- FOR_EACH(i, 0, nsteps) {
- ssp_ran_triangle_uniform(rng, A, B, C, samps[0]);
- if(samps[0][0] < 0.0)
- counter[0] += 1;
- else
- counter[1] += 1;
- }
- CHECK(labs((long)(counter[1] - counter[0])) < 100, 1);
-
- counter[0] = counter[1] = 0;
- FOR_EACH(i, 0, nsteps) {
- ssp_ran_triangle_uniform(rng, A, B, C, samps[0]);
- if(samps[0][1] < 1 - 1/sqrt(2))
- counter[0] += 1;
- else
- counter[1] += 1;
- }
- CHECK(labs((long)(counter[1] - counter[0])) < 100, 1);
-
- ssp_rng_ref_put(rng);
- check_memory_allocator(&allocator);
- mem_shutdown_proxy_allocator(&allocator);
- CHECK(mem_allocated_size(), 0);
-
+ test_float();
+ test_double();
return 0;
}
-
diff --git a/src/test_ssp_ran_triangle.h b/src/test_ssp_ran_triangle.h
@@ -0,0 +1,170 @@
+/* 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_TRIANGLE_H
+#define TEST_SSP_RAN_TRIANGLE_H
+
+#include "ssp.h"
+#include "test_ssp_utils.h"
+
+#include <rsys/float4.h>
+
+#define NSAMPS 128
+
+#endif /* TEST_SSP_RAN_TRIANGLE_H */
+
+#if TYPE_FLOAT==0
+#define REAL double
+#define TEST test_double
+#define RNG_UNIFORM_R ssp_rng_uniform_double
+#define RAN_TRIANGLE_UNIFORM ssp_ran_triangle_uniform
+#define RAN_TRIANGLE_UNIFORM_PDF ssp_ran_triangle_uniform_pdf
+#define EQ_EPS_R eq_eps
+#define R3 d3
+#define R3_DOT d3_dot
+#define R3_SUB d3_sub
+#define R3_MINUS d3_minus
+#define R3_CROSS d3_cross
+#define R3_LEN d3_len
+
+#elif TYPE_FLOAT==1
+#define REAL float
+#define TEST test_float
+#define RNG_UNIFORM_R ssp_rng_uniform_float
+#define RAN_TRIANGLE_UNIFORM ssp_ran_triangle_uniform_float
+#define RAN_TRIANGLE_UNIFORM_PDF ssp_ran_triangle_uniform_float_pdf
+#define EQ_EPS_R eq_epsf
+#define R3 f3
+#define R3_DOT f3_dot
+#define R3_SUB f3_sub
+#define R3_MINUS f3_minus
+#define R3_CROSS f3_cross
+#define R3_LEN f3_len
+#else
+#error "TYPE_SUFFIX must be defined either 0 or 1"
+#endif
+
+static void
+TEST()
+{
+ struct ssp_rng* rng;
+ struct mem_allocator allocator;
+ REAL samps[NSAMPS][3];
+ REAL A[3], B[3], C[3];
+ REAL v0[3], v1[3], v2[3], v3[3], v4[3], v5[3];
+ REAL plane[4];
+ REAL pdf;
+ size_t counter[2];
+ size_t nsteps;
+ size_t i;
+
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+
+ CHECK(ssp_rng_create(&allocator, &ssp_rng_mt19937_64, &rng), RES_OK);
+
+ FOR_EACH(i, 0, 3) {
+ A[i] = RNG_UNIFORM_R(rng, 0, 1);
+ B[i] = RNG_UNIFORM_R(rng, 0, 1);
+ C[i] = RNG_UNIFORM_R(rng, 0, 1);
+ }
+
+ R3_SUB(v0, B, A);
+ R3_SUB(v1, C, A);
+ R3_SUB(v2, C, B);
+ R3_MINUS(v3, v0);
+ R3_MINUS(v4, v1);
+ R3_MINUS(v5, v2);
+ R3_CROSS(plane, v0, v1);
+ plane[3] = -R3_DOT(plane, C);
+
+ FOR_EACH(i, 0, NSAMPS) {
+ REAL tmp0[3], tmp1[3];
+ REAL dot = 0;
+ REAL area = 0;
+
+ CHECK(RAN_TRIANGLE_UNIFORM(rng, A, B, C, samps[i], &pdf), samps[i]);
+ CHECK(EQ_EPS_R(R3_DOT(plane, samps[i]), -plane[3], (REAL)1.e-6), 1);
+
+ R3_SUB(tmp0, samps[i], A);
+ dot = R3_DOT(R3_CROSS(tmp0, tmp0, v0), R3_CROSS(tmp1, v1, v0));
+ CHECK(sign(dot), 1);
+ R3_SUB(tmp0, samps[i], B);
+ dot = R3_DOT(R3_CROSS(tmp0, tmp0, v2), R3_CROSS(tmp1, v3, v2));
+ CHECK(sign(dot), 1);
+ R3_SUB(tmp0, samps[i], C);
+ dot = R3_DOT(R3_CROSS(tmp0, tmp0, v4), R3_CROSS(tmp1, v5, v4));
+ CHECK(sign(dot), 1);
+
+ area = R3_LEN(tmp1) * 0.5f;
+ CHECK(EQ_EPS_R(1 / area, RAN_TRIANGLE_UNIFORM_PDF(A, B, C), (REAL)1.e-6), 1);
+ }
+
+ nsteps = 10000;
+ counter[0] = counter[1] = 0;
+ R3(A, -1, 0, 0);
+ R3(B, 1, 0, 0);
+ R3(C, 0, 1, 0);
+ FOR_EACH(i, 0, nsteps) {
+ RAN_TRIANGLE_UNIFORM(rng, A, B, C, samps[0], NULL);
+ if(samps[0][0] < 0.0)
+ counter[0] += 1;
+ else
+ counter[1] += 1;
+ }
+ CHECK(labs((long)(counter[1] - counter[0])) < 100, 1);
+
+ counter[0] = counter[1] = 0;
+ FOR_EACH(i, 0, nsteps) {
+ RAN_TRIANGLE_UNIFORM(rng, A, B, C, samps[0], NULL);
+ if(samps[0][1] < 1 - 1/sqrt(2))
+ counter[0] += 1;
+ else
+ counter[1] += 1;
+ }
+ CHECK(labs((long)(counter[1] - counter[0])) < 100, 1);
+
+ ssp_rng_ref_put(rng);
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHECK(mem_allocated_size(), 0);
+}
+
+
+#undef REAL
+#undef TEST
+#undef RNG_UNIFORM_R
+#undef RAN_TRIANGLE_UNIFORM
+#undef RAN_TRIANGLE_UNIFORM_PDF
+#undef EQ_EPS_R
+#undef R3
+#undef R3_DOT
+#undef R3_SUB
+#undef R3_MINUS
+#undef R3_CROSS
+#undef R3_LEN
+#undef TYPE_FLOAT
diff --git a/src/test_ssp_ran_uniform_disk.c b/src/test_ssp_ran_uniform_disk.c
@@ -28,67 +28,17 @@
* 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"
-#include <rsys/math.h>
+#define TYPE_FLOAT 1
+#include "test_ssp_ran_uniform_disk.h"
-#include <string.h>
-
-#define NBS 1000000
+#define TYPE_FLOAT 0
+#include "test_ssp_ran_uniform_disk.h"
int
main(int argc, char** argv)
{
- struct ssp_rng_proxy* proxy;
- struct ssp_rng* rng;
- struct mem_allocator allocator;
- int i;
- unsigned r, c, nb = 0;
- double pt[2];
- double x = 0, x2 = 0;
- double mean, std;
- double exp_mean = 20 * 20 * NBS / (PI * 100 * 100 ), exp_std = 0;
- unsigned counts[10][10];
(void)argc, (void)argv;
-
- mem_init_proxy_allocator(&allocator, &mem_default_allocator);
-
- CHECK(ssp_rng_proxy_create(&allocator, &ssp_rng_threefry, 1, &proxy), RES_OK);
- CHECK(ssp_rng_proxy_create_rng(proxy, 0, &rng), RES_OK);
-
- memset(counts, 0, sizeof(counts));
- for (i = 0; i < NBS; i++) {
- ssp_ran_uniform_disk(rng, 100, pt);
- /* locate pt in a 10x10 grid */
- r = (unsigned)((100 + pt[0]) / 20);
- c = (unsigned)((100 + pt[1]) / 20);
- ++counts[r][c];
- }
-
- for (r = 0; r < 10; r++) {
- unsigned _x = (r >= 5 ? r - 4 : r - 5) * 20;
- for (c = 0; c < 10; c++) {
- unsigned _y = (c >= 5 ? c - 4 : c - 5) * 20;
- unsigned r2 = _x * _x + _y * _y;
- if (r2 > 100 * 100)
- /* this square is not (fully) in the disk */
- continue;
- ++nb;
- x += counts[r][c];
- x2 += counts[r][c] * counts[r][c];
- }
- }
- mean = x/nb;
- std = sqrt(x2/nb - x/nb*x/nb);
- printf("%g %g\n", mean, std);
- CHECK(fabs(mean - exp_mean) < 10, 1);
- CHECK(fabs(std - exp_std) < 100, 1);
-
- CHECK(ssp_rng_ref_put(rng), RES_OK);
- CHECK(ssp_rng_proxy_ref_put(proxy), RES_OK);
-
- check_memory_allocator(&allocator);
- mem_shutdown_proxy_allocator(&allocator);
- CHECK(mem_allocated_size(), 0);
+ test_float();
+ test_double();
return 0;
}
diff --git a/src/test_ssp_ran_uniform_disk.h b/src/test_ssp_ran_uniform_disk.h
@@ -0,0 +1,152 @@
+/* Copyright (C) |Meso|Star> 2015-2017 (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
+ * 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_UNIFORM_DISK_H
+#define TEST_SSP_RAN_UNIFORM_DISK_H
+
+#include "ssp.h"
+#include "test_ssp_utils.h"
+#include <rsys/math.h>
+
+#include <string.h>
+
+#define NBS 1000000
+
+#endif /* TEST_SSP_RAN_UNIFORM_DISK_H */
+
+#if TYPE_FLOAT==0
+#define REAL double
+#define TEST test_double
+#define RNG_UNIFORM_R ssp_rng_uniform_double
+#define RNG_UNIFORM_DISK_LOCAL ssp_ran_uniform_disk_local
+#define RNG_UNIFORM_DISK ssp_ran_uniform_disk
+#define R33_BASIS d33_basis
+#define R33_MULR3 d33_muld3
+#define R3_EQ_EPS d3_eq_eps
+#define R3_NORMALIZE d3_normalize
+#define SQRT sqrt
+
+#elif TYPE_FLOAT==1
+#define REAL float
+#define TEST test_float
+#define RNG_UNIFORM_R ssp_rng_uniform_float
+#define RNG_UNIFORM_DISK_LOCAL ssp_ran_uniform_disk_float_local
+#define RNG_UNIFORM_DISK ssp_ran_uniform_disk_float
+#define R33_BASIS f33_basis
+#define R33_MULR3 f33_mulf3
+#define R3_EQ_EPS f3_eq_eps
+#define R3_NORMALIZE f3_normalize
+#define SQRT sqrtf
+#else
+#error "TYPE_SUFFIX must be defined either 0 or 1"
+#endif
+
+static void
+TEST()
+{
+ struct ssp_rng* rng0, *rng1;
+ struct mem_allocator allocator;
+ int i;
+ int r, c, nb = 0;
+ REAL pt[3], pt2[3], up[3];
+ REAL frame[9];
+ REAL x_sum = 0, x2_sum = 0;
+ REAL mean, std;
+ REAL exp_mean = NBS * (20 * 20 / ((REAL)PI * 100 * 100)), exp_std = 0;
+ unsigned counts[10][10];
+ uint64_t seed = 1234;
+
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+
+ CHECK(ssp_rng_create(&allocator, &ssp_rng_threefry, &rng0), RES_OK);
+ CHECK(ssp_rng_create(&allocator, &ssp_rng_threefry, &rng1), RES_OK);
+
+ up[0] = RNG_UNIFORM_R(rng1, -1, 1);
+ up[1] = RNG_UNIFORM_R(rng1, -1, 1);
+ up[2] = RNG_UNIFORM_R(rng1, -1, 1);
+ R3_NORMALIZE(up, up);
+ R33_BASIS(frame, up);
+
+ ssp_rng_set(rng0, seed);
+ ssp_rng_set(rng1, seed);
+
+ memset(counts, 0, sizeof(counts));
+ for(i = 0; i < NBS; i++) {
+ REAL tmp[3];
+ RNG_UNIFORM_DISK_LOCAL(rng0, 100, pt, NULL);
+ RNG_UNIFORM_DISK(rng1, 100, up, pt2, NULL);
+ R33_MULR3(tmp, frame, pt);
+ CHECK(R3_EQ_EPS(tmp, pt2, (REAL)1.e-6), 1);
+ ASSERT(pt[2] == 0);
+ /* locate pt in a 10x10 grid */
+ r = (int)((100 + pt[0]) / 20);
+ c = (int)((100 + pt[1]) / 20);
+ ++counts[r][c];
+ }
+
+ for(r = 0; r < 10; r++) {
+ int x = (r >= 5 ? r - 4 : r - 5) * 20;
+ for(c = 0; c < 10; c++) {
+ int y = (c >= 5 ? c - 4 : c - 5) * 20;
+ int r2 = x * x + y * y;
+ if(r2 > 100 * 100)
+ /* this square is not (fully) in the disk */
+ continue;
+ ++nb;
+ x_sum += (REAL)counts[r][c];
+ x2_sum += (REAL)(counts[r][c] * counts[r][c]);
+ }
+ }
+ mean = x_sum / (REAL)nb;
+ std = x2_sum / (REAL)nb - mean*mean;
+ std = std > 0 ? SQRT(std) : 0;
+ printf("%g %g\n", mean, std);
+ CHECK(fabs(mean - exp_mean) < 10, 1);
+ CHECK(fabs(std - exp_std) < 200, 1);
+
+ CHECK(ssp_rng_ref_put(rng0), RES_OK);
+ CHECK(ssp_rng_ref_put(rng1), RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHECK(mem_allocated_size(), 0);
+}
+
+#undef REAL
+#undef TEST
+#undef RNG_UNIFORM_R
+#undef RNG_UNIFORM_DISK_LOCAL
+#undef RNG_UNIFORM_DISK
+#undef R33_BASIS
+#undef R33_MULR3
+#undef R3_EQ_EPS
+#undef R3_NORMALIZE
+#undef SQRT
+#undef TYPE_FLOAT
diff --git a/src/test_ssp_rng.h b/src/test_ssp_rng.h
@@ -0,0 +1,32 @@
+/* Copyright (C) |Meso|Star> 2015-2017 (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
+* 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. */
+
+
+/* This file is intentionally empty */
diff --git a/src/test_ssp_rng_proxy.h b/src/test_ssp_rng_proxy.h
@@ -0,0 +1,32 @@
+/* Copyright (C) |Meso|Star> 2015-2017 (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
+* 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. */
+
+
+/* This file is intentionally empty */
