star-sp

commit 7b4b0c0d8ec8df467dfaee37e5410fb0ec31afe1
parent cc3fc4074314b7f03b350f3f3e744ed9f93b048b
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Thu,  9 Apr 2015 21:01:01 +0200

Add and test the ranlux48 PRNG

Diffstat:
M
cmake/CMakeLists.txt
 | 
1
+
M
src/ssp.h
 | 
2
++
M
src/ssp_rng.c
 | 
92
+++++++++++++++++++++++++++++++++++++++++++++++++------------------------------
M
src/test_ssp_rng.c
 | 
21
++++++++++++++++++++-

4 files changed, 80 insertions(+), 36 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -101,6 +101,7 @@ if(NOT NO_TEST)
   build_test(test_ssp_rng)
   add_test(test_ssp_rng_kiss test_ssp_rng kiss)
   add_test(test_ssp_rng_mt19937_64 test_ssp_rng mt19937_64)
+  add_test(test_ssp_rng_mt19937_64 test_ssp_rng ranlux48)
 
   new_test(test_ssp_ran_hemisphere m)
 endif()
diff --git a/src/ssp.h b/src/ssp.h
@@ -84,6 +84,8 @@ BEGIN_DECLS
 SSP_API struct ssp_rng_type ssp_rng_kiss;
 /* 64-bits Mersenne Twister builtin PRNG type of Matsumoto and Nishimura, 2000 */
 SSP_API struct ssp_rng_type ssp_rng_mt19937_64;
+/* 48-bits RANLUX builtin PRNG type of Lusher and James, 1994 */
+SSP_API struct ssp_rng_type ssp_rng_ranlux48;
 
 /*******************************************************************************
  * API declaration
diff --git a/src/ssp_rng.c b/src/ssp_rng.c
@@ -131,80 +131,102 @@ struct ssp_rng_type ssp_rng_kiss = {
 };
 
 /*******************************************************************************
- * 64-bits Mersenne Twister PRNG
+ * C++11 PRNG
  ******************************************************************************/
+template<typename RNG>
 static void
-rng_mt19937_64_set(void* data, const unsigned long seed)
+rng_cxx_set(void* data, const unsigned long seed)
 {
-  std::mt19937_64* mt = (std::mt19937_64*)data;
-  ASSERT(mt);
-  mt->seed(seed);
+  RNG* rng = (RNG*)data;
+  ASSERT(rng);
+  rng->seed(seed);
 }
 
+template<typename RNG>
 static unsigned long
-rng_mt19937_64_get(void* data)
+rng_cxx_get(void* data)
 {
-  std::mt19937_64* mt = (std::mt19937_64*)data;
-  ASSERT(mt);
-  return (*mt)();
+  RNG* rng = (RNG*)data;
+  ASSERT(rng);
+  return (*rng)();
 }
 
+template<typename RNG>
 static unsigned long
-rng_mt19937_64_uniform_int
+rng_cxx_uniform_int
   (void* data, const unsigned long lower, const unsigned long upper)
 {
-  std::mt19937_64* mt = (std::mt19937_64*)data;
-  ASSERT(mt && lower <= upper);
-  return std::uniform_int_distribution<unsigned long>(lower, upper)(*mt);
+  RNG* rng = (RNG*)data;
+  ASSERT(rng && lower <= upper);
+  return std::uniform_int_distribution<unsigned long>(lower, upper)(*rng);
 }
 
+template<typename RNG>
 static double
-rng_mt19937_64_uniform_double
+rng_cxx_uniform_double
   (void* data, const double lower, const double upper)
 {
-  std::mt19937_64* mt = (std::mt19937_64*)data;
-  ASSERT(mt && lower <= upper);
-  return std::uniform_real_distribution<>(lower, upper)(*mt);
+  RNG* rng = (RNG*)data;
+  ASSERT(rng && lower <= upper);
+  return std::uniform_real_distribution<>(lower, upper)(*rng);
 }
 
+template<typename RNG>
 static double
-rng_mt19937_64_canonical(void* data)
+rng_cxx_canonical(void* data)
 {
-  std::mt19937_64* mt = (std::mt19937_64*)data;
-  ASSERT(mt);
-  return std::generate_canonical<double, 48>(*mt);
+  RNG* rng = (RNG*)data;
+  ASSERT(rng);
+  return std::generate_canonical<double, 48>(*rng);
 }
 
+template<typename RNG>
 static void
