Skip to content

Port new Range implementation and only have one uniform float distribution #274

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 8 commits into from
Mar 3, 2018
Merged

Port new Range implementation and only have one uniform float distribution #274

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
741c310
Port new Range implementation and only have one uniform float distrib…
pitdicker Feb 22, 2018
526667e
Use `Range` in `Distribution<char>`
pitdicker Feb 28, 2018
68bb4aa
Generate bool using sign test
pitdicker Feb 28, 2018
4a9c465
Update `Uniform` documentation
pitdicker Feb 28, 2018
616861b
Documentation fixes
pitdicker Mar 3, 2018
c90fee7
Correct Range compatability code
pitdicker Mar 3, 2018
783277a
Remove unnecessary range tests
pitdicker Mar 3, 2018
7b968e8
Fix spelling
pitdicker Mar 3, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 4 additions & 9 deletions benches/distributions.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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));
Expand Down
229 changes: 56 additions & 173 deletions src/distributions/float.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,203 +14,86 @@ 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
//
// * <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/dSFMT.pdf>
// * <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/dSFMT-slide-e.pdf>
impl Distribution<f32> for Uniform {
fn sample<R: Rng + ?Sized>(&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<f64> for Uniform {
fn sample<R: Rng + ?Sized>(&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
/// [2<sup>0</sup>..2<sup>1</sup>), 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<R: Rng + ?Sized>(&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<R: Rng + ?Sized>(&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<R: Rng + ?Sized>(&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)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If 1+ε is the smallest representable number above 1, then 1-ε/2 is representable; ok. This is the same adjustment as the Open01 removed here but in a single number. Looks fine and functionally identical.

}
}
}
}
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;

#[test]
fn floating_point_edge_cases() {
let mut zeros = StepRng::new(0, 0);
assert_eq!(zeros.gen::<f32>(), 0.0);
assert_eq!(zeros.gen::<f64>(), 0.0);

let mut one = StepRng::new(1, 0);
assert_eq!(one.gen::<f32>(), EPSILON32);
assert_eq!(one.gen::<f64>(), EPSILON64);

let mut max = StepRng::new(!0, 0);
assert_eq!(max.gen::<f32>(), 1.0 - EPSILON32);
assert_eq!(max.gen::<f64>(), 1.0 - EPSILON64);
}
assert_eq!(zeros.gen::<f32>(), 0.0 + EPSILON32 / 2.0);
assert_eq!(zeros.gen::<f64>(), 0.0 + EPSILON64 / 2.0);

#[test]
fn fp_closed_edge_cases() {
let mut zeros = StepRng::new(0, 0);
assert_eq!(zeros.sample::<f32, _>(Closed01), 0.0);
assert_eq!(zeros.sample::<f64, _>(Closed01), 0.0);

let mut one = StepRng::new(1, 0);
let one32 = one.sample::<f32, _>(Closed01);
let one64 = one.sample::<f64, _>(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::<f32, _>(Closed01), 1.0);
assert_eq!(max.sample::<f64, _>(Closed01), 1.0);
}

#[test]
fn fp_open_edge_cases() {
let mut zeros = StepRng::new(0, 0);
assert_eq!(zeros.sample::<f32, _>(Open01), 0.0 + EPSILON32 / 2.0);
assert_eq!(zeros.sample::<f64, _>(Open01), 0.0 + EPSILON64 / 2.0);

let mut one = StepRng::new(1, 0);
let one32 = one.sample::<f32, _>(Open01);
let one64 = one.sample::<f64, _>(Open01);
let mut one = StepRng::new(1 << 9, 0);
let one32 = one.gen::<f32>();
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::<f32, _>(Open01), 1.0 - EPSILON32 / 2.0);
assert_eq!(max.sample::<f64, _>(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::<f64>();
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::<f32>(), 1.0 - EPSILON32 / 2.0);
assert_eq!(max.gen::<f64>(), 1.0 - EPSILON64 / 2.0);
}
}
6 changes: 3 additions & 3 deletions src/distributions/gamma.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
///
Expand Down Expand Up @@ -144,7 +144,7 @@ impl Distribution<f64> for Gamma {
}
impl Distribution<f64> for GammaSmallShape {
fn sample<R: Rng + ?Sized>(&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)
}
Expand All @@ -159,7 +159,7 @@ impl Distribution<f64> 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 ||
Expand Down
Loading