Random: Implemented fastrand_double().

diff --git a/ucw/random-fast.c b/ucw/random-fast.c
index 672ccd883f4a260aaca26d31593c8fc1ff4c0ef0..9306b484d16e60f5e89239dc5e3b89a6afb2bd07 100644 (file)
--- a/ucw/random-fast.c
+++ b/ucw/random-fast.c
@@ -14,7 +14,9 @@
 #include <ucw/threads.h>
 #include <ucw/unaligned.h>
 
+#include <float.h>
 #include <limits.h>
+#include <math.h>
 #include <stdlib.h>
 #include <unistd.h>
 #include <sys/time.h>
@@ -607,3 +609,38 @@ void fastrand_mem(struct fastrand *c, void *ptr, size_t size)
        }
     }
 }
+
+double fastrand_double(struct fastrand *c)
+{
+  // Generate 52 random bits
+  FASTRAND_CHECK_SEED(c);
+  u64 r = fastrand_one(c);
+#if FASTRAND_ONE_BITS >= 52
+  r >>= FASTRAND_ONE_BITS - 52;
+#else
+  // Do we want to generate doubles with full 52 precision?
+#if 0
+  r = (r << (52 - FASTRAND_ONE_BITS)) | (fastrand_one(c) >> (2 * FASTRAND_ONE_BITS - 52));
+#else
+  r <<= 52 - FASTRAND_ONE_BITS;
+#endif
+#endif
+
+  // With standard IEC 60559 / IEEE 754 format, we can directly store @r as a double
+  // from range [1.0, 2.0) and subtract 1.0.
+  // But the alternative with multiplication is often faster anyway.
+#if 0 && defined(__STDC_IEC_559__)
+  union {
+    u64 u;
+    double d;
+  } tmp;
+  ASSERT(sizeof(double) == 8 && sizeof(tmp) == 8);
+  tmp.u = (0x3ffLLU << 52) | r;
+  return tmp.d - 1.0;
+#elif 1
+  static const double mul = 1.0 / (1LLU << 52);
+  return r * mul;
+#else
+  return ldexp(r, -52);
+#endif
+}

diff --git a/ucw/random-test.c b/ucw/random-test.c
index 5f76423f6c66f9eebc0ff2d6c9457809e27fd7b1..1e51aa74fb97ce7599b5e9cc915398ece02446cb 100644 (file)
--- a/ucw/random-test.c
+++ b/ucw/random-test.c
@@ -13,6 +13,7 @@
 #include <ucw/time.h>
 
 #include <errno.h>
+#include <float.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <time.h>
@@ -69,6 +70,10 @@ static void test_fast(void)
       TEST((y = fastrand_max(r, x)) < x);
       x = fastrand_u64(r) ? : 1;
       TEST((y = fastrand_max_u64(r, x)) < x);
+
+      double d;
+      d = fastrand_double(r);
+      TEST(d >= 0.0 - DBL_EPSILON && d < 1.0 + DBL_EPSILON);
     }
   fastrand_delete(r);
   printf("OK\n");
@@ -173,6 +178,7 @@ static void bench_fast(void)
   BENCH_SLOW(4, fastrand_max_u64(r, (1LLU << 63) + 1));
   BENCH_SLOW(2, fastrand_max(r, 1 + fastrand_max(r, 1000000)));
   BENCH_SLOW(4, fastrand_max_u64(r, 1 + fastrand_max_u64(r, 1000000000000LL)));
+  BENCH_SLOW(4, fastrand_double(r));
 
   byte buf[1000];
   BENCH(fastrand_mem(r, buf, 1));

diff --git a/ucw/random.h b/ucw/random.h
index 703b1b09a7660bc5002359122123a7ac9d048ca1..0a287d189f7885486f17b6b5f27c869b13a54ee7 100644 (file)
--- a/ucw/random.h
+++ b/ucw/random.h
@@ -24,6 +24,7 @@
 #define fastrand_max ucw_fastrand_max
 #define fastrand_max_u64 ucw_fastrand_max_u64
 #define fastrand_mem ucw_fastrand_mem
+#define fastrand_double ucw_fastrand_double
 #endif
 
 /* random-fast.c */
@@ -94,4 +95,7 @@ u64 fastrand_max_u64(struct fastrand *c, u64 bound);
 /** Generate @size random bytes. **/
 void fastrand_mem(struct fastrand *c, void *ptr, size_t size);
 
+/** Uniformly generate a random number from range [0.0, 1.0) + possible inaccuracies from basic floating point operations. **/
+double fastrand_double(struct fastrand *c);
+
 #endif