From ef604b4dc939382985659b3c91599c28b7e9946c Mon Sep 17 00:00:00 2001
From: Pavel Charvat <pchar@ucw.cz>
Date: Tue, 22 Mar 2022 15:46:30 +0000
Subject: [PATCH] Random: Implemented fastrand_double().

---
 ucw/random-fast.c | 37 +++++++++++++++++++++++++++++++++++++
 ucw/random-test.c |  6 ++++++
 ucw/random.h      |  4 ++++
 3 files changed, 47 insertions(+)

diff --git a/ucw/random-fast.c b/ucw/random-fast.c
index 672ccd88..9306b484 100644
--- 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 5f76423f..1e51aa74 100644
--- 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 703b1b09..0a287d18 100644
--- 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
-- 
2.39.2