Implement HighPrecision01 distribution

Author: pitdicker

benches/distributions.rs

@@ -109,6 +109,8 @@ distr_float!(distr_open01_f32, f32, Open01);
 distr_float!(distr_open01_f64, f64, Open01);
 distr_float!(distr_openclosed01_f32, f32, OpenClosed01);
 distr_float!(distr_openclosed01_f64, f64, OpenClosed01);
+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

@@ -97,6 +97,35 @@ impl<F: HPFloatHelper> HighPrecision<F> {
     }
 }
 
+/// Generate a floating point number in the half-open interval `[0, 1)` with a
+/// uniform distribution.
+///
+/// This is different from `Uniform` in that 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`).
+///
+/// The smallest interval between values that can be generated is 2^-32
+/// (2.3283064e-10) for `f32`, and 2^-64 (5.421010862427522e-20) for `f64`.
+/// But this interval increases further away from zero because of limitations of
+/// the floating point format. Close to 1.0 the interval is 2^-24 (5.9604645e-8)
+/// for `f32`, and 2^-53 (1.1102230246251565) for `f64`. Compare this with
+/// `Uniform`, which has a fixed interval of 2^23 and 2^-52 respectively.
+///
+/// Note: in the future this may change change to request even more bits from
+/// the RNG if the value gets very close to 0.0, so it always has as many digits
+/// of precision as the float can represent.
+///
+/// # 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;
 
@@ -203,6 +232,64 @@ float_impls! { f64x8, u64x8, f64, u64, 52, 1023 }
 
 macro_rules! high_precision_float_impls {
     ($ty:ty, $uty:ty, $ity:ty, $fraction_bits:expr, $exponent_bits:expr) => {
+        impl Distribution<$ty> for HighPrecision01 {
+            /// Generate a floating point number in the half-open interval
+            /// `[0, 1)` with a uniform distribution.
+            ///
+            /// This is different from `Uniform` in that 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: $uty = rng.gen();
+
+                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)
+            }
+        }
+
         impl HPFloatHelper for $ty {
             type SignedInt = $ity;
             type SignedIntDistribution = <$ity as SampleUniform>::Sampler;
@@ -630,4 +717,18 @@ mod tests {
                  ]);
         }
     }
+
+    #[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);
+    }
 }

src/distributions/mod.rs

@@ -175,7 +175,7 @@ use Rng;
 
 #[doc(inline)] pub use self::other::Alphanumeric;
 #[doc(inline)] pub use self::uniform::Uniform;
-#[doc(inline)] pub use self::float::{OpenClosed01, Open01, HighPrecision};
+#[doc(inline)] pub use self::float::{OpenClosed01, Open01, HighPrecision, HighPrecision01};
 #[cfg(feature="alloc")]
 #[doc(inline)] pub use self::weighted::WeightedIndex;
 #[cfg(feature="std")]