diff --git a/benches/distributions.rs b/benches/distributions.rs index 6d15aeb3f26..de13e9d3a7a 100644 --- a/benches/distributions.rs +++ b/benches/distributions.rs @@ -95,6 +95,8 @@ distr!(distr_uniform_codepoint, char, Standard); distr_float!(distr_uniform_f32, f32, Standard); distr_float!(distr_uniform_f64, f64, Standard); +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)); diff --git a/src/distributions/float.rs b/src/distributions/float.rs index b1b76852fcd..d383fbdbf67 100644 --- a/src/distributions/float.rs +++ b/src/distributions/float.rs @@ -10,10 +10,37 @@ //! Basic floating-point number distributions -use core::mem; +use core::{cmp, mem}; use Rng; use distributions::{Distribution, Standard}; +/// Generate a floating point number in the half-open interval `[0, 1)` with a +/// uniform distribution, with as much precision as the floating-point type +/// can represent, including sub-normals. +/// +/// Technically 0 is representable, but the probability of occurrence is +/// remote (1 in 2^149 for `f32` or 1 in 2^1074 for `f64`). +/// +/// This is different from `Uniform` in that it uses as many random bits as +/// required to get high precision close to 0. Normally only a single call to +/// the source RNG is required (32 bits for `f32` or 64 bits for `f64`); 1 in +/// 2^9 (`f32`) or 2^12 (`f64`) samples need an extra call; of these 1 in 2^32 +/// or 1 in 2^64 require a third call, etc.; i.e. even for `f32` a third call is +/// almost impossible to observe with an unbiased RNG. Due to the extra logic +/// there is some performance overhead relative to `Uniform`; this is more +/// significant for `f32` than for `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); +/// ``` +#[derive(Clone, Copy, Debug)] +pub struct HighPrecision01; + pub(crate) trait IntoFloat { type F; @@ -54,6 +81,66 @@ macro_rules! float_impls { fraction.into_float_with_exponent(0) - (1.0 - EPSILON / 2.0) } } + + impl Distribution<$ty> for HighPrecision01 { + /// Generate a floating point number in the half-open interval + /// `[0, 1)` with a uniform distribution. See [`HighPrecision01`]. + /// + /// # 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 + /// zeros in the random bits. In the case where all bits are zero, + /// we simply generate a new random number and add the number of + /// trailing zeros to the previous count (up to maximum exponent). + /// + /// Each exponent is responsible for a piece of the distribution + /// between [0,1). We take the above exponent, add 1 and negate; + /// thus with probability 1/2 we have exponent -1 which fills the + /// range [0.5,1); with probability 1/4 we have exponent -2 which + /// fills the range [0.25,0.5), etc. If the exponent reaches the + /// minimum allowed, the floating-point format drops the implied + /// fraction bit, thus allowing numbers down to 0 to be sampled. + /// + /// [`HighPrecision01`]: struct.HighPrecision01.html + #[inline] + fn sample(&self, rng: &mut R) -> $ty { + // Unusual case. Separate function to allow inlining of rest. + #[inline(never)] + fn fallback(mut exp: i32, fraction: $uty, rng: &mut R) -> $ty { + // Performance impact of code here is negligible. + let bits = rng.gen::<$uty>(); + exp += bits.trailing_zeros() as i32; + // If RNG were guaranteed unbiased we could skip the + // check against exp; unfortunately it may be. + // Worst case ("zeros" RNG) has recursion depth 16. + if bits == 0 && exp < $exponent_bias { + return fallback(exp, fraction, rng); + } + exp = cmp::min(exp, $exponent_bias); + fraction.into_float_with_exponent(-exp) + } + + let fraction_mask = (1 << $fraction_bits) - 1; + let value = rng.$next_u(); + + let fraction = value & fraction_mask; + let remaining = value >> $fraction_bits; + if remaining == 0 { + // exp is compile-time constant so this reduces to a function call: + let size_bits = (mem::size_of::<$ty>() * 8) as i32; + let exp = (size_bits - $fraction_bits as i32) + 1; + return fallback(exp, fraction, rng); + } + + // Usual case: exponent from -1 to -9 (f32) or -12 (f64) + let exp = remaining.trailing_zeros() as i32 + 1; + fraction.into_float_with_exponent(-exp) + } + } } } float_impls! { f32, u32, 23, 127, next_u32 } @@ -62,7 +149,8 @@ float_impls! { f64, u64, 52, 1023, next_u64 } #[cfg(test)] mod tests { - use Rng; + use {Rng}; + use distributions::HighPrecision01; use mock::StepRng; const EPSILON32: f32 = ::core::f32::EPSILON; @@ -86,4 +174,45 @@ mod tests { assert_eq!(max.gen::(), 1.0 - EPSILON32 / 2.0); assert_eq!(max.gen::(), 1.0 - EPSILON64 / 2.0); } + + #[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::(HighPrecision01), 0.0); + assert_eq!(zeros.sample::(HighPrecision01), 0.0); + + let mut ones = StepRng::new(0xffff_ffff_ffff_ffff, 0); + assert_eq!(ones.sample::(HighPrecision01), 0.99999994); + assert_eq!(ones.sample::(HighPrecision01), 0.9999999999999999); + } + + #[cfg(feature="std")] mod mean { + use {Rng, SmallRng, SeedableRng, thread_rng}; + use distributions::{Uniform, HighPrecision01}; + + macro_rules! test_mean { + ($name:ident, $ty:ty, $distr:expr) => { + #[test] + fn $name() { + // TODO: no need to &mut here: + let mut rng = SmallRng::from_rng(&mut thread_rng()).unwrap(); + let mut total: $ty = 0.0; + const N: u32 = 1_000_000; + for _ in 0..N { + total += rng.sample::<$ty, _>($distr); + } + let avg = total / (N as $ty); + //println!("average over {} samples: {}", N, avg); + assert!(0.499 < avg && avg < 0.501); + } + } } + + test_mean!(test_mean_f32, f32, Uniform); + test_mean!(test_mean_f64, f64, Uniform); + test_mean!(test_mean_high_f32, f32, HighPrecision01); + test_mean!(test_mean_high_f64, f64, HighPrecision01); + } } diff --git a/src/distributions/mod.rs b/src/distributions/mod.rs index 6afbf923709..3107692301e 100644 --- a/src/distributions/mod.rs +++ b/src/distributions/mod.rs @@ -18,6 +18,7 @@ use Rng; pub use self::other::Alphanumeric; +pub use self::float::HighPrecision01; pub use self::range::Range; #[cfg(feature="std")] pub use self::gamma::{Gamma, ChiSquared, FisherF, StudentT};