38 #include <visp3/core/vpConfig.h>
40 #ifdef VISP_HAVE_CATCH2
41 #define CATCH_CONFIG_RUNNER
44 #include <visp3/core/vpGaussRand.h>
45 #include <visp3/core/vpMath.h>
46 #include <visp3/core/vpTime.h>
48 #ifdef ENABLE_VISP_NAMESPACE
67 x = a * (x - k * q) - k * r;
78 const long modulo = ntab - 2;
87 for (j = 0; j < ntab; j++) {
99 double ans = (double)y / normalizer;
113 VP_EXPLICIT vpUniRandOld(
const long seed = 0)
114 : a(16807), m(2147483647), q(127773), r(2836), normalizer(2147484721.0), x((seed) ? seed : 739806647)
118 virtual ~vpUniRandOld() { };
121 double operator()() {
return draw1(); }
125 TEST_CASE(
"Check Gaussian draw",
"[visp_rand]")
127 std::vector<double> vec(100000);
128 const double sigma = 5.0, mean = -7.5;
133 for (
size_t i = 0; i < vec.size(); i++) {
138 std::cout << vec.size() <<
" Gaussian draw performed in " << chrono.
getDurationMs() <<
" ms" << std::endl;
141 std::cout <<
"Calculated sigma: " << calculated_sigma << std::endl;
142 std::cout <<
"Calculated mean: " << calculated_mean << std::endl;
144 CHECK(calculated_sigma == Approx(sigma).epsilon(0.01));
145 CHECK(calculated_mean == Approx(mean).epsilon(0.01));
148 TEST_CASE(
"Check Gaussian draw independance",
"[visp_rand]")
150 const double sigma = 5.0, mean = -7.5;
152 SECTION(
"Two simultaneous vpGaussRand instances with the same seed should produce the same results")
156 for (
int i = 0; i < 10; i++) {
157 CHECK(rng1() == Approx(rng2()).margin(1e-6));
160 SECTION(
"Two vpGaussRand instances with the same seed should produce the same results")
162 std::vector<double> vec1, vec2;
165 for (
int i = 0; i < 10; i++) {
166 vec1.push_back(rng());
171 for (
int i = 0; i < 10; i++) {
172 vec2.push_back(rng());
175 REQUIRE(vec1.size() == vec2.size());
177 for (
size_t i = 0; i < vec1.size(); i++) {
178 CHECK(vec1[i] == Approx(vec2[i]).margin(1e-6));
183 TEST_CASE(
"Check uniform draw",
"[visp_rand]")
185 const int niters = 500000;
194 for (
int i = 0; i < niters; i++) {
198 if (sqrt(x * x + y * y) <= 1.0) {
202 double pi = 4.0 * inside / niters;
205 double pi_error = pi - M_PI;
206 std::cout <<
"vpUniRand calculated pi: " << pi <<
" in " << chrono.
getDurationMs() <<
" ms" << std::endl;
207 std::cout <<
"pi error: " << pi_error << std::endl;
209 CHECK(pi == Approx(M_PI).margin(0.005));
212 SECTION(
"C++ rand()")
219 for (
int i = 0; i < niters; i++) {
220 double x =
static_cast<double>(rand()) / RAND_MAX;
221 double y =
static_cast<double>(rand()) / RAND_MAX;
223 if (sqrt(x * x + y * y) <= 1.0) {
227 double pi = 4.0 * inside / niters;
230 double pi_error = pi - M_PI;
231 std::cout <<
"C++ rand() calculated pi: " << pi <<
" in " << chrono.
getDurationMs() <<
" ms" << std::endl;
232 std::cout <<
"pi error: " << pi_error << std::endl;
234 CHECK(pi == Approx(M_PI).margin(0.01));
237 SECTION(
"Old ViSP vpUniRand implementation")
244 for (
int i = 0; i < niters; i++) {
248 if (sqrt(x * x + y * y) <= 1.0) {
252 double pi = 4.0 * inside / niters;
255 double pi_error = pi - M_PI;
256 std::cout <<
"Old ViSP vpUniRand implementation calculated pi: " << pi <<
" in " << chrono.
getDurationMs() <<
" ms"
258 std::cout <<
"pi error: " << pi_error << std::endl;
260 CHECK(pi == Approx(M_PI).margin(0.005));
264 TEST_CASE(
"Check uniform draw range",
"[visp_rand]")
266 const int niters = 1000;
269 SECTION(
"Check[0.0, 1.0) range")
271 for (
int i = 0; i < niters; i++) {
278 SECTION(
"Check[-7, 10) range")
280 const int a = -7, b = 10;
281 for (
int i = 0; i < niters; i++) {
288 SECTION(
"Check[-4.5f, 105.3f) range")
290 const float a = -4.5f, b = 105.3f;
291 for (
int i = 0; i < niters; i++) {
298 SECTION(
"Check[14.6, 56.78) range")
300 const double a = 14.6, b = 56.78;
301 for (
int i = 0; i < niters; i++) {
309 TEST_CASE(
"Check uniform draw independance",
"[visp_rand]")
311 SECTION(
"Two simultaneous vpUniRand instances with the same seed should produce the same results")
316 for (
int i = 0; i < 10; i++) {
323 for (
int i = 0; i < 10; i++) {
324 CHECK(rng1.
uniform(-1.0, 1.0) == Approx(rng2.
uniform(-1.0, 1.0)).margin(1e-6));
328 SECTION(
"Two vpUniRand instances with the same seed should produce the same results")
331 std::vector<uint32_t> vec1, vec2;
334 for (
int i = 0; i < 10; i++) {
335 vec1.push_back(rng.
next());
340 for (
int i = 0; i < 10; i++) {
341 vec2.push_back(rng.
next());
344 REQUIRE(vec1.size() == vec2.size());
346 for (
size_t i = 0; i < vec1.size(); i++) {
347 CHECK(vec1[i] == vec2[i]);
351 std::vector<double> vec1, vec2;
354 for (
int i = 0; i < 10; i++) {
355 vec1.push_back(rng.
uniform(-1.0, 2.0));
360 for (
int i = 0; i < 10; i++) {
361 vec2.push_back(rng.
uniform(-1.0, 2.0));
364 REQUIRE(vec1.size() == vec2.size());
366 for (
size_t i = 0; i < vec1.size(); i++) {
367 CHECK(vec1[i] == Approx(vec2[i]).margin(1e-6));
373 int main(
int argc,
char *argv[])
375 Catch::Session session;
378 session.applyCommandLine(argc, argv);
380 int numFailed = session.run();
388 int main() {
return EXIT_SUCCESS; }
void start(bool reset=true)
Class for generating random number with normal probability density.
static double getStdev(const std::vector< double > &v, bool useBesselCorrection=false)
static double getMean(const std::vector< double > &v)
Class for generating random numbers with uniform probability density.
int uniform(int a, int b)