Implement HighPrecision01 distribution

benches/distributions.rs

@@ -95,6 +95,8 @@ distr!(distr_uniform_codepoint, char, Standard);
 
 distr_float!(distr_uniform_f32, f32, Standard);
 distr_float!(distr_uniform_f64, f64, Standard);
+distr_float!(distr_high_precision_f32, f32, HighPrecision01);
+distr_float!(distr_high_precision_f64, f64, HighPrecision01);
 
 // distributions
 distr_float!(distr_exp, f64, Exp::new(1.23 * 4.56));

src/distributions/float.rs

@@ -10,10 +10,37 @@
 
 //! Basic floating-point number distributions
 
-use core::mem;
+use core::{cmp, mem};
 use Rng;
 use distributions::{Distribution, Standard};
 
+/// Generate a floating point number in the half-open interval `[0, 1)` with a
+/// uniform distribution, with as much precision as the floating-point type
+/// can represent, including sub-normals.
+///
+/// Technically 0 is representable, but the probability of occurrence is
+/// remote (1 in 2^149 for `f32` or 1 in 2^1074 for `f64`).
+///
+/// This is different from `Uniform` in that it uses as many random bits as
+/// required to get high precision close to 0. Normally only a single call to
+/// the source RNG is required (32 bits for `f32` or 64 bits for `f64`); 1 in
+/// 2^9 (`f32`) or 2^12 (`f64`) samples need an extra call; of these 1 in 2^32
+/// or 1 in 2^64 require a third call, etc.; i.e. even for `f32` a third call is
+/// almost impossible to observe with an unbiased RNG. Due to the extra logic
+/// there is some performance overhead relative to `Uniform`; this is more
+/// significant for `f32` than for `f64`.
+///
+/// # Example
+/// ```rust
+/// use rand::{NewRng, SmallRng, Rng};
+/// use rand::distributions::HighPrecision01;
+///
+/// let val: f32 = SmallRng::new().sample(HighPrecision01);
+/// println!("f32 from [0,1): {}", val);
+/// ```
+#[derive(Clone, Copy, Debug)]
+pub struct HighPrecision01;
+
 pub(crate) trait IntoFloat {
     type F;
 
@@ -54,6 +81,66 @@ macro_rules! float_impls {
                 fraction.into_float_with_exponent(0) - (1.0 - EPSILON / 2.0)
             }
         }
+
+        impl Distribution<$ty> for HighPrecision01 {
+            /// Generate a floating point number in the half-open interval
+            /// `[0, 1)` with a uniform distribution. See [`HighPrecision01`].
+            ///
+            /// # Algorithm
+            /// (Note: this description used values that apply to `f32` to
+            /// illustrate the algorithm).
+            ///
+            /// The trick to generate a uniform distribution over [0,1) is to
+            /// set the exponent to the -log2 of the remaining random bits. A
+            /// simpler alternative to -log2 is to count the number of trailing
+            /// zeros in the random bits. In the case where all bits are zero,
+            /// we simply generate a new random number and add the number of
+            /// trailing zeros to the previous count (up to maximum exponent).
+            ///
+            /// Each exponent is responsible for a piece of the distribution
+            /// between [0,1). We take the above exponent, add 1 and negate;
+            /// thus with probability 1/2 we have exponent -1 which fills the
+            /// range [0.5,1); with probability 1/4 we have exponent -2 which
+            /// fills the range [0.25,0.5), etc. If the exponent reaches the
+            /// minimum allowed, the floating-point format drops the implied
+            /// fraction bit, thus allowing numbers down to 0 to be sampled.
+            ///
+            /// [`HighPrecision01`]: struct.HighPrecision01.html
+            #[inline]
+            fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> $ty {
+                // Unusual case. Separate function to allow inlining of rest.
+                #[inline(never)]
+                fn fallback<R: Rng + ?Sized>(mut exp: i32, fraction: $uty, rng: &mut R) -> $ty {
+                    // Performance impact of code here is negligible.
+                    let bits = rng.gen::<$uty>();
+                    exp += bits.trailing_zeros() as i32;
+                    // If RNG were guaranteed unbiased we could skip the
+                    // check against exp; unfortunately it may be.
+                    // Worst case ("zeros" RNG) has recursion depth 16.
+                    if bits == 0 && exp < $exponent_bias {
+                        return fallback(exp, fraction, rng);
+                    }
+                    exp = cmp::min(exp, $exponent_bias);
+                    fraction.into_float_with_exponent(-exp)
+                }
+
+                let fraction_mask = (1 << $fraction_bits) - 1;
+                let value = rng.$next_u();
+
+                let fraction = value & fraction_mask;
+                let remaining = value >> $fraction_bits;
+                if remaining == 0 {
+                    // exp is compile-time constant so this reduces to a function call:
+                    let size_bits = (mem::size_of::<$ty>() * 8) as i32;
+                    let exp = (size_bits - $fraction_bits as i32) + 1;
+                    return fallback(exp, fraction, rng);
+                }
+
+                // Usual case: exponent from -1 to -9 (f32) or -12 (f64)
+                let exp = remaining.trailing_zeros() as i32 + 1;
+                fraction.into_float_with_exponent(-exp)
+            }
+        }
     }
 }
 float_impls! { f32, u32, 23, 127, next_u32 }
