diff --git a/benches/distributions.rs b/benches/distributions.rs index 568d5024992..2ab624ae70a 100644 --- a/benches/distributions.rs +++ b/benches/distributions.rs @@ -38,8 +38,8 @@ distr!(distr_range_i64, i64, Range::new(3i64, 12345678901234)); #[cfg(feature = "i128_support")] distr!(distr_range_i128, i128, Range::new(-12345678901234i128, 12345678901234567890)); -distr!(distr_range_float32, f32, Range::new(2.26f32, 2.319)); -distr!(distr_range_float, f64, Range::new(2.26f64, 2.319)); +distr!(distr_range_f32, f32, Range::new(2.26f32, 2.319)); +distr!(distr_range_f64, f64, Range::new(2.26f64, 2.319)); // uniform distr!(distr_uniform_i8, i8, Uniform); @@ -52,13 +52,8 @@ distr!(distr_uniform_i128, i128, Uniform); distr!(distr_uniform_bool, bool, Uniform); distr!(distr_uniform_codepoint, char, Uniform); -distr!(distr_uniform01_float32, f32, Uniform); -distr!(distr_closed01_float32, f32, Closed01); -distr!(distr_open01_float32, f32, Open01); - -distr!(distr_uniform01_float, f64, Uniform); -distr!(distr_closed01_float, f64, Closed01); -distr!(distr_open01_float, f64, Open01); +distr!(distr_uniform_f32, f32, Uniform); +distr!(distr_uniform_f64, f64, Uniform); // distributions distr!(distr_exp, f64, Exp::new(2.71828 * 3.14159)); diff --git a/src/distributions/float.rs b/src/distributions/float.rs index 4f927fbd164..09b27e409bc 100644 --- a/src/distributions/float.rs +++ b/src/distributions/float.rs @@ -14,125 +14,66 @@ use core::mem; use Rng; use distributions::{Distribution, Uniform}; - -/// A distribution to sample floating point numbers uniformly in the open -/// interval `(0, 1)` (not including either endpoint). -/// -/// See also: [`Closed01`] for the closed `[0, 1]`; [`Uniform`] for the -/// half-open `[0, 1)`. -/// -/// # Example -/// ```rust -/// use rand::{weak_rng, Rng}; -/// use rand::distributions::Open01; -/// -/// let val: f32 = weak_rng().sample(Open01); -/// println!("f32 from (0,1): {}", val); -/// ``` -/// -/// [`Uniform`]: struct.Uniform.html -/// [`Closed01`]: struct.Closed01.html -#[derive(Clone, Copy, Debug)] -pub struct Open01; - -/// A distribution to sample floating point numbers uniformly in the closed -/// interval `[0, 1]` (including both endpoints). -/// -/// See also: [`Open01`] for the open `(0, 1)`; [`Uniform`] for the half-open -/// `[0, 1)`. -/// -/// # Example -/// ```rust -/// use rand::{weak_rng, Rng}; -/// use rand::distributions::Closed01; -/// -/// let val: f32 = weak_rng().sample(Closed01); -/// println!("f32 from [0,1]: {}", val); -/// ``` -/// -/// [`Uniform`]: struct.Uniform.html -/// [`Open01`]: struct.Open01.html -#[derive(Clone, Copy, Debug)] -pub struct Closed01; - - -// Return the next random f32 selected from the half-open -// interval `[0, 1)`. -// -// This uses a technique described by Saito and Matsumoto at -// MCQMC'08. Given that the IEEE floating point numbers are -// uniformly distributed over [1,2), we generate a number in -// this range and then offset it onto the range [0,1). Our -// choice of bits (masking v. shifting) is arbitrary and -// should be immaterial for high quality generators. For low -// quality generators (ex. LCG), prefer bitshifting due to -// correlation between sequential low order bits. -// -// See: -// A PRNG specialized in double precision floating point numbers using -// an affine transition -// -// * -// * -impl Distribution for Uniform { - fn sample(&self, rng: &mut R) -> f32 { - const UPPER_MASK: u32 = 0x3F800000; - const LOWER_MASK: u32 = 0x7FFFFF; - let tmp = UPPER_MASK | (rng.next_u32() & LOWER_MASK); - let result: f32 = unsafe { mem::transmute(tmp) }; - result - 1.0 - } -} -impl Distribution for Uniform { - fn sample(&self, rng: &mut R) -> f64 { - const UPPER_MASK: u64 = 0x3FF0000000000000; - const LOWER_MASK: u64 = 0xFFFFFFFFFFFFF; - let tmp = UPPER_MASK | (rng.next_u64() & LOWER_MASK); - let result: f64 = unsafe { mem::transmute(tmp) }; - result - 1.0 - } +pub(crate) trait IntoFloat { + type F; + + /// Helper method to combine the fraction and a contant exponent into a + /// float. + /// + /// Only the least significant bits of `self` may be set, 23 for `f32` and + /// 52 for `f64`. + /// The resulting value will fall in a range that depends on the exponent. + /// As an example the range with exponent 0 will be + /// [20..21), which is [1..2). + #[inline(always)] + fn into_float_with_exponent(self, exponent: i32) -> Self::F; } macro_rules! float_impls { - ($mod_name:ident, $ty:ty, $mantissa_bits:expr) => { - mod $mod_name { - use Rng; - use distributions::{Distribution}; - use super::{Open01, Closed01}; - - const SCALE: $ty = (1u64 << $mantissa_bits) as $ty; - - impl Distribution<$ty> for Open01 { - #[inline] - fn sample(&self, rng: &mut R) -> $ty { - // add 0.5 * epsilon, so that smallest number is - // greater than 0, and largest number is still - // less than 1, specifically 1 - 0.5 * epsilon. - let x: $ty = rng.gen(); - x + 0.5 / SCALE - } + ($ty:ty, $uty:ty, $fraction_bits:expr, $exponent_bias:expr, + $next_u:ident) => { + 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) } } - impl Distribution<$ty> for Closed01 { - #[inline] - fn sample(&self, rng: &mut R) -> $ty { - // rescale so that 1.0 - epsilon becomes 1.0 - // precisely. - let x: $ty = rng.gen(); - x * SCALE / (SCALE - 1.0) - } + } + + impl Distribution<$ty> for Uniform { + /// Generate a floating point number in the open interval `(0, 1)` + /// (not including either endpoint) with a uniform distribution. + /// + /// # Example + /// ```rust + /// use rand::{weak_rng, Rng}; + /// use rand::distributions::Uniform; + /// + /// let val: f32 = weak_rng().sample(Uniform); + /// println!("f32 from (0,1): {}", val); + /// ``` + fn sample(&self, rng: &mut R) -> $ty { + const EPSILON: $ty = 1.0 / (1u64 << $fraction_bits) as $ty; + let float_size = mem::size_of::<$ty>() * 8; + + let value = rng.$next_u(); + let fraction = value >> (float_size - $fraction_bits); + fraction.into_float_with_exponent(0) - (1.0 - EPSILON / 2.0) } } } } -float_impls! { f64_rand_impls, f64, 52 } -float_impls! { f32_rand_impls, f32, 23 } +float_impls! { f32, u32, 23, 127, next_u32 } +float_impls! { f64, u64, 52, 1023, next_u64 } #[cfg(test)] mod tests { use Rng; use mock::StepRng; - use distributions::{Open01, Closed01}; const EPSILON32: f32 = ::core::f32::EPSILON; const EPSILON64: f64 = ::core::f64::EPSILON; @@ -140,77 +81,19 @@ mod tests { #[test] fn floating_point_edge_cases() { let mut zeros = StepRng::new(0, 0); - assert_eq!(zeros.gen::(), 0.0); - assert_eq!(zeros.gen::(), 0.0); - - let mut one = StepRng::new(1, 0); - assert_eq!(one.gen::(), EPSILON32); - assert_eq!(one.gen::(), EPSILON64); - - let mut max = StepRng::new(!0, 0); - assert_eq!(max.gen::(), 1.0 - EPSILON32); - assert_eq!(max.gen::(), 1.0 - EPSILON64); - } + assert_eq!(zeros.gen::(), 0.0 + EPSILON32 / 2.0); + assert_eq!(zeros.gen::(), 0.0 + EPSILON64 / 2.0); - #[test] - fn fp_closed_edge_cases() { - let mut zeros = StepRng::new(0, 0); - assert_eq!(zeros.sample::(Closed01), 0.0); - assert_eq!(zeros.sample::(Closed01), 0.0); - - let mut one = StepRng::new(1, 0); - let one32 = one.sample::(Closed01); - let one64 = one.sample::(Closed01); - assert!(EPSILON32 < one32 && one32 < EPSILON32 * 1.01); - assert!(EPSILON64 < one64 && one64 < EPSILON64 * 1.01); - - let mut max = StepRng::new(!0, 0); - assert_eq!(max.sample::(Closed01), 1.0); - assert_eq!(max.sample::(Closed01), 1.0); - } - - #[test] - fn fp_open_edge_cases() { - let mut zeros = StepRng::new(0, 0); - assert_eq!(zeros.sample::(Open01), 0.0 + EPSILON32 / 2.0); - assert_eq!(zeros.sample::(Open01), 0.0 + EPSILON64 / 2.0); - - let mut one = StepRng::new(1, 0); - let one32 = one.sample::(Open01); - let one64 = one.sample::(Open01); + let mut one = StepRng::new(1 << 9, 0); + let one32 = one.gen::(); assert!(EPSILON32 < one32 && one32 < EPSILON32 * 2.0); - assert!(EPSILON64 < one64 && one64 < EPSILON64 * 2.0); - - let mut max = StepRng::new(!0, 0); - assert_eq!(max.sample::(Open01), 1.0 - EPSILON32 / 2.0); - assert_eq!(max.sample::(Open01), 1.0 - EPSILON64 / 2.0); - } - #[test] - fn rand_open() { - // this is unlikely to catch an incorrect implementation that - // generates exactly 0 or 1, but it keeps it sane. - let mut rng = ::test::rng(510); - for _ in 0..1_000 { - // strict inequalities - let f: f64 = rng.sample(Open01); - assert!(0.0 < f && f < 1.0); - - let f: f32 = rng.sample(Open01); - assert!(0.0 < f && f < 1.0); - } - } - - #[test] - fn rand_closed() { - let mut rng = ::test::rng(511); - for _ in 0..1_000 { - // strict inequalities - let f: f64 = rng.sample(Closed01); - assert!(0.0 <= f && f <= 1.0); + let mut one = StepRng::new(1 << 12, 0); + let one64 = one.gen::(); + assert!(EPSILON64 < one64 && one64 < EPSILON64 * 2.0); - let f: f32 = rng.sample(Closed01); - assert!(0.0 <= f && f <= 1.0); - } + let mut max = StepRng::new(!0, 0); + assert_eq!(max.gen::(), 1.0 - EPSILON32 / 2.0); + assert_eq!(max.gen::(), 1.0 - EPSILON64 / 2.0); } } diff --git a/src/distributions/gamma.rs b/src/distributions/gamma.rs index 361e12ed6a2..7522b2fa117 100644 --- a/src/distributions/gamma.rs +++ b/src/distributions/gamma.rs @@ -17,7 +17,7 @@ use self::ChiSquaredRepr::*; use {Rng}; use distributions::normal::StandardNormal; -use distributions::{Distribution, Exp, Open01}; +use distributions::{Distribution, Exp, Uniform}; /// The Gamma distribution `Gamma(shape, scale)` distribution. /// @@ -144,7 +144,7 @@ impl Distribution for Gamma { } impl Distribution for GammaSmallShape { fn sample(&self, rng: &mut R) -> f64 { - let u: f64 = rng.sample(Open01); + let u: f64 = rng.sample(Uniform); self.large_shape.sample(rng) * u.powf(self.inv_shape) } @@ -159,7 +159,7 @@ impl Distribution for GammaLargeShape { } let v = v_cbrt * v_cbrt * v_cbrt; - let u: f64 = rng.sample(Open01); + let u: f64 = rng.sample(Uniform); let x_sqr = x * x; if u < 1.0 - 0.0331 * x_sqr * x_sqr || diff --git a/src/distributions/mod.rs b/src/distributions/mod.rs index 3ff5700a376..0e8d3237eb3 100644 --- a/src/distributions/mod.rs +++ b/src/distributions/mod.rs @@ -17,7 +17,6 @@ use Rng; -pub use self::float::{Open01, Closed01}; pub use self::range::Range; #[cfg(feature="std")] pub use self::gamma::{Gamma, ChiSquared, FisherF, StudentT}; @@ -39,6 +38,8 @@ mod integer; mod other; #[cfg(feature="std")] mod ziggurat_tables; +#[cfg(feature="std")] +use distributions::float::IntoFloat; /// Types that can be used to create a random instance of `Support`. #[deprecated(since="0.5.0", note="use Distribution instead")] @@ -53,9 +54,6 @@ pub trait Sample { /// Since no state is recorded, each sample is (statistically) /// independent of all others, assuming the `Rng` used has this /// property. -// FIXME maybe having this separate is overkill (the only reason is to -// take &self rather than &mut self)? or maybe this should be the -// trait called `Sample` and the other should be `DependentSample`. #[allow(deprecated)] #[deprecated(since="0.5.0", note="use Distribution instead")] pub trait IndependentSample: Sample { @@ -74,7 +72,7 @@ mod impls { use distributions::gamma::{Gamma, ChiSquared, FisherF, StudentT}; #[cfg(feature="std")] use distributions::normal::{Normal, LogNormal}; - use distributions::range::{Range, SampleRange}; + use distributions::range::{Range, RangeImpl}; impl<'a, T: Clone> Sample for WeightedChoice<'a, T> { fn sample(&mut self, rng: &mut R) -> T { @@ -87,13 +85,13 @@ mod impls { } } - impl Sample for Range { - fn sample(&mut self, rng: &mut R) -> Sup { + impl Sample for Range { + fn sample(&mut self, rng: &mut R) -> T::X { Distribution::sample(self, rng) } } - impl IndependentSample for Range { - fn ind_sample(&self, rng: &mut R) -> Sup { + impl IndependentSample for Range { + fn ind_sample(&self, rng: &mut R) -> T::X { Distribution::sample(self, rng) } } @@ -135,8 +133,8 @@ impl<'a, T, D: Distribution> Distribution for &'a D { /// A generic random value distribution. Generates values for various types /// with numerically uniform distribution. /// -/// For floating-point numbers, this generates values from the half-open range -/// `[0, 1)` (excluding 1). See also [`Open01`] and [`Closed01`] for alternatives. +/// For floating-point numbers, this generates values from the open range +/// `(0, 1)` (i.e. excluding 0.0 and 1.0). /// /// ## Built-in Implementations /// @@ -148,13 +146,11 @@ impl<'a, T, D: Distribution> Distribution for &'a D { /// over all values of the type. /// * `char`: Uniformly distributed over all Unicode scalar values, i.e. all /// code points in the range `0...0x10_FFFF`, except for the range -/// `0xD800...0xDFFF` (the surrogate code points). This includes +/// `0xD800...0xDFFF` (the surrogate code points). This includes /// unassigned/reserved code points. /// * `bool`: Generates `false` or `true`, each with probability 0.5. /// * Floating point types (`f32` and `f64`): Uniformly distributed in the -/// half-open range `[0, 1)`. (The [`Open01`], [`Closed01`], [`Exp1`], and -/// [`StandardNormal`] distributions produce floating point numbers with -/// alternative ranges or distributions.) +/// open range `(0, 1)`. /// /// The following aggregate types also implement the distribution `Uniform` as /// long as their component types implement it: @@ -185,8 +181,6 @@ impl<'a, T, D: Distribution> Distribution for &'a D { /// println!("f32 from [0,1): {}", val); /// ``` /// -/// [`Open01`]: struct.Open01.html -/// [`Closed01`]: struct.Closed01.html /// [`Exp1`]: struct.Exp1.html /// [`StandardNormal`]: struct.StandardNormal.html #[derive(Debug)] @@ -236,7 +230,7 @@ pub struct Weighted { #[derive(Debug)] pub struct WeightedChoice<'a, T:'a> { items: &'a mut [Weighted], - weight_range: Range + weight_range: Range>, } impl<'a, T: Clone> WeightedChoice<'a, T> { @@ -344,21 +338,26 @@ fn ziggurat( mut pdf: P, mut zero_case: Z) -> f64 where P: FnMut(f64) -> f64, Z: FnMut(&mut R, f64) -> f64 { - const SCALE: f64 = (1u64 << 53) as f64; loop { - // reimplement the f64 generation as an optimisation suggested - // by the Doornik paper: we have a lot of precision-space - // (i.e. there are 11 bits of the 64 of a u64 to use after - // creating a f64), so we might as well reuse some to save - // generating a whole extra random number. (Seems to be 15% - // faster.) - let bits: u64 = rng.gen(); - let i = (bits & 0xff) as usize; - let f = (bits >> 11) as f64 / SCALE; - - // u is either U(-1, 1) or U(0, 1) depending on if this is a - // symmetric distribution or not. - let u = if symmetric {2.0 * f - 1.0} else {f}; + // As an optimisation we re-implement the conversion to a f64. + // From the remaining 12 most significant bits we use 8 to construct `i`. + // This saves us generating a whole extra random number, while the added + // precision of using 64 bits for f64 does not buy us much. + let bits = rng.next_u64(); + let i = bits as usize & 0xff; + + let u = if symmetric { + // Convert to a value in the range [2,4) and substract to get [-1,1) + // We can't convert to an open range directly, that would require + // substracting `3.0 - EPSILON`, which is not representable. + // It is possible with an extra step, but an open range does not + // seem neccesary for the ziggurat algorithm anyway. + (bits >> 12).into_float_with_exponent(1) - 3.0 + } else { + // Convert to a value in the range [1,2) and substract to get (0,1) + (bits >> 12).into_float_with_exponent(0) + - (1.0 - ::core::f64::EPSILON / 2.0) + }; let x = u * x_tab[i]; let test_x = if symmetric { x.abs() } else {x}; @@ -386,17 +385,22 @@ mod tests { #[test] fn test_weighted_choice() { // this makes assumptions about the internal implementation of - // WeightedChoice, specifically: it doesn't reorder the items, - // it doesn't do weird things to the RNG (so 0 maps to 0, 1 to - // 1, internally; modulo a modulo operation). + // WeightedChoice. It may fail when the implementation in + // `distributions::range::RangeInt changes. macro_rules! t { ($items:expr, $expected:expr) => {{ let mut items = $items; + let mut total_weight = 0; + for item in &items { total_weight += item.weight; } + let wc = WeightedChoice::new(&mut items); let expected = $expected; - let mut rng = StepRng::new(0, 1); + // Use extremely large steps between the random numbers, because + // we test with small ranges and RangeInt is designed to prefer + // the most significant bits. + let mut rng = StepRng::new(0, !0 / (total_weight as u64)); for &val in expected.iter() { assert_eq!(wc.sample(&mut rng), val) @@ -411,12 +415,12 @@ mod tests { Weighted { weight: 2, item: 21}, Weighted { weight: 0, item: 22}, Weighted { weight: 1, item: 23}], - [21,21, 23]); + [21, 21, 23]); // different weights t!([Weighted { weight: 4, item: 30}, Weighted { weight: 3, item: 31}], - [30,30,30,30, 31,31,31]); + [30, 31, 30, 31, 30, 31, 30]); // check that we're binary searching // correctly with some vectors of odd @@ -434,7 +438,7 @@ mod tests { Weighted { weight: 1, item: 54}, Weighted { weight: 1, item: 55}, Weighted { weight: 1, item: 56}], - [50, 51, 52, 53, 54, 55, 56]); + [50, 54, 51, 55, 52, 56, 53]); } #[test] diff --git a/src/distributions/normal.rs b/src/distributions/normal.rs index 089cb9046d5..95d7b29f587 100644 --- a/src/distributions/normal.rs +++ b/src/distributions/normal.rs @@ -11,7 +11,7 @@ //! The normal and derived distributions. use {Rng}; -use distributions::{ziggurat, ziggurat_tables, Distribution, Open01}; +use distributions::{ziggurat, ziggurat_tables, Distribution, Uniform}; /// Samples floating-point numbers according to the normal distribution /// `N(0, 1)` (a.k.a. a standard normal, or Gaussian). This is equivalent to @@ -55,8 +55,8 @@ impl Distribution for StandardNormal { let mut y = 0.0f64; while -2.0 * y < x * x { - let x_: f64 = rng.sample(Open01); - let y_: f64 = rng.sample(Open01); + let x_: f64 = rng.sample(Uniform); + let y_: f64 = rng.sample(Uniform); x = x_.ln() / ziggurat_tables::ZIG_NORM_R; y = y_.ln(); diff --git a/src/distributions/other.rs b/src/distributions/other.rs index 873043addac..2c433a6856d 100644 --- a/src/distributions/other.rs +++ b/src/distributions/other.rs @@ -13,19 +13,17 @@ use core::char; use {Rng}; -use distributions::{Distribution, Uniform}; +use distributions::{Distribution, Uniform, Range}; impl Distribution for Uniform { #[inline] fn sample(&self, rng: &mut R) -> char { - // a char is 21 bits - const CHAR_MASK: u32 = 0x001f_ffff; + let range = Range::new(0u32, 0x11_0000); loop { - // Rejection sampling. About 0.2% of numbers with at most - // 21-bits are invalid codepoints (surrogates), so this - // will succeed first go almost every time. - match char::from_u32(rng.next_u32() & CHAR_MASK) { + match char::from_u32(range.sample(rng)) { Some(c) => return c, + // About 0.2% of numbers in the range 0..0x110000 are invalid + // codepoints (surrogates). None => {} } } @@ -35,7 +33,11 @@ impl Distribution for Uniform { impl Distribution for Uniform { #[inline] fn sample(&self, rng: &mut R) -> bool { - rng.gen::() & 1 == 1 + // We can compare against an arbitrary bit of an u32 to get a bool. + // Because the least significant bits of a lower quality RNG can have + // simple patterns, we compare against the most significant bit. This is + // easiest done using a sign test. + (rng.next_u32() as i32) < 0 } } diff --git a/src/distributions/range.rs b/src/distributions/range.rs index b5cbae4c50e..5349fc05506 100644 --- a/src/distributions/range.rs +++ b/src/distributions/range.rs @@ -1,4 +1,4 @@ -// Copyright 2013 The Rust Project Developers. See the COPYRIGHT +// Copyright 2017 The Rust Project Developers. See the COPYRIGHT // file at the top-level directory of this distribution and at // https://rust-lang.org/COPYRIGHT. // @@ -8,28 +8,33 @@ // option. This file may not be copied, modified, or distributed // except according to those terms. -//! Generating numbers between two others. - -// this is surprisingly complicated to be both generic & correct - -use core::num::Wrapping as w; +//! A distribution generating numbers within a given range. use Rng; -use distributions::Distribution; +use distributions::{Distribution, Uniform}; +use distributions::float::IntoFloat; /// Sample values uniformly between two bounds. /// -/// This gives a uniform distribution (assuming the RNG used to sample -/// it is itself uniform & the `SampleRange` implementation for the -/// given type is correct), even for edge cases like `low = 0u8`, -/// `high = 170u8`, for which a naive modulo operation would return -/// numbers less than 85 with double the probability to those greater -/// than 85. +/// `Range::new` and `Range::new_inclusive` will set up a `Range`, which does +/// some preparations up front to make sampling values faster. +/// `Range::sample_single` is optimized for sampling values once or only a +/// limited number of times from a range. +/// +/// If you need to sample many values from a range, consider using `new` or +/// `new_inclusive`. This is also the best choice if the range is constant, +/// because then the preparations can be evaluated at compile-time. +/// Otherwise `sample_single` may be the best choice. +/// +/// Sampling uniformly from a range can be surprisingly complicated to be both +/// generic and correct. Consider for example edge cases like `low = 0u8`, +/// `high = 170u8`, for which a naive modulo operation would return numbers less +/// than 85 with double the probability to those greater than 85. /// -/// Types should attempt to sample in `[low, high)`, i.e., not -/// including `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. +/// Types should attempt to sample in `[low, high)` for `Range::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. /// /// # Example /// @@ -47,80 +52,266 @@ use distributions::Distribution; /// } /// ``` #[derive(Clone, Copy, Debug)] -pub struct Range { - low: X, - range: X, - accept_zone: X +pub struct Range { + inner: T, } -impl Range { - /// Create a new `Range` instance that samples uniformly from - /// `[low, high)`. Panics if `low >= high`. - pub fn new(low: X, high: X) -> Range { +/// Ignore the `RangeInt` parameterisation; these non-member functions are +/// available to all range types. (Rust requires a type, and won't accept +/// `T: RangeImpl`. Consider this a minor language issue.) +impl Range> { + /// Create a new `Range` instance which samples uniformly from the half + /// open range `[low, high)` (excluding `high`). Panics if `low >= high`. + pub fn new(low: X, high: X) -> Range { + assert!(low < high, "Range::new called with `low >= high`"); + Range { inner: RangeImpl::new(low, high) } + } + + /// Create a new `Range` instance which samples uniformly from the closed + /// range `[low, high]` (inclusive). Panics if `low >= high`. + pub fn new_inclusive(low: X, high: X) -> Range { assert!(low < high, "Range::new called with `low >= high`"); - SampleRange::construct_range(low, high) + Range { inner: RangeImpl::new_inclusive(low, high) } + } + + /// Sample a single value uniformly from `[low, high)`. + /// Panics if `low >= high`. + pub fn sample_single(low: X, high: X, rng: &mut R) -> X { + assert!(low < high, "Range::sample_single called with low >= high"); + X::T::sample_single(low, high, rng) } } -impl Distribution for Range { - fn sample(&self, rng: &mut R) -> Sup { - SampleRange::sample_range(self, rng) +impl Distribution for Range { + fn sample(&self, rng: &mut R) -> T::X { + self.inner.sample(rng) } } -/// The helper trait for types that have a sensible way to sample -/// uniformly between two values. This should not be used directly, -/// and is only to facilitate `Range`. -pub trait SampleRange : Sized { - /// Construct the `Range` object that `sample_range` - /// requires. This should not ever be called directly, only via - /// `Range::new`, which will check that `low < high`, so this - /// function doesn't have to repeat the check. - fn construct_range(low: Self, high: Self) -> Range; - - /// Sample a value from the given `Range` with the given `Rng` as - /// a source of randomness. - fn sample_range(r: &Range, rng: &mut R) -> Self; +/// Helper trait for creating objects using the correct implementation of +/// `RangeImpl` for the sampling type; this enables `Range::new(a, b)` to work. +pub trait SampleRange: PartialOrd+Sized { + type T: RangeImpl; } -macro_rules! integer_impl { - ($ty:ty, $unsigned:ident) => { +/// Helper trait handling actual range sampling. +/// +/// If you want to implement `Range` sampling for your own type, then +/// implement both this trait and `SampleRange`: +/// +/// ```rust +/// use rand::{Rng, thread_rng}; +/// use rand::distributions::Distribution; +/// use rand::distributions::range::{Range, SampleRange, RangeImpl, RangeFloat}; +/// +/// #[derive(Clone, Copy, PartialEq, PartialOrd)] +/// struct MyF32(f32); +/// +/// #[derive(Clone, Copy, Debug)] +/// struct RangeMyF32 { +/// inner: RangeFloat, +/// } +/// impl RangeImpl for RangeMyF32 { +/// type X = MyF32; +/// fn new(low: Self::X, high: Self::X) -> Self { +/// RangeMyF32 { +/// inner: RangeFloat::::new(low.0, high.0), +/// } +/// } +/// fn new_inclusive(low: Self::X, high: Self::X) -> Self { +/// RangeImpl::new(low, high) +/// } +/// fn sample(&self, rng: &mut R) -> Self::X { +/// MyF32(self.inner.sample(rng)) +/// } +/// } +/// +/// impl SampleRange for MyF32 { +/// type T = RangeMyF32; +/// } +/// +/// let (low, high) = (MyF32(17.0f32), MyF32(22.0f32)); +/// let range = Range::new(low, high); +/// let x = range.sample(&mut thread_rng()); +/// ``` +pub trait RangeImpl: Sized { + /// The type sampled by this implementation. + type X: PartialOrd; + + /// Construct self, with inclusive lower bound and exclusive upper bound + /// `[low, high)`. + /// + /// Usually users should not call this directly but instead use + /// `Range::new`, which asserts that `low < high` before calling this. + fn new(low: Self::X, high: Self::X) -> Self; + + /// Construct self, with inclusive bounds `[low, high]`. + /// + /// Usually users should not call this directly but instead use + /// `Range::new_inclusive`, which asserts that `low < high` before calling + /// this. + fn new_inclusive(low: Self::X, high: Self::X) -> Self; + + /// Sample a value. + fn sample(&self, rng: &mut R) -> Self::X; + + /// Sample a single value uniformly from a range with inclusive lower bound + /// and exclusive upper bound `[low, high)`. + /// + /// Usually users should not call this directly but instead use + /// `Range::sample_single`, which asserts that `low < high` before calling + /// this. + /// + /// Via this method range implementations can provide a method optimized for + /// sampling only a limited number of values from range. The default + /// implementation just sets up a range with `RangeImpl::new` and samples + /// from that. + fn sample_single(low: Self::X, high: Self::X, rng: &mut R) + -> Self::X + { + let range: Self = RangeImpl::new(low, high); + range.sample(rng) + } +} + +/// Implementation of `RangeImpl` for integer types. +#[derive(Clone, Copy, Debug)] +pub struct RangeInt { + low: X, + range: X, + zone: X, +} + +macro_rules! range_int_impl { + ($ty:ty, $signed:ty, $unsigned:ident, + $i_large:ident, $u_large:ident) => { impl SampleRange for $ty { - // we play free and fast with unsigned vs signed here + type T = RangeInt<$ty>; + } + + impl RangeImpl for RangeInt<$ty> { + // We play free and fast with unsigned vs signed here // (when $ty is signed), but that's fine, since the // contract of this macro is for $ty and $unsigned to be - // "bit-equal", so casting between them is a no-op & a - // bijection. + // "bit-equal", so casting between them is a no-op. + + type X = $ty; - #[inline] - fn construct_range(low: $ty, high: $ty) -> Range<$ty> { - let range = (w(high as $unsigned) - w(low as $unsigned)).0; - let unsigned_max: $unsigned = ::core::$unsigned::MAX; + #[inline] // if the range is constant, this helps LLVM to do the + // calculations at compile-time. + fn new(low: Self::X, high: Self::X) -> Self { + RangeImpl::new_inclusive(low, high - 1) + } + + #[inline] // if the range is constant, this helps LLVM to do the + // calculations at compile-time. + fn new_inclusive(low: Self::X, high: Self::X) -> Self { + // For a closed range the number of possible numbers we should + // generate is `range = (high - low + 1)`. It is not possible to + // end up with a uniform distribution if we map _all_ the random + // integers that can be generated to this range. We have to map + // integers from a `zone` that is a multiple of the range. The + // rest of the integers, that cause a bias, are rejected. + // + // The problem with `range` is that to cover the full range of + // the type, it has to store `unsigned_max + 1`, which can't be + // represented. But if the range covers the full range of the + // type, no modulus is needed. A range of size 0 can't exist, so + // we use that to represent this special case. Wrapping + // arithmetic even makes representing `unsigned_max + 1` as 0 + // simple. + // + // We don't calculate zone directly, but first calculate the + // number of integers to reject. To handle `unsigned_max + 1` + // not fitting in the type, we use: + // ints_to_reject = (unsigned_max + 1) % range; + // ints_to_reject = (unsigned_max - range + 1) % range; + // + // The smallest integer prngs generate is u32. That is why for + // small integer sizes (i8/u8 and i16/u16) there is an + // optimisation: don't pick the largest zone that can fit in the + // small type, but pick the largest zone that can fit in an u32. + // This improves the chance to get a random integer that fits in + // the zone to 998 in 1000 in the worst case. + // + // There is a problem however: we can't store such a large range + // in `RangeInt`, that can only hold values of the size of $ty. + // `ints_to_reject` is always less than half the size of the + // small integer. For an u8 it only ever uses 7 bits. This means + // that all but the last 7 bits of `zone` are always 1's (or 15 + // in the case of u16). So nothing is lost by trucating `zone`. + // + // An alternative to using a modulus is widening multiply: + // After a widening multiply by `range`, the result is in the + // high word. Then comparing the low word against `zone` makes + // sure our distribution is uniform. + let unsigned_max: $u_large = ::core::$u_large::MAX; - // this is the largest number that fits into $unsigned - // that `range` divides evenly, so, if we've sampled - // `n` uniformly from this region, then `n % range` is - // uniform in [0, range) - let zone = unsigned_max - unsigned_max % range; + let range = (high as $u_large) + .wrapping_sub(low as $u_large) + .wrapping_add(1); + let ints_to_reject = + if range > 0 { + (unsigned_max - range + 1) % range + } else { + 0 + }; + let zone = unsigned_max - ints_to_reject; - Range { + RangeInt { low: low, + // These are really $unsigned values, but store as $ty: range: range as $ty, - accept_zone: zone as $ty + zone: zone as $ty + } + } + + fn sample(&self, rng: &mut R) -> Self::X { + let range = self.range as $unsigned as $u_large; + if range > 0 { + // Some casting to recover the trucated bits of `zone`: + // First bit-cast to a signed int. Next sign-extend to the + // larger type. Then bit-cast to unsigned. + // For types that already have the right size, all the + // casting is a no-op. + let zone = self.zone as $signed as $i_large as $u_large; + loop { + let v: $u_large = Uniform.sample(rng); + let (hi, lo) = v.wmul(range); + if lo <= zone { + return self.low.wrapping_add(hi as $ty); + } + } + } else { + // Sample from the entire integer range. + Uniform.sample(rng) } } - #[inline] - fn sample_range(r: &Range<$ty>, rng: &mut R) -> $ty { + fn sample_single(low: Self::X, + high: Self::X, + rng: &mut R) -> Self::X + { + let range = (high as $u_large) + .wrapping_sub(low as $u_large); + let zone = + if ::core::$unsigned::MAX <= ::core::u16::MAX as $unsigned { + // Using a modulus is faster than the approximation for + // i8 and i16. I suppose we trade the cost of one + // modulus for near-perfect branch prediction. + let unsigned_max: $u_large = ::core::$u_large::MAX; + let ints_to_reject = (unsigned_max - range + 1) % range; + unsigned_max - ints_to_reject + } else { + // conservative but fast approximation + range << range.leading_zeros() + }; + loop { - // rejection sample - let v = rng.gen::<$unsigned>(); - // until we find something that fits into the - // region which r.range evenly divides (this will - // be uniformly distributed) - if v < r.accept_zone as $unsigned { - // and return it, with some adjustments - return (w(r.low) + w((v % r.range as $unsigned) as $ty)).0; + let v: $u_large = Uniform.sample(rng); + let (hi, lo) = v.wmul(range); + if lo <= zone { + return low.wrapping_add(hi as $ty); } } } @@ -128,51 +319,166 @@ macro_rules! integer_impl { } } -integer_impl! { i8, u8 } -integer_impl! { i16, u16 } -integer_impl! { i32, u32 } -integer_impl! { i64, u64 } +range_int_impl! { i8, i8, u8, i32, u32 } +range_int_impl! { i16, i16, u16, i32, u32 } +range_int_impl! { i32, i32, u32, i32, u32 } +range_int_impl! { i64, i64, u64, i64, u64 } +#[cfg(feature = "i128_support")] +range_int_impl! { i128, i128, u128, u128, u128 } +range_int_impl! { isize, isize, usize, isize, usize } +range_int_impl! { u8, i8, u8, i32, u32 } +range_int_impl! { u16, i16, u16, i32, u32 } +range_int_impl! { u32, i32, u32, i32, u32 } +range_int_impl! { u64, i64, u64, i64, u64 } +range_int_impl! { usize, isize, usize, isize, usize } +#[cfg(feature = "i128_support")] +range_int_impl! { u128, u128, u128, i128, u128 } + + +trait WideningMultiply { + type Output; + + fn wmul(self, x: RHS) -> Self::Output; +} + +macro_rules! wmul_impl { + ($ty:ty, $wide:ty, $shift:expr) => { + impl WideningMultiply for $ty { + type Output = ($ty, $ty); + + #[inline(always)] + fn wmul(self, x: $ty) -> Self::Output { + let tmp = (self as $wide) * (x as $wide); + ((tmp >> $shift) as $ty, tmp as $ty) + } + } + } +} + +wmul_impl! { u8, u16, 8 } +wmul_impl! { u16, u32, 16 } +wmul_impl! { u32, u64, 32 } #[cfg(feature = "i128_support")] -integer_impl! { i128, u128 } -integer_impl! { isize, usize } -integer_impl! { u8, u8 } -integer_impl! { u16, u16 } -integer_impl! { u32, u32 } -integer_impl! { u64, u64 } +wmul_impl! { u64, u128, 64 } + +// This code is a translation of the __mulddi3 function in LLVM's +// compiler-rt. It is an optimised variant of the common method +// `(a + b) * (c + d) = ac + ad + bc + bd`. +// +// For some reason LLVM can optimise the C version very well, but +// keeps shuffeling registers in this Rust translation. +macro_rules! wmul_impl_large { + ($ty:ty, $half:expr) => { + impl WideningMultiply for $ty { + type Output = ($ty, $ty); + + #[inline(always)] + fn wmul(self, b: $ty) -> Self::Output { + const LOWER_MASK: $ty = !0 >> $half; + let mut low = (self & LOWER_MASK).wrapping_mul(b & LOWER_MASK); + let mut t = low >> $half; + low &= LOWER_MASK; + t += (self >> $half).wrapping_mul(b & LOWER_MASK); + low += (t & LOWER_MASK) << $half; + let mut high = t >> $half; + t = low >> $half; + low &= LOWER_MASK; + t += (b >> $half).wrapping_mul(self & LOWER_MASK); + low += (t & LOWER_MASK) << $half; + high += t >> $half; + high += (self >> $half).wrapping_mul(b >> $half); + + (high, low) + } + } + } +} + +#[cfg(not(feature = "i128_support"))] +wmul_impl_large! { u64, 32 } #[cfg(feature = "i128_support")] -integer_impl! { u128, u128 } -integer_impl! { usize, usize } +wmul_impl_large! { u128, 64 } + -macro_rules! float_impl { +macro_rules! wmul_impl_usize { ($ty:ty) => { + impl WideningMultiply for usize { + type Output = (usize, usize); + + #[inline(always)] + fn wmul(self, x: usize) -> Self::Output { + let (high, low) = (self as $ty).wmul(x as $ty); + (high as usize, low as usize) + } + } + } +} + +#[cfg(target_pointer_width = "32")] +wmul_impl_usize! { u32 } +#[cfg(target_pointer_width = "64")] +wmul_impl_usize! { u64 } + + + +/// Implementation of `RangeImpl` for float types. +#[derive(Clone, Copy, Debug)] +pub struct RangeFloat { + scale: X, + offset: X, +} + +macro_rules! range_float_impl { + ($ty:ty, $bits_to_discard:expr, $next_u:ident) => { impl SampleRange for $ty { - fn construct_range(low: $ty, high: $ty) -> Range<$ty> { - Range { - low: low, - range: high - low, - accept_zone: 0.0 // unused + type T = RangeFloat<$ty>; + } + + impl RangeImpl for RangeFloat<$ty> { + type X = $ty; + + fn new(low: Self::X, high: Self::X) -> Self { + let scale = high - low; + let offset = low - scale; + RangeFloat { + scale: scale, + offset: offset, } } - fn sample_range(r: &Range<$ty>, rng: &mut R) -> $ty { - r.low + r.range * rng.gen::<$ty>() + + fn new_inclusive(low: Self::X, high: Self::X) -> Self { + // Same as `new`, because the boundaries of a floats range are + // (at least currently) not exact due to rounding errors. + RangeImpl::new(low, high) + } + + 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); + // Doing multiply before addition allows some architectures to + // use a single instruction. + value1_2 * self.scale + self.offset } } } } -float_impl! { f32 } -float_impl! { f64 } +range_float_impl! { f32, 32 - 23, next_u32 } +range_float_impl! { f64, 64 - 52, next_u64 } + #[cfg(test)] mod tests { - use distributions::Distribution; - use super::Range as Range; + use Rng; + use distributions::range::{Range, RangeImpl, RangeFloat, SampleRange}; #[should_panic] #[test] fn test_range_bad_limits_equal() { Range::new(10, 10); } + #[should_panic] #[test] fn test_range_bad_limits_flipped() { @@ -189,23 +495,19 @@ mod tests { (10, 127), (::core::$ty::MIN, ::core::$ty::MAX)]; for &(low, high) in v.iter() { - let sampler: Range<$ty> = Range::new(low, high); + let my_range = Range::new(low, high); for _ in 0..1000 { - let v = sampler.sample(&mut rng); - assert!(low <= v && v < high); - let v = sampler.sample(&mut rng); + let v: $ty = rng.sample(my_range); assert!(low <= v && v < high); } } )* }} } - #[cfg(not(feature = "i128_support"))] t!(i8, i16, i32, i64, isize, u8, u16, u32, u64, usize); #[cfg(feature = "i128_support")] - t!(i8, i16, i32, i64, i128, isize, - u8, u16, u32, u64, u128, usize); + t!(i128, u128) } #[test] @@ -219,11 +521,9 @@ mod tests { (1e-35, 1e-25), (-1e35, 1e35)]; for &(low, high) in v.iter() { - let sampler: Range<$ty> = Range::new(low, high); + let my_range = Range::new(low, high); for _ in 0..1000 { - let v = sampler.sample(&mut rng); - assert!(low <= v && v < high); - let v = sampler.sample(&mut rng); + let v: $ty = rng.sample(my_range); assert!(low <= v && v < high); } } @@ -233,4 +533,40 @@ mod tests { t!(f32, f64) } + #[test] + fn test_custom_range() { + #[derive(Clone, Copy, PartialEq, PartialOrd)] + struct MyF32 { + x: f32, + } + #[derive(Clone, Copy, Debug)] + struct RangeMyF32 { + inner: RangeFloat, + } + impl RangeImpl for RangeMyF32 { + type X = MyF32; + fn new(low: Self::X, high: Self::X) -> Self { + RangeMyF32 { + inner: RangeFloat::::new(low.x, high.x), + } + } + fn new_inclusive(low: Self::X, high: Self::X) -> Self { + RangeImpl::new(low, high) + } + fn sample(&self, rng: &mut R) -> Self::X { + MyF32 { x: self.inner.sample(rng) } + } + } + impl SampleRange for MyF32 { + type T = RangeMyF32; + } + + let (low, high) = (MyF32{ x: 17.0f32 }, MyF32{ x: 22.0f32 }); + let range = Range::new(low, high); + let mut rng = ::test::rng(804); + for _ in 0..100 { + let x: MyF32 = rng.sample(range); + assert!(low <= x && x < high); + } + } } diff --git a/src/lib.rs b/src/lib.rs index ad7612118e1..c6d546d604c 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -155,6 +155,7 @@ //! ``` //! use rand::Rng; //! use rand::distributions::{Distribution, Range}; +//! use rand::distributions::range::RangeInt; //! //! struct SimulationResult { //! win: bool, @@ -162,7 +163,7 @@ //! } //! //! // Run a single simulation of the Monty Hall problem. -//! fn simulate(random_door: &Range, rng: &mut R) +//! fn simulate(random_door: &Range>, rng: &mut R) //! -> SimulationResult { //! let car = random_door.sample(rng); //! @@ -203,7 +204,7 @@ //! let num_simulations = 10000; //! //! let mut rng = rand::thread_rng(); -//! let random_door = Range::new(0, 3); +//! let random_door = Range::new(0u32, 3); //! //! let (mut switch_wins, mut switch_losses) = (0, 0); //! let (mut keep_wins, mut keep_losses) = (0, 0); @@ -630,8 +631,7 @@ pub trait Rng: RngCore + Sized { /// println!("{}", m); /// ``` fn gen_range(&mut self, low: T, high: T) -> T { - assert!(low < high, "Rng::gen_range called with low >= high"); - Range::new(low, high).sample(self) + Range::sample_single(low, high, self) } /// Return a bool with a 1 in n chance of true