Implement HighPrecision01 distribution

Author: pitdicker

benches/distributions.rs

@@ -55,6 +55,8 @@ distr!(distr_uniform_codepoint, char, Uniform);
 
 distr!(distr_uniform_f32, f32, Uniform);
 distr!(distr_uniform_f64, f64, Uniform);
+distr!(distr_high_precision_f32, f32, HighPrecision01);
+distr!(distr_high_precision_f64, f64, HighPrecision01);
 
 // distributions
 distr!(distr_exp, f64, Exp::new(2.71828 * 3.14159));

src/distributions/float.rs

@@ -14,6 +14,9 @@ use core::mem;
 use Rng;
 use distributions::{Distribution, Uniform};
 
+#[derive(Clone, Copy, Debug)]
+pub struct HighPrecision01;
+
 pub(crate) trait IntoFloat {
     type F;
 
@@ -64,6 +67,64 @@ 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 open interval `(0, 1)`
+            /// (not including either endpoint) with a uniform distribution.
+            ///
+            /// This is different from `Uniform` in that it it uses all 32 bits
+            /// of an RNG for a `f32`, instead of only 23, the number of bits
+            /// that fit in a floats fraction (or 64 instead of 52 bits for a
+            /// `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);
+            /// ```
+            ///
+            /// # 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
+            /// zero's of the random bits.
+            ///
+            /// Each exponent is responsible for a piece of the distribution
+            /// between [0,1). The exponent -1 fills the part [0.5,1). -2 fills
+            /// [0.25,0.5). The lowest exponent we can get is -10. So a problem
+            /// with this method is that we can not fill the part between zero
+            /// and the part from -10. The solution is to treat numbers with an
+            /// exponent of -10 as if they have -9 as exponent, and substract
+            /// 2^-9 (implemented in the `fallback` function).
+            #[inline]
+            fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> $ty {
+                #[inline(never)]
+                fn fallback(fraction: $uty) -> $ty {
+                    let float_size = (mem::size_of::<$ty>() * 8) as i32;
+                    let min_exponent = $fraction_bits as i32 - float_size;
+                    let adjust = // 2^MIN_EXPONENT
+                        (0 as $uty).into_float_with_exponent(min_exponent);
+                    fraction.into_float_with_exponent(min_exponent) - adjust
+                }
+
+                let fraction_mask = (1 << $fraction_bits) - 1;
+                let value = rng.$next_u();
+
+                let fraction = value & fraction_mask;
+                let remaining = value >> $fraction_bits;
+                // If `remaing ==0` we end up in the lowest exponent, which
+                // needs special treatment.
+                if remaining == 0 { return fallback(fraction) }
+                let exp = remaining.trailing_zeros() as i32 + 1;
+                fraction.into_float_with_exponent(-exp)
+            }
+        }
     }
 }
 float_impls! { f32, u32, 23, 127, next_u32 }

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};