From e49be7aae3b33ba9afa1e517022fab3558af1aba Mon Sep 17 00:00:00 2001
From: Huon Wilson <dbau.pp+github@gmail.com>
Date: Sun, 2 Nov 2014 22:47:19 +1100
Subject: [PATCH] rand: Add next_f64/f32 to Rng.

Some random number generates output floating point numbers directly, so
by providing these methods all the functionality in librand is available
with high-performance for these things.

An example of such an is dSFMT (Double precision SIMD-oriented Fast
Mersenne Twister).

The choice to use the open interval [0, 1) has backing elsewhere,
e.g. GSL (GNU Scientific Library) uses this range, and dSFMT supports
generating this natively (I believe the most natural range for that
library is [1, 2), but that is not totally sensible from a user
perspective, and would trip people up).

Fixes https://github.com/rust-lang/rfcs/issues/425.
---
 src/librand/distributions/mod.rs |  7 ++++++
 src/librand/lib.rs               | 40 ++++++++++++++++++++++++++++++++
 src/librand/rand_impls.rs        | 23 +++++++-----------
 3 files changed, 56 insertions(+), 14 deletions(-)

diff --git a/src/librand/distributions/mod.rs b/src/librand/distributions/mod.rs
index 06bd04814c0..66060319891 100644
--- a/src/librand/distributions/mod.rs
+++ b/src/librand/distributions/mod.rs
@@ -227,6 +227,13 @@ fn ziggurat<R:Rng>(
         // creating a f64), so we might as well reuse some to save
         // generating a whole extra random number. (Seems to be 15%
         // faster.)
+        //
+        // This unfortunately misses out on the benefits of direct
+        // floating point generation if an RNG like dSMFT is
+        // used. (That is, such RNGs create floats directly, highly
+        // efficiently and overload next_f32/f64, so by not calling it
+        // this may be slower than it would be otherwise.)
+        // FIXME: investigate/optimise for the above.
         let bits: u64 = rng.gen();
         let i = (bits & 0xff) as uint;
         let f = (bits >> 11) as f64 / SCALE;
diff --git a/src/librand/lib.rs b/src/librand/lib.rs
index 405b70492a3..3c528c564a7 100644
--- a/src/librand/lib.rs
+++ b/src/librand/lib.rs
@@ -78,6 +78,46 @@ pub trait Rng {
         (self.next_u32() as u64 << 32) | (self.next_u32() as u64)
     }
 
+    /// Return the next random f32 selected from the half-open
+    /// interval `[0, 1)`.
+    ///
+    /// By default this is implemented in terms of `next_u32`, but a
+    /// random number generator which can generate numbers satisfying
+    /// the requirements directly can overload this for performance.
+    /// It is required that the return value lies in `[0, 1)`.
+    ///
+    /// See `Closed01` for the closed interval `[0,1]`, and
+    /// `Open01` for the open interval `(0,1)`.
+    fn next_f32(&mut self) -> f32 {
+        const MANTISSA_BITS: uint = 24;
+        const IGNORED_BITS: uint = 8;
+        const SCALE: f32 = (1u64 << MANTISSA_BITS) as f32;
+
+        // using any more than `MANTISSA_BITS` bits will
+        // cause (e.g.) 0xffff_ffff to correspond to 1
+        // exactly, so we need to drop some (8 for f32, 11
+        // for f64) to guarantee the open end.
+        (self.next_u32() >> IGNORED_BITS) as f32 / SCALE
+    }
+
+    /// Return the next random f64 selected from the half-open
+    /// interval `[0, 1)`.
+    ///
+    /// By default this is implemented in terms of `next_u64`, but a
+    /// random number generator which can generate numbers satisfying
+    /// the requirements directly can overload this for performance.
+    /// It is required that the return value lies in `[0, 1)`.
+    ///
+    /// See `Closed01` for the closed interval `[0,1]`, and
+    /// `Open01` for the open interval `(0,1)`.
+    fn next_f64(&mut self) -> f64 {
+        const MANTISSA_BITS: uint = 53;
+        const IGNORED_BITS: uint = 11;
+        const SCALE: f64 = (1u64 << MANTISSA_BITS) as f64;
+
+        (self.next_u64() >> IGNORED_BITS) as f64 / SCALE
+    }
+
     /// Fill `dest` with random data.
     ///
     /// This has a default implementation in terms of `next_u64` and
diff --git a/src/librand/rand_impls.rs b/src/librand/rand_impls.rs
index 77433877ec6..96f40bcc156 100644
--- a/src/librand/rand_impls.rs
+++ b/src/librand/rand_impls.rs
@@ -96,11 +96,11 @@ impl Rand for u64 {
 }
 
 macro_rules! float_impls {
-    ($mod_name:ident, $ty:ty, $mantissa_bits:expr, $method_name:ident, $ignored_bits:expr) => {
+    ($mod_name:ident, $ty:ty, $mantissa_bits:expr, $method_name:ident) => {
         mod $mod_name {
             use {Rand, Rng, Open01, Closed01};
 
-            static SCALE: $ty = (1u64 << $mantissa_bits) as $ty;
+            const SCALE: $ty = (1u64 << $mantissa_bits) as $ty;
 
             impl Rand for $ty {
                 /// Generate a floating point number in the half-open
@@ -110,11 +110,7 @@ macro_rules! float_impls {
                 /// and `Open01` for the open interval `(0,1)`.
                 #[inline]
                 fn rand<R: Rng>(rng: &mut R) -> $ty {
-                    // using any more than `mantissa_bits` bits will
-                    // cause (e.g.) 0xffff_ffff to correspond to 1
-                    // exactly, so we need to drop some (8 for f32, 11
-                    // for f64) to guarantee the open end.
-                    (rng.$method_name() >> $ignored_bits) as $ty / SCALE
+                    rng.$method_name()
                 }
             }
             impl Rand for Open01<$ty> {
@@ -124,23 +120,22 @@ macro_rules! float_impls {
                     // the precision of f64/f32 at 1.0), so that small
                     // numbers are larger than 0, but large numbers
                     // aren't pushed to/above 1.
-                    Open01(((rng.$method_name() >> $ignored_bits) as $ty + 0.25) / SCALE)
+                    Open01(rng.$method_name() + 0.25 / SCALE)
                 }
             }
             impl Rand for Closed01<$ty> {
                 #[inline]
                 fn rand<R: Rng>(rng: &mut R) -> Closed01<$ty> {
-                    // divide by the maximum value of the numerator to
-                    // get a non-zero probability of getting exactly
-                    // 1.0.
-                    Closed01((rng.$method_name() >> $ignored_bits) as $ty / (SCALE - 1.0))
+                    // rescale so that 1.0 - epsilon becomes 1.0
+                    // precisely.
+                    Closed01(rng.$method_name() * SCALE / (SCALE - 1.0))
                 }
             }
         }
     }
 }
-float_impls! { f64_rand_impls, f64, 53, next_u64, 11 }
-float_impls! { f32_rand_impls, f32, 24, next_u32, 8 }
+float_impls! { f64_rand_impls, f64, 53, next_f64 }
+float_impls! { f32_rand_impls, f32, 24, next_f32 }
 
 impl Rand for char {
     #[inline]