star-sp

commit 0f117c384b39500bcf13c67551226e7d963d98e8
parent abd4afca65ec25ca19208670afd1f8723fcfa7db
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri, 21 Oct 2016 12:29:49 +0200

Speed up the canonical uniform distributions

Pre-compute RNG constants at the construction of the RNG. Avoid the
evaluation of a log2 at each invocation of the canonical distribution.

Diffstat:
M
src/ssp_rng.c
 | 
36
+++++++++++++++++++++++++++++-------
M
src/ssp_rng_c.h
 | 
6
++++++

2 files changed, 35 insertions(+), 7 deletions(-)
diff --git a/src/ssp_rng.c b/src/ssp_rng.c
@@ -69,7 +69,6 @@ struct rng_kiss { uint32_t x, y, z, c; };
 
 static res_T
 rng_kiss_set(void* data, const size_t seed)
-
 {
   struct rng_kiss* kiss = (struct rng_kiss*)data;
   RAN_NAMESPACE::mt19937 rng_mt(seed % UINT32_MAX);
@@ -462,6 +461,7 @@ ssp_rng_create
 {
   struct mem_allocator* allocator;
   struct ssp_rng* rng = NULL;
+  size_t log2r;
   res_T res = RES_OK;
 
   if(!type || !out_rng || !check_type(type)) {
@@ -484,12 +484,15 @@ ssp_rng_create
     goto error;
   }
   res = type->init(allocator, rng->state);
-  if (res != RES_OK)
-  {
-    goto error;
-  }
+  if(res != RES_OK) goto error;
   rng->type = *type;
 
+  /* Precompute some values for the canonical distribution */
+  rng->r = (long double)rng->type.max - (long double)rng->type.min + 1.0L;
+  log2r = (size_t)(std::log(rng->r) / std::log(2.0L));
+  rng->dbl_k = std::max<size_t>(1UL, (53 + log2r - 1UL) / log2r);
+  rng->flt_k = std::max<size_t>(1UL, (24 + log2r - 1UL) / log2r);
+
 exit:
   if(out_rng) *out_rng = rng;
   return res;
@@ -563,9 +566,21 @@ double
 ssp_rng_canonical(struct ssp_rng* rng)
 {
   if(!rng) FATAL("The Random Number Generator is NULL\n");
+#if 0
   rng_cxx rng_cxx(*rng);
   return RAN_NAMESPACE::generate_canonical
     <double, std::numeric_limits<double>::digits>(rng_cxx);
+#else
+  /* Optimize version */
+  size_t k = rng->dbl_k;
+  double sum = 0;
+  double tmp = 1;
+  for(; k !=0; --k) {
+    sum += (double)(rng->type.get(rng->state) - rng->type.min) * tmp;
+    tmp = (double)(tmp*rng->r);
+  }
+  return sum/tmp;
+#endif
 }
 
 float
@@ -582,8 +597,15 @@ ssp_rng_canonical_float(struct ssp_rng* rng)
 #else
   float r;
   do {
-    r = RAN_NAMESPACE::generate_canonical
-    <float, std::numeric_limits<float>::digits>(rng_cxx);
+    /* Optimise version of the generate_canonical function */
+    size_t k = rng->flt_k;
+    float sum = 0;
+    float tmp = 1;
+    for(; k !=0; --k) {
+      sum += (float)(rng->type.get(rng->state) - rng->type.min) * tmp;
+      tmp = (float)(tmp*rng->r);
+    }
+    r = sum/tmp;
   } while(r >= 1);
   return r;
 #endif
diff --git a/src/ssp_rng_c.h b/src/ssp_rng_c.h
@@ -60,6 +60,12 @@ struct ssp_rng {
   struct ssp_rng_type type;
   void* state;
   struct mem_allocator* allocator;
+
+  /* Precomputed RNG constants used to speed up the canonical generations */
+  long double r; /* max - min + 1 */
+  size_t dbl_k; /* max(1, (#bits_mantisse_double + log2(r)) / log2(r)) */
+  size_t flt_k; /* max(1, (#bits_mantisse_float  + log2(r)) / log2(r)) */
+
   ref_T ref;
 };