~ chicken-core (chicken-5) 6a3a8253525816e876ca48398cbb2b3e8cf3dd11

commit 6a3a8253525816e876ca48398cbb2b3e8cf3dd11
Author:     Kooda <kooda@upyum.com>
AuthorDate: Sat Nov 4 13:10:21 2017 +0100
Commit:     Kooda <kooda@upyum.com>
CommitDate: Sat Nov 4 13:12:20 2017 +0100

    Add pseudo-random-real to the chicken.random module.
    
    Based on Taylor R. Campbell algorithm explained on
    https://mumble.net/~campbell/2014/04/28/uniform-random-float

diff --git a/chicken.h b/chicken.h
index 357928b9..b0a7de30 100644
--- a/chicken.h
+++ b/chicken.h
@@ -2095,6 +2095,7 @@ C_fctexport void *C_lookup_procedure_ptr(C_char *id);
 
 C_fctexport C_word C_random_fixnum(C_word n) C_regparm;
 C_fctexport C_word C_fcall C_s_a_u_i_random_int(C_word **ptr, C_word n, C_word rn) C_regparm;
+C_fctexport C_word C_fcall C_a_i_random_real(C_word **ptr, C_word n) C_regparm;
 C_fctexport C_word C_random_bytes(C_word buf, C_word size);
 C_fctexport C_word C_set_random_seed(C_word buf, C_word n);
 
diff --git a/extras.scm b/extras.scm
index 10e1bf0c..0eaacd32 100644
--- a/extras.scm
+++ b/extras.scm
@@ -644,7 +644,7 @@
 ;;; Random numbers:
 
 (module chicken.random
-  (set-pseudo-random-seed! pseudo-random-integer random-bytes)
+  (set-pseudo-random-seed! pseudo-random-integer pseudo-random-real random-bytes)
 
 (import scheme chicken chicken.time chicken.io foreign)
 
@@ -668,6 +668,9 @@
         (else
           (##core#inline_allocate ("C_s_a_u_i_random_int" 2) n))))
 
+(define (pseudo-random-real)
+  (##core#inline_allocate ("C_a_i_random_real" 2)))
+
 (define random-bytes
   (let ((nstate (foreign-value "C_RANDOM_STATE_SIZE" unsigned-int)))
     (lambda (#!optional buf size)
diff --git a/runtime.c b/runtime.c
index 9d312478..5dc17466 100644
--- a/runtime.c
+++ b/runtime.c
@@ -12676,6 +12676,38 @@ C_s_a_u_i_random_int(C_word **ptr, C_word n, C_word rn)
   return C_bignum_simplify(result);
 }
 
+/*
+ * C_a_i_random_real: Generate a stream of bits uniformly at random and
+ * interpret it as the fractional part of the binary expansion of a
+ * number in [0, 1], 0.00001010011111010100...; then round it.
+ * More information on https://mumble.net/~campbell/2014/04/28/uniform-random-float
+ */
+#define random64() ((((C_u64) random_word()) << 32) | ((C_u64) random_word()))
+#define	clz64	__builtin_clzll		/* XXX GCCism */
+C_regparm C_word C_fcall
+C_a_i_random_real(C_word **ptr, C_word n) {
+  int exponent = -64;
+  uint64_t significand;
+  unsigned shift;
+
+  while (C_unlikely((significand = random64()) == 0)) {
+    exponent -= 64;
+    if (C_unlikely(exponent < -1074))
+      return 0;
+  }
+
+  shift = clz64(significand);
+  if (shift != 0) {
+    exponent -= shift;
+    significand <<= shift;
+    significand |= (random64() >> (64 - shift));
+  }
+
+  significand |= 1;
+  return C_flonum(ptr, ldexp((double)significand, exponent));
+}
+#undef random64
+
 
 C_word C_set_random_seed(C_word buf, C_word n)
 {