Higher precision for floats in range [0,1)

benches/distributions.rs

@@ -23,10 +23,10 @@ fn distr_baseline(b: &mut Bencher) {
 
     b.iter(|| {
         for _ in 0..::RAND_BENCH_N {
-            f64::rand(&mut rng, Default);
+            u64::rand(&mut rng, Default);
         }
     });
-    b.bytes = size_of::<f64>() as u64 * ::RAND_BENCH_N;
+    b.bytes = size_of::<u64>() as u64 * ::RAND_BENCH_N;
 }
 
 

benches/misc.rs

@@ -10,7 +10,7 @@ use test::{black_box, Bencher};
 use rand::StdRng;
 use rand::prng::XorShiftRng;
 use rand::sequences::{sample, Shuffle};
-use rand::distributions::{Rand, Uniform, Uniform01};
+use rand::distributions::{Rand, Uniform, Uniform01, Closed01, Open01};
 
 #[bench]
 fn misc_baseline_32(b: &mut Bencher) {
@@ -35,7 +35,7 @@ fn misc_baseline_64(b: &mut Bencher) {
 }
 
 #[bench]
-fn misc_convert_f32(b: &mut Bencher) {
+fn misc_uniform01_f32(b: &mut Bencher) {
     let mut rng = StdRng::new().unwrap();
     b.iter(|| {
         for _ in 0..RAND_BENCH_N {
@@ -46,7 +46,7 @@ fn misc_convert_f32(b: &mut Bencher) {
 }
 
 #[bench]
-fn misc_convert_f64(b: &mut Bencher) {
+fn misc_uniform01_f64(b: &mut Bencher) {
     let mut rng = StdRng::new().unwrap();
     b.iter(|| {
         for _ in 0..RAND_BENCH_N {
@@ -56,6 +56,50 @@ fn misc_convert_f64(b: &mut Bencher) {
     b.bytes = size_of::<f64>() as u64 * RAND_BENCH_N;
 }
 
+#[bench]
+fn misc_closed01_f32(b: &mut Bencher) {
+    let mut rng = StdRng::new().unwrap();
+    b.iter(|| {
+        for _ in 0..RAND_BENCH_N {
+            black_box(f32::rand(&mut rng, Closed01));
+        }
+    });
+    b.bytes = size_of::<f32>() as u64 * RAND_BENCH_N;
+}
+
+#[bench]
+fn misc_closed01_f64(b: &mut Bencher) {
+    let mut rng = StdRng::new().unwrap();
+    b.iter(|| {
+        for _ in 0..RAND_BENCH_N {
+            black_box(f64::rand(&mut rng, Closed01));
+        }
+    });
+    b.bytes = size_of::<f64>() as u64 * RAND_BENCH_N;
+}
+
+#[bench]
+fn misc_open01_f32(b: &mut Bencher) {
+    let mut rng = StdRng::new().unwrap();
+    b.iter(|| {
+        for _ in 0..RAND_BENCH_N {
+            black_box(f32::rand(&mut rng, Open01));
+        }
+    });
+    b.bytes = size_of::<f32>() as u64 * RAND_BENCH_N;
+}
+
+#[bench]
+fn misc_open01_f64(b: &mut Bencher) {
+    let mut rng = StdRng::new().unwrap();
+    b.iter(|| {
+        for _ in 0..RAND_BENCH_N {
+            black_box(f64::rand(&mut rng, Open01));
+        }
+    });
+    b.bytes = size_of::<f64>() as u64 * RAND_BENCH_N;
+}
+
 #[bench]
 fn misc_shuffle_100(b: &mut Bencher) {
     let mut rng = XorShiftRng::new().unwrap();

src/distributions/mod.rs

@@ -23,6 +23,7 @@ pub use self::default::Default;
 pub use self::uniform::{uniform, codepoint, ascii_word_char};
 pub use self::uniform::{Uniform, Uniform01, Open01, Closed01, AsciiWordChar};
 pub use self::range::{Range};
+use utils::FloatConversions;
 
 #[cfg(feature="std")]
 pub use self::gamma::{Gamma, ChiSquared, FisherF, StudentT};
@@ -158,28 +159,25 @@ fn ziggurat<R: Rng+?Sized, P, Z>(
             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.)
-        //
-        // This unfortunately misses out on the benefits of direct
-        // floating point generation if an RNG like dSMFT is
-        // used. (That is, such RNGs create floats directly, highly
-        // efficiently and overload next_f32/f64, so by not calling it
-        // this may be slower than it would be otherwise.)
-        // FIXME: investigate/optimise for the above.
+        // As an optimisation convert the random u64 to a f64 using only
+        // 53 bits, as many as will fit in the float's fraction.
+        // Of the remaining 11 least 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: u64 = uniform(rng);
-        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};
+        // FIXME: the distribution is not open, but closed-open.
+        //        Can that cause problems or a bias?
+        let u = if symmetric {
+                    bits.closed_open11_fixed()
+                } else {
+                    bits.closed_open01_fixed()
+                };
+        let i = (bits & 0xff) as usize;
+
         let x = u * x_tab[i];
 
         let test_x = if symmetric { x.abs() } else {x};

src/distributions/range2.rs

@@ -24,7 +24,8 @@
 use core::num::Wrapping as w;
 
 use Rng;
-use distributions::{Distribution, Uniform01, Rand};
+use distributions::Distribution;
+use utils::FloatConversions;
 
 /// Sample values uniformly between two bounds.
 ///
@@ -177,37 +178,116 @@ range_int_impl! { usize, usize }
 
 /// Implementation of `RangeImpl` for float types.
 #[derive(Clone, Copy, Debug)]
-pub struct RangeFloat<X> {
-    low: X,
-    range: X,
+pub enum RangeFloat<X> {
+    // A range that is completely positive
+    Positive { offset: X, scale: X },
+    // A range that is completely negative
+    Negative { offset: X, scale: X },
+    // A range that is goes from negative to positive, but has more positive
+    // than negative numbers
+    PosAndNeg { cutoff: X, scale: X },
+    // A range that is goes from negative to positive, but has more negative
+    // than positive numbers
+    NegAndPos { cutoff: X, scale: X },
 }
 
 macro_rules! range_float_impl {
-    ($ty:ty) => {
+    ($ty:ty, $next_u:path) => {
         impl SampleRange for $ty {
             type T = RangeFloat<$ty>;
         }
-        
+
         impl RangeImpl for RangeFloat<$ty> {
+            // We can distinguish between two different kinds of ranges:
+            //
+            // Completely positive or negative.
+            // A characteristic of these ranges is that they get more bits of
+            // precision as they get closer to 0.0. For positive ranges it is as
+            // simple as applying a scale and offset to get the right range.
+            // For a negative range, we apply a negative scale to transform the
+            // range [0,1) to (-x,0]. Because this results in a range with it
+            // closed and open sides reversed, we do not sample from
+            // `closed_open01` but from `open_closed01`.
+            //
+            // A range that is both negative and positive.
+            // This range has as characteristic that it does not have most of
+            // its bits of precision near its low or high end, but in the middle
+            // near 0.0. As the basis we use the [-1,1) range from
+            // `closed_open11`.
+            // How it works is easiest to explain with an example.
+            // Suppose we want to sample from the range [-20, 100). We are
+            // first going to scale the range from [-1,1) to (-60,60]. Note the
+            // range changes from closed-open to open-closed because the scale
+            // is negative.
+            // If the sample happens to fall in the part (-60,-20), move it to
+            // (-100,60). Then multiply by -1.
+            // Scaling keeps the bits of precision intact. Moving the assymetric
+            // part also does not waste precision, as the more you move away
+            // from zero the less precision can be stored.
+
             type X = $ty;
-            
+
             fn new(low: Self::X, high: Self::X) -> Self {
-                RangeFloat {
-                    low: low,
-                    range: high - low,
+                if low >= 0.0 {
+                    RangeFloat::Positive {
+                        offset: low,
+                        scale: high - low,
+                    }
+                } else if high <= 0.0 {
+                    RangeFloat::Negative {
+                        offset: high,
+                        scale: low - high,
+                    }
+                } else if -low <= high {
+                    RangeFloat::PosAndNeg {
+                        cutoff: low,
+                        scale: (high + low) / -2.0 + low,
+                    }
+                } else {
+                    RangeFloat::NegAndPos {
+                        cutoff: high,
+                        scale: (high + low) / 2.0 - high,
+                    }
                 }
             }
-            
+
             fn sample<R: Rng+?Sized>(&self, rng: &mut R) -> Self::X {
-                let x: $ty = Rand::rand(rng, Uniform01);
-                self.low + self.range * x
+                let rnd = $next_u(rng);
+                match *self {
+                    RangeFloat::Positive { offset, scale } => {
+                        let x: $ty = rnd.closed_open01();
+                        offset + scale * x
+                    }
+                    RangeFloat::Negative { offset, scale } => {
+                        let x: $ty = rnd.open_closed01();
+                        offset + scale * x
+                    }
+                    RangeFloat::PosAndNeg { cutoff, scale } => {
+                        let x: $ty = rnd.closed_open11();
+                        let x = x * scale;
+                        if x < cutoff {
+                            (cutoff - scale) - x
+                        } else {
+                            x
+                        }
+                    }
+                    RangeFloat::NegAndPos { cutoff, scale } => {
+                        let x: $ty = rnd.closed_open11();
+                        let x = x * scale;
+                        if x >= cutoff {
+                            (cutoff + scale) - x
+                        } else {
+                            x
+                        }
+                    }
+                }
             }
         }
     }
 }
 
-range_float_impl! { f32 }
-range_float_impl! { f64 }
+range_float_impl! { f32, Rng::next_u32 }
+range_float_impl! { f64, Rng::next_u64 }
 
 #[cfg(test)]
 mod tests {