Skip to content

rand: Add next_f64/f32 to Rng. #18534

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.

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 1 commit into from
Nov 4, 2014
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
7 changes: 7 additions & 0 deletions src/librand/distributions/mod.rs
Original file line number Diff line number Diff line change
Expand Up @@ -227,6 +227,13 @@ fn ziggurat<R:Rng>(
// 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.
let bits: u64 = rng.gen();
let i = (bits & 0xff) as uint;
let f = (bits >> 11) as f64 / SCALE;
Expand Down
40 changes: 40 additions & 0 deletions src/librand/lib.rs
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,46 @@ pub trait Rng {
(self.next_u32() as u64 << 32) | (self.next_u32() as u64)
}

/// Return the next random f32 selected from the half-open
/// interval `[0, 1)`.
///
/// By default this is implemented in terms of `next_u32`, but a
/// random number generator which can generate numbers satisfying
/// the requirements directly can overload this for performance.
/// It is required that the return value lies in `[0, 1)`.
///
/// See `Closed01` for the closed interval `[0,1]`, and
/// `Open01` for the open interval `(0,1)`.
fn next_f32(&mut self) -> f32 {
const MANTISSA_BITS: uint = 24;
const IGNORED_BITS: uint = 8;
const SCALE: f32 = (1u64 << MANTISSA_BITS) as f32;

// using any more than `MANTISSA_BITS` bits will
// cause (e.g.) 0xffff_ffff to correspond to 1
// exactly, so we need to drop some (8 for f32, 11
// for f64) to guarantee the open end.
(self.next_u32() >> IGNORED_BITS) as f32 / SCALE
}

/// Return the next random f64 selected from the half-open
/// interval `[0, 1)`.
///
/// By default this is implemented in terms of `next_u64`, but a
/// random number generator which can generate numbers satisfying
/// the requirements directly can overload this for performance.
/// It is required that the return value lies in `[0, 1)`.
///
/// See `Closed01` for the closed interval `[0,1]`, and
/// `Open01` for the open interval `(0,1)`.
fn next_f64(&mut self) -> f64 {
const MANTISSA_BITS: uint = 53;
const IGNORED_BITS: uint = 11;
const SCALE: f64 = (1u64 << MANTISSA_BITS) as f64;

(self.next_u64() >> IGNORED_BITS) as f64 / SCALE
}

/// Fill `dest` with random data.
///
/// This has a default implementation in terms of `next_u64` and
Expand Down
23 changes: 9 additions & 14 deletions src/librand/rand_impls.rs
Original file line number Diff line number Diff line change
Expand Up @@ -96,11 +96,11 @@ impl Rand for u64 {
}

macro_rules! float_impls {
($mod_name:ident, $ty:ty, $mantissa_bits:expr, $method_name:ident, $ignored_bits:expr) => {
($mod_name:ident, $ty:ty, $mantissa_bits:expr, $method_name:ident) => {
mod $mod_name {
use {Rand, Rng, Open01, Closed01};

static SCALE: $ty = (1u64 << $mantissa_bits) as $ty;
const SCALE: $ty = (1u64 << $mantissa_bits) as $ty;

impl Rand for $ty {
/// Generate a floating point number in the half-open
Expand All @@ -110,11 +110,7 @@ macro_rules! float_impls {
/// and `Open01` for the open interval `(0,1)`.
#[inline]
fn rand<R: Rng>(rng: &mut R) -> $ty {
// using any more than `mantissa_bits` bits will
// cause (e.g.) 0xffff_ffff to correspond to 1
// exactly, so we need to drop some (8 for f32, 11
// for f64) to guarantee the open end.
(rng.$method_name() >> $ignored_bits) as $ty / SCALE
rng.$method_name()
}
}
impl Rand for Open01<$ty> {
Expand All @@ -124,23 +120,22 @@ macro_rules! float_impls {
// the precision of f64/f32 at 1.0), so that small
// numbers are larger than 0, but large numbers
// aren't pushed to/above 1.
Open01(((rng.$method_name() >> $ignored_bits) as $ty + 0.25) / SCALE)
Open01(rng.$method_name() + 0.25 / SCALE)
Copy link

Choose a reason for hiding this comment

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

I'm not knowledgeable about the topic but does this in any way skew random streams that are otherwise evenly distributed?

Copy link
Member Author

Choose a reason for hiding this comment

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

It does; there is precisely one value that doubles in frequency: 0.25. Let eps be the floating point precision, e.g. 2-53 for f64 and x = eps/4 be the value added here, then 0.25 = 0.25 + x = 0.25' + x where 0.25' = 0.25 - x is the largest float less than 0.25.

However, this is doubling from 1 in 2-53 to 1 in 2-52 (or 24 → 23 for f32) which is rather small. It's exceedingly unlikely that a consumer of f64's will see 0.25 even once (but rather less unlikely for f32's).

I think I investigated this before for the original implementation to see what was done elsewhere, but I have forgot it, and did not record my findings.

(Note that this line isn't actually changing the functionality.)

Copy link
Member Author

Choose a reason for hiding this comment

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

dSFMT operates on numbers between 1.0 and 2.0, just ors with 1 in the smallest bit and then subtracts 1.0 to lie in the open interval (0.0, 1.0) (e.g. 0 becomes 2-53). That behaviour halves the output space, but it does not lead to a bias.

GSL does rejection sampling (calling the RNG until 0 is not encountered, which is rather likely). I guess we should consider doing that.

Copy link

Choose a reason for hiding this comment

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

Thanks @huonw!

}
}
impl Rand for Closed01<$ty> {
#[inline]
fn rand<R: Rng>(rng: &mut R) -> Closed01<$ty> {
// divide by the maximum value of the numerator to
// get a non-zero probability of getting exactly
// 1.0.
Closed01((rng.$method_name() >> $ignored_bits) as $ty / (SCALE - 1.0))
// rescale so that 1.0 - epsilon becomes 1.0
// precisely.
Closed01(rng.$method_name() * SCALE / (SCALE - 1.0))
}
}
}
}
}
float_impls! { f64_rand_impls, f64, 53, next_u64, 11 }
float_impls! { f32_rand_impls, f32, 24, next_u32, 8 }
float_impls! { f64_rand_impls, f64, 53, next_f64 }
float_impls! { f32_rand_impls, f32, 24, next_f32 }

impl Rand for char {
#[inline]
Expand Down