Skip to content

Commit

Permalink
Ensure that Uniform<f32/f64> and gen_range<f32/f64> always return val…
Browse files Browse the repository at this point in the history 
…ues in the requested range. Resolves rust-random#476
  • Loading branch information
@sicking
sicking committed May 25, 2018
1 parent 753a76a commit 839ebdb
Showing 1 changed file with 136 additions and 57 deletions.
193 changes: 136 additions & 57 deletions src/distributions/uniform.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,10 +122,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
///
Expand Down Expand Up @@ -508,9 +507,9 @@ 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.
/// Currently there is very little difference between [`new`] and
/// [`new_inclusive`], because the boundaries of a floats range are a bit of a
/// fuzzy concept due to rounding errors.

This comment has been minimized.

Sign in to view

Copy link
@sicking

sicking May 25, 2018

Author Owner

This whole paragraph should likely be removed. Thought I had done that before pushing.
///
/// [`UniformSampler`]: trait.UniformSampler.html
/// [`new`]: trait.UniformSampler.html#tymethod.new
Expand All @@ -519,12 +518,13 @@ wmul_impl_usize! { u64 }
/// [`Standard`]: ../struct.Standard.html
#[derive(Clone, Copy, Debug)]
pub struct UniformFloat<X> {
low: X,
high: X,
scale: X,
offset: X,
}

macro_rules! uniform_float_impl {
($ty:ty, $bits_to_discard:expr, $next_u:ident) => {
($ty:ty, $bits_to_discard:expr, $next_u:ident, $max_fin:expr) => {
impl SampleUniform for $ty {
type Sampler = UniformFloat<$ty>;
}
Expand All @@ -534,59 +534,120 @@ 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 mut scale = high - low;

if scale.is_infinite() {
scale = $max_fin;
}

UniformFloat { low, high, 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 mut adjusted_high;
if high == 0.0 {
// This handles high == -0.0 as well
adjusted_high = <$ty>::from_bits(1);
} else if high > 0.0 {
// If high == $max_fin, this results in
// adjusted_high == infinity. However that's ok since
// high == $max_fin means that all finite numbers are
// inside the range. We just need to be careful when
// calculating `scale` below.
adjusted_high = <$ty>::from_bits(high.to_bits() + 1);
} else {
adjusted_high = <$ty>::from_bits(high.to_bits() - 1);
// XXX this isn't actually needed. But avoiding -0 seems nice.
if adjusted_high == 0.0 {
adjusted_high = 0.0
}
}

debug_assert!(adjusted_high > high ||
(<$ty>::from_bits(high.to_bits() + 1) as f64).is_infinite());

let mut scale = if adjusted_high.is_infinite() {
high - low
} else {
adjusted_high - low
};
if scale.is_infinite() {
scale = $max_fin;
}

UniformFloat { low, high: adjusted_high, 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);
// 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
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;

// 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.
let res = value0_1 * self.scale + self.low;

debug_assert!(self.low <= res);
if res < self.high {
return res;
}
}
}

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, 32 - 23, next_u32, f32::from_bits(0x7f7f_ffff) }
uniform_float_impl! { f64, 64 - 52, next_u64, f64::from_bits(0x7fef_ffff_ffff_ffff) }

/// The back-end implementing [`UniformSampler`] for `Duration`.
///
Expand Down Expand Up @@ -766,24 +827,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]
Expand Down Expand Up @@ -848,7 +926,8 @@ 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.high, 7.0);
assert_eq!(r.inner.scale, 5.0);
}
}

0 comments on commit 839ebdb

Please sign in to comment.