-rng_mt19937_64_init(struct mem_allocator* allocator, void* data)
+rng_cxx_init(struct mem_allocator* allocator, void* data)
 {
   (void)allocator;
   ASSERT(data);
-  new (data) std::mt19937_64;
+  new (data) RNG;
 }
 
+
+template<typename RNG>
 static void
-rng_mt19937_64_release(void* data)
+rng_cxx_release(void* data)
 {
-  std::mt19937_64* mt = (std::mt19937_64*)data;
-  ASSERT(mt);
-  mt->~mersenne_twister_engine();
+  RNG* rng = (RNG*)data;
+  ASSERT(rng);
+  rng->~RNG();
 }
 
-/* Exported type */
+/* 64-bits Mersenne Twister PRNG */
 struct ssp_rng_type ssp_rng_mt19937_64 = {
-  rng_mt19937_64_init,
-  rng_mt19937_64_release,
-  rng_mt19937_64_set,
-  rng_mt19937_64_get,
-  rng_mt19937_64_uniform_int,
-  rng_mt19937_64_uniform_double,
-  rng_mt19937_64_canonical,
+  rng_cxx_init<std::mt19937_64>,
+  rng_cxx_release<std::mt19937_64>,
+  rng_cxx_set<std::mt19937_64>,
+  rng_cxx_get<std::mt19937_64>,
+  rng_cxx_uniform_int<std::mt19937_64>,
+  rng_cxx_uniform_double<std::mt19937_64>,
+  rng_cxx_canonical<std::mt19937_64>,
   std::mt19937_64::min(),
   std::mt19937_64::max(),
   sizeof(std::mt19937_64)
 };
 
+/* 48-bits RANLUX PRNG */
+struct ssp_rng_type ssp_rng_ranlux48 = {
+  rng_cxx_init<std::ranlux48>,
+  rng_cxx_release<std::ranlux48>,
+  rng_cxx_set<std::ranlux48>,
+  rng_cxx_get<std::ranlux48>,
+  rng_cxx_uniform_int<std::ranlux48>,
+  rng_cxx_uniform_double<std::ranlux48>,
+  rng_cxx_canonical<std::ranlux48>,
+  std::ranlux48::min(),
+  std::ranlux48::max(),
+  sizeof(std::ranlux48)
+};
+
 /*******************************************************************************
  * Helper functions
  ******************************************************************************/
diff --git a/src/test_ssp_rng.c b/src/test_ssp_rng.c
@@ -31,6 +31,8 @@
 #include "ssp.h"
 #include "test_ssp_utils.h"
 
+#include <rsys/clock_time.h>
+
 #include <string.h>
 
 #define NRAND 1024
@@ -40,9 +42,11 @@ test_rng(struct ssp_rng_type* type)
 {
   struct ssp_rng* rng;
   struct mem_allocator allocator;
+  struct time t0, t1;
   unsigned long datai0[NRAND];
   unsigned long datai1[NRAND];
   double dataf[NRAND];
+  char buf[512];
   int i, j;
 
   mem_init_proxy_allocator(&allocator, &mem_default_allocator);
@@ -115,6 +119,15 @@ test_rng(struct ssp_rng_type* type)
     }
   }
 
+  time_current(&t0);
+  FOR_EACH(i, 0, 1000000) {
+    ssp_rng_get(rng);
+  }
+  time_current(&t1);
+  time_sub(&t0, &t1, &t0);
+  time_dump(&t0, TIME_SEC|TIME_MSEC|TIME_USEC, NULL, buf, sizeof(buf));
+  printf("1,000,000 random numbers in %s\n", buf);
+
   CHECK(ssp_rng_ref_put(rng), RES_OK);
   check_memory_allocator(&allocator);
   mem_shutdown_proxy_allocator(&allocator);
@@ -124,14 +137,20 @@ int
 main(int argc, char** argv)
 {
   if(argc <= 1) {
-    fprintf(stderr, "Usage: %s <kiss|mt19937_64>\n", argv[0]);
+    fprintf(stderr, "Usage: %s <kiss|mt19937_64|ranlux48>\n", argv[0]);
     exit(0);
   }
   if(!strcmp(argv[1], "kiss")) {
     test_rng(&ssp_rng_kiss);
   } else if(!strcmp(argv[1], "mt19937_64")) {
     test_rng(&ssp_rng_mt19937_64);
+  } else if(!strcmp(argv[1], "ranlux48")) {
+    test_rng(&ssp_rng_ranlux48);
+  } else {
+    fprintf(stderr, "Unknown RNG `%s'\n", argv[1]);
+    CHECK(0, 1);
   }
+
   CHECK(mem_allocated_size(), 0);
   return 0;
 }