Ensure that Uniform<f32/f64> and gen_range<f32/f64> always return val…

…ues in the requested range. Resolves #476

benches/distributions.rs

@@ -10,7 +10,7 @@ const RAND_BENCH_N: u64 = 1000;
 use std::mem::size_of;
 use test::Bencher;
 
-use rand::{Rng, FromEntropy, XorShiftRng};
+use rand::{Rng, FromEntropy, XorShiftRng, StdRng};
 use rand::distributions::*;
 
 macro_rules! distr_int {
@@ -144,6 +144,33 @@ gen_range_int!(gen_range_i64, i64, 3i64, 123_456_789_123);
 #[cfg(feature = "i128_support")]
 gen_range_int!(gen_range_i128, i128, -12345678901234i128, 123_456_789_123_456_789);
 
+// construct and sample from a floating-point range
+macro_rules! gen_range_float {
+    ($fnn:ident, $ty:ident, $low:expr, $high:expr) => {
+        #[bench]
+        fn $fnn(b: &mut Bencher) {
+            let mut rng = StdRng::from_entropy();
+
+            b.iter(|| {
+                let mut high = $high;
+                let mut low = $low;
+                let mut accum: $ty = 0.0;
+                for _ in 0..::RAND_BENCH_N {
+                    accum += rng.gen_range(low, high);
+                    // force recalculation of range each time
+                    low += 0.9;
+                    high += 1.1;
+                }
+                accum
+            });
+            b.bytes = size_of::<$ty>() as u64 * ::RAND_BENCH_N;
+        }
+    }
+}
+
+gen_range_float!(gen_range_f32, f32, -20000.0f32, 100000.0);
+gen_range_float!(gen_range_f64, f64, 123.456f64, 7890.12);
+
 #[bench]
 fn dist_iter(b: &mut Bencher) {
     let mut rng = XorShiftRng::from_entropy();

src/distributions/float.rs

@@ -82,16 +82,62 @@ pub(crate) trait IntoFloat {
     fn into_float_with_exponent(self, exponent: i32) -> Self::F;
 }
 
+#[cfg(not(std))]
+pub(crate) trait Float : Sized {
+    type Bits;
+
+    fn is_nan(self) -> bool;
+
+    fn is_infinite(self) -> bool;
+
+    fn is_finite(self) -> bool;
+
+    fn to_bits(self) -> Self::Bits;
+
+    fn from_bits(v: Self::Bits) -> Self;
+}
+
 macro_rules! float_impls {
-    ($ty:ty, $uty:ty, $fraction_bits:expr, $exponent_bias:expr) => {
+    ($ty:ty, $ty_ident:ident, $uty:ty, $fraction_bits:expr, $exponent_bias:expr) => {
+        #[cfg(not(std))]
+        impl Float for $ty {
+            type Bits = $uty;
+
+            #[inline]
+            fn is_nan(self) -> bool {
+                self != self
+            }
+
+            #[inline]
+            fn is_infinite(self) -> bool {
+                self == ::core::$ty_ident::INFINITY || self == ::core::$ty_ident::NEG_INFINITY
+            }
+
+            #[inline]
+            fn is_finite(self) -> bool {
+                !(self.is_nan() || self.is_infinite())
+            }
+
+            #[inline]
+            fn to_bits(self) -> Self::Bits {
+                unsafe { mem::transmute(self) }
+            }
+
+            #[inline]
+            fn from_bits(v: Self::Bits) -> Self {
+                // It turns out the safety issues with sNaN were overblown! Hooray!
+                unsafe { mem::transmute(v) }
+            }
+        }
+
         impl IntoFloat for $uty {
             type F = $ty;
             #[inline(always)]
             fn into_float_with_exponent(self, exponent: i32) -> $ty {
                 // The exponent is encoded using an offset-binary representation
                 let exponent_bits =
                     (($exponent_bias + exponent) as $uty) << $fraction_bits;
-                unsafe { mem::transmute(self | exponent_bits) }
+                <$ty>::from_bits(self | exponent_bits)
             }
         }
 
@@ -140,8 +186,8 @@ macro_rules! float_impls {
         }
     }
 }
-float_impls! { f32, u32, 23, 127 }
-float_impls! { f64, u64, 52, 1023 }
+float_impls! { f32, f32, u32, 23, 127 }
+float_impls! { f64, f64, u64, 52, 1023 }
 
 
 #[cfg(test)]

src/distributions/uniform.rs

@@ -100,6 +100,10 @@
 #[cfg(feature = "std")]
 use std::time::Duration;
 
+#[cfg(not(feature = "std"))]
+#[allow(unused_imports)] // rustc doesn't detect that this is actually used
+use distributions::float::Float;
+
 use Rng;
 use distributions::Distribution;
 use distributions::float::IntoFloat;
@@ -122,10 +126,9 @@ use distributions::float::IntoFloat;
 /// generated by the RNG than the low bits, since with some RNGs the low-bits
 /// are of lower quality than the high bits.
 ///
-/// Implementations should attempt to sample in `[low, high)` for
-/// `Uniform::new(low, high)`, i.e., excluding `high`, but this may be very
-/// difficult. All the primitive integer types satisfy this property, and the
-/// float types normally satisfy it, but rounding may mean `high` can occur.
+/// Implementations must sample in `[low, high)` range for
+/// `Uniform::new(low, high)`, i.e., excluding `high`. In particular care must
+/// be taken to ensure that rounding never results values `< low` or `>= high`.
 ///
 /// # Example
 ///
@@ -508,23 +511,19 @@ wmul_impl_usize! { u64 }
 /// multiply and addition. Values produced this way have what equals 22 bits of
 /// random digits for an `f32`, and 52 for an `f64`.
 ///
-/// Currently there is no difference between [`new`] and [`new_inclusive`],
-/// because the boundaries of a floats range are a bit of a fuzzy concept due to
-/// rounding errors.
-///
 /// [`UniformSampler`]: trait.UniformSampler.html
 /// [`new`]: trait.UniformSampler.html#tymethod.new
 /// [`new_inclusive`]: trait.UniformSampler.html#tymethod.new_inclusive
 /// [`Uniform`]: struct.Uniform.html
 /// [`Standard`]: ../struct.Standard.html
 #[derive(Clone, Copy, Debug)]
 pub struct UniformFloat<X> {
+    low: X,
     scale: X,
-    offset: X,
 }
 
 macro_rules! uniform_float_impl {
-    ($ty:ty, $bits_to_discard:expr, $next_u:ident) => {
+    ($ty:ty, $uty:ident, $bits_to_discard:expr, $next_u:ident, $max_fin:expr) => {
         impl SampleUniform for $ty {
             type Sampler = UniformFloat<$ty>;
         }
@@ -534,59 +533,100 @@ macro_rules! uniform_float_impl {
 
             fn new(low: Self::X, high: Self::X) -> Self {
                 assert!(low < high, "Uniform::new called with `low >= high`");
-                let scale = high - low;
-                let offset = low - scale;
-                UniformFloat {
-                    scale: scale,
-                    offset: offset,
+                assert!(low.is_finite() && high.is_finite(), "Uniform::new called with non-finite boundaries");
+
+                let max_rand = (::core::$uty::MAX >> $bits_to_discard)
+                               .into_float_with_exponent(0) - 1.0;
+
+                let mut scale = high - low;
+
+                while scale * max_rand + low >= high {
+                    scale = <$ty>::from_bits(scale.to_bits() - 1);
                 }
+
+                debug_assert!(scale >= 0.0);
+
+                UniformFloat { low, scale }
             }
 
             fn new_inclusive(low: Self::X, high: Self::X) -> Self {
-                assert!(low <= high,
-                        "Uniform::new_inclusive called with `low > high`");
-                let scale = high - low;
-                let offset = low - scale;
-                UniformFloat {
-                    scale: scale,
-                    offset: offset,
+                assert!(low <= high, "Uniform::new_inclusive called with `low > high`");
+                assert!(low.is_finite() && high.is_finite(), "Uniform::new_inclusive called with non-finite boundaries");
+
+                let max_rand = (::core::$uty::MAX >> $bits_to_discard)
+                               .into_float_with_exponent(0) - 1.0;
+
+                let mut scale = (high - low) / max_rand;
+
+                // Adjust scale lower until the max value we can't generate
+                // values above `high`
+                while scale * max_rand + low > high {
+                    scale = <$ty>::from_bits(scale.to_bits() - 1);
                 }
+
+                debug_assert!(scale >= 0.0);
+
+                UniformFloat { low, scale }
             }
 
             fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::X {
                 // Generate a value in the range [1, 2)
                 let value1_2 = (rng.$next_u() >> $bits_to_discard)
                                .into_float_with_exponent(0);
+
+                // Get a value in the range [0, 1) in order to avoid
+                // overflowing into infinity when multiplying with scale
+                let value0_1 = value1_2 - 1.0;
+
                 // We don't use `f64::mul_add`, because it is not available with
                 // `no_std`. Furthermore, it is slower for some targets (but
                 // faster for others). However, the order of multiplication and
                 // addition is important, because on some platforms (e.g. ARM)
                 // it will be optimized to a single (non-FMA) instruction.
-                value1_2 * self.scale + self.offset
+                value0_1 * self.scale + self.low
             }
 
             fn sample_single<R: Rng + ?Sized>(low: Self::X,
                                               high: Self::X,
                                               rng: &mut R) -> Self::X {
                 assert!(low < high,
                         "Uniform::sample_single called with low >= high");
-                let scale = high - low;
-                let offset = low - scale;
-                // Generate a value in the range [1, 2)
-                let value1_2 = (rng.$next_u() >> $bits_to_discard)
-                               .into_float_with_exponent(0);
-                // Doing multiply before addition allows some architectures to
-                // use a single instruction.
-                value1_2 * scale + offset
+                let mut scale = high - low;
+
+                loop {
+                    // Generate a value in the range [1, 2)
+                    let value1_2 = (rng.$next_u() >> $bits_to_discard)
+                                   .into_float_with_exponent(0);
+
+                    // Get a value in the range [0, 1) in order to avoid
+                    // overflowing into infinity when multiplying with scale
+                    let value0_1 = value1_2 - 1.0;
+
+                    // Doing multiply before addition allows some architectures
+                    // to use a single instruction.
+                    let res = value0_1 * scale + low;
+
+                    debug_assert!(low <= res);
+                    if res < high {
+                        return res;
+                    }
+
+                    // This technically should be outside and before the loop.
+                    // If `scale` is infinite then the first iteration through
+                    // the loop is guaranteed to fail. However the common case
+                    // is that scale is finite and the first iteration
+                    // succeeds, so optimize for that case.
+                    if scale.is_infinite() {
+                        scale = $max_fin;
+                    }
+                }
             }
         }
     }
 }
 
-uniform_float_impl! { f32, 32 - 23, next_u32 }
-uniform_float_impl! { f64, 64 - 52, next_u64 }
-
-
+uniform_float_impl! { f32, u32, 32 - 23, next_u32, f32::from_bits(0x7f7f_ffff) }
+uniform_float_impl! { f64, u64, 64 - 52, next_u64, f64::from_bits(0x7fef_ffff_ffff_ffff) }
 
 /// The back-end implementing [`UniformSampler`] for `Duration`.
 ///
@@ -768,24 +808,41 @@ mod tests {
     fn test_floats() {
         let mut rng = ::test::rng(252);
         macro_rules! t {
-            ($($ty:ty),*) => {{
-                $(
-                   let v: &[($ty, $ty)] = &[(0.0, 100.0),
-                                            (-1e35, -1e25),
-                                            (1e-35, 1e-25),
-                                            (-1e35, 1e35)];
-                   for &(low, high) in v.iter() {
-                        let my_uniform = Uniform::new(low, high);
-                        for _ in 0..1000 {
-                            let v: $ty = rng.sample(my_uniform);
-                            assert!(low <= v && v < high);
-                        }
+            ($ty:ty, $max_fin:expr) => {{
+                let v: &[($ty, $ty)] = &[(0.0, 100.0),
+                                         (-1e35, -1e25),
+                                         (1e-35, 1e-25),
+                                         (-1e35, 1e35),
+                                         (<$ty>::from_bits(0), <$ty>::from_bits(3)),
+                                         (-<$ty>::from_bits(10), -<$ty>::from_bits(1)),
+                                         (-<$ty>::from_bits(5), 0.0),
+                                         (-<$ty>::from_bits(7), -0.0),
+                                         (10.0, $max_fin),
+                                         (-10.0, $max_fin),
+                                         (-$max_fin, $max_fin),
+                                        ];
+                for &(low, high) in v.iter() {
+                    let my_uniform = Uniform::new(low, high);
+                    let my_incl_uniform = Uniform::new_inclusive(low, high);
+                    for _ in 0..1000 {
+                        let v = rng.sample(my_uniform);
+                        assert!(low <= v && v < high);
+                        let v = rng.sample(my_incl_uniform);
+                        assert!(low <= v && v <= high);
+                        let v = rng.gen_range(low, high);
+                        assert!(low <= v && v < high);
                     }
-                 )*
+                }
+
+                let edge_case_inclusive = Uniform::new_inclusive($max_fin, $max_fin);
+                let v = rng.sample(edge_case_inclusive);
+                assert_eq!(v, $max_fin);
+
             }}
         }
 
-        t!(f32, f64)
+        t!(f32, f32::from_bits(0x7f7f_ffff));
+        t!(f64, f64::from_bits(0x7fef_ffff_ffff_ffff));
     }
 
     #[test]
@@ -850,7 +907,7 @@ mod tests {
         assert_eq!(r.inner.low, 2);
         assert_eq!(r.inner.range, 5);
         let r = Uniform::from(2.0f64..7.0);
-        assert_eq!(r.inner.offset, -3.0);
+        assert_eq!(r.inner.low, 2.0);
         assert_eq!(r.inner.scale, 5.0);
     }
 }