diff --git a/benches/distributions.rs b/benches/distributions.rs index 478619c5112..c443faca2ef 100644 --- a/benches/distributions.rs +++ b/benches/distributions.rs @@ -145,6 +145,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 = XorShiftRng::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(); diff --git a/src/distributions/float.rs b/src/distributions/float.rs index 0058122443c..0e512d2851b 100644 --- a/src/distributions/float.rs +++ b/src/distributions/float.rs @@ -82,8 +82,54 @@ 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)] @@ -91,7 +137,7 @@ macro_rules! float_impls { // 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)] diff --git a/src/distributions/uniform.rs b/src/distributions/uniform.rs index 98bca748d01..05150779615 100644 --- a/src/distributions/uniform.rs +++ b/src/distributions/uniform.rs @@ -113,6 +113,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; @@ -135,10 +139,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 /// @@ -564,10 +567,6 @@ 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 @@ -575,12 +574,12 @@ wmul_impl_usize! { u64 } /// [`Standard`]: ../struct.Standard.html #[derive(Clone, Copy, Debug)] pub struct UniformFloat { + 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>; } @@ -595,12 +594,20 @@ macro_rules! uniform_float_impl { let low = *low_b.borrow(); let high = *high_b.borrow(); 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_b: B1, high_b: B2) -> Self @@ -609,26 +616,40 @@ macro_rules! uniform_float_impl { { let low = *low_b.borrow(); let high = *high_b.borrow(); - 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(&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(low_b: B1, high_b: B2, rng: &mut R) @@ -640,23 +661,42 @@ macro_rules! uniform_float_impl { let high = *high_b.borrow(); 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`. /// @@ -860,24 +900,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] @@ -949,7 +1006,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); } }