@@ -62,7 +149,8 @@ float_impls! { f64, u64, 52, 1023, next_u64 }
 
 #[cfg(test)]
 mod tests {
-    use Rng;
+    use {Rng};
+    use distributions::HighPrecision01;
     use mock::StepRng;
 
     const EPSILON32: f32 = ::core::f32::EPSILON;
@@ -86,4 +174,45 @@ mod tests {
         assert_eq!(max.gen::<f32>(), 1.0 - EPSILON32 / 2.0);
         assert_eq!(max.gen::<f64>(), 1.0 - EPSILON64 / 2.0);
     }
+
+    #[test]
+    fn high_precision_01_edge_cases() {
+        // Test that the distribution is a half-open range over [0,1).
+        // These constants happen to generate the lowest and highest floats in
+        // the range.
+        let mut zeros = StepRng::new(0, 0);
+        assert_eq!(zeros.sample::<f32, _>(HighPrecision01), 0.0);
+        assert_eq!(zeros.sample::<f64, _>(HighPrecision01), 0.0);
+
+        let mut ones = StepRng::new(0xffff_ffff_ffff_ffff, 0);
+        assert_eq!(ones.sample::<f32, _>(HighPrecision01), 0.99999994);
+        assert_eq!(ones.sample::<f64, _>(HighPrecision01), 0.9999999999999999);
+    }
+
+    #[cfg(feature="std")] mod mean {
+        use {Rng, SmallRng, SeedableRng, thread_rng};
+        use distributions::{Uniform, HighPrecision01};
+
+        macro_rules! test_mean {
+            ($name:ident, $ty:ty, $distr:expr) => {
+        #[test]
+        fn $name() {
+            // TODO: no need to &mut here:
+            let mut rng = SmallRng::from_rng(&mut thread_rng()).unwrap();
+            let mut total: $ty = 0.0;
+            const N: u32 = 1_000_000;
+            for _ in 0..N {
+                total += rng.sample::<$ty, _>($distr);
+            }
+            let avg = total / (N as $ty);
+            //println!("average over {} samples: {}", N, avg);
+            assert!(0.499 < avg && avg < 0.501);
+        }
+        } }
+
+        test_mean!(test_mean_f32, f32, Uniform);
+        test_mean!(test_mean_f64, f64, Uniform);
+        test_mean!(test_mean_high_f32, f32, HighPrecision01);
+        test_mean!(test_mean_high_f64, f64, HighPrecision01);
+    }
 }

src/distributions/mod.rs

@@ -18,6 +18,7 @@
 use Rng;
 
 pub use self::other::Alphanumeric;
+pub use self::float::HighPrecision01;
 pub use self::range::Range;
 #[cfg(feature="std")]
 pub use self::gamma::{Gamma, ChiSquared, FisherF, StudentT};