Skip to content
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

Tweak Open01 and Closed01 #237

Closed
wants to merge 3 commits into from
Closed

Tweak Open01 and Closed01 #237

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
cd830e7
Fix off-by-one error in mantissa_bits passed to float_impls
tspiteri Jan 15, 2018
cd18f48
remove bias in Open01
tspiteri Jan 15, 2018
a5da781
add tests for floating-point edge cases
tspiteri Feb 2, 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
49 changes: 39 additions & 10 deletions src/rand_impls.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -114,6 +114,7 @@ macro_rules! float_impls {
mod $mod_name {
use {Rand, Rng, Open01, Closed01};

// 1.0 / epsilon
const SCALE: $ty = (1u64 << $mantissa_bits) as $ty;

impl Rand for $ty {
Expand All @@ -130,11 +131,10 @@ macro_rules! float_impls {
impl Rand for Open01<$ty> {
#[inline]
fn rand<R: Rng>(rng: &mut R) -> Open01<$ty> {
// add a small amount (specifically 2 bits below
// 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() + 0.25 / SCALE)
Copy link
Contributor

Choose a reason for hiding this comment

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

This method was correct.
0.0 + 0.25 / SCALE will change the range to have 0.25 / SCALE as lower bound.
(1.0 - 1.0 / SCALE) + 0.25 / SCALE is still (1.0 - 1.0 / SCALE), thanks to the more limited precision near 1.0. At least, it works when the floating point rounding mode is set to nearest, and not up.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

This is a little more subjective than for the other change (that's why I put it in a second commit). The change in the first commit (changing 24, 53 into 23, 52) is to make the code work as was intended. The change from 0.25 / SCALE (ɛ/4) into 0.5 / SCALE (ɛ/2) is to change the intended behavior. Adding ɛ/4 will change from the half-open range [0,1) into an open range (0,1), but it is not symmetric about 0.5, so it is biased. Also, the spacing between all the possible outputs is not equal any more, as I'll try to show below.

The value obtained from rng can have any of the possible value:

0, ɛ, 2ɛ, …, 0.5−ɛ, 0.5, 0.5+ɛ, …, 1−2ɛ, 1−ɛ.

Adding ɛ/4 (the intended current behavior) would give the following possible outputs:

0.25ɛ, 1.25ɛ, 2.25ɛ, …, 0.5−0.75ɛ, 0.5+0.25ɛ, 0.5+1.25ɛ, …, 1−1.75ɛ, 1−0.75ɛ,

but rounded. Rounding has the effect that 0.25ɛ is added to all values < 0.5, but nothing is added to values ≥ 0.5, as then the added bit does not fit in the precision and is rounded off, so the actual outputs are:

0.25ɛ, 1.25ɛ, 2.25ɛ, …, 0.5−0.75ɛ, 0.5, 0.5+ɛ, …, 1−2ɛ, 1−ɛ.

Adding ɛ/2 to the value from rng would give the following possible outputs instead:

0.5ɛ, 1.5ɛ, 2.5ɛ, …, 0.5−0.5ɛ, 0.5+0.5ɛ, 0.5+1.5ɛ, …, 1−1.5ɛ, 1−0.5ɛ.

All these values can be represented exactly, so no rounding takes place. I see two advantages to this change:

  1. The expected value (average) is now exactly 0.5, which is what I would expect from a random number in the open range (0,1).
  2. The gap between possible values is always equal to ɛ, unlike the current intended behavior where the gap is almost always equal to ɛ, but is 0.75ɛ in the middle (between 0.5−0.75ɛ and 0.5).

// 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.
Open01(rng.$method_name() + 0.5 / SCALE)
}
}
impl Rand for Closed01<$ty> {
Expand All @@ -148,8 +148,8 @@ macro_rules! float_impls {
}
}
}
float_impls! { f64_rand_impls, f64, 53, next_f64 }
float_impls! { f32_rand_impls, f32, 24, next_f32 }
Copy link
Contributor

Choose a reason for hiding this comment

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

Floats only have the number of mantissa bits as you corrected, and one implicit bit. But does the logic here not rely on that? Are you sure the SCALE calculated in the macro is correct now?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

The comment for Closed01 is "rescale so that 1.0 - epsilon becomes 1.0 precisely". The rescaling is done by multiplying by SCALE / (SCALE - 1.0). Now for SCALE / (SCALE - 1.0) to rescale 1.0 - EPSILON into 1.0, we need SCALE / (SCALE - 1.0) == 1.0 / (1.0 - EPSILON), which is true when SCALE == 1.0 / EPSILON. For f32 ɛ is 2−23, so 1/ɛ (SCALE) should be equal to 1<<23, not 1<<24. For f64 ɛ is 2−52, so 1/ɛ (SCALE) should be equal to 1<<52, not 1<<53.

float_impls! { f64_rand_impls, f64, 52, next_f64 }
float_impls! { f32_rand_impls, f32, 23, next_f32 }

impl Rand for char {
#[inline]
Expand Down Expand Up @@ -269,9 +269,38 @@ mod tests {

#[test]
fn floating_point_edge_cases() {
// the test for exact equality is correct here.
assert!(ConstantRng(0xffff_ffff).gen::<f32>() != 1.0);
assert!(ConstantRng(0xffff_ffff_ffff_ffff).gen::<f64>() != 1.0);
const EPSILON32: f32 = 1.0 / (1u32 << 23) as f32;
const EPSILON64: f64 = 1.0 / (1u64 << 52) as f64;
Copy link
Member

Choose a reason for hiding this comment

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

Why not use core::f32::EPSILON?

ε is defined as the difference between 1.0 and the next largest representable number; since all numbers output (excepting Closed01) are strictly less than 1.0 the exponent will be 1 less than for 1.0, hence the largest representable number less than 1 should be 1 - ε / 2.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I first used std::f32::EPSILON but that failed the continuous integration tests and I forgot about core, so I'll change this.


let mut zeros = ConstantRng(0);
let mut ones = ConstantRng(!0);

let zero32 = zeros.gen::<f32>();
let zero64 = zeros.gen::<f64>();
let one32 = ones.gen::<f32>();
let one64 = ones.gen::<f64>();
assert_eq!(zero32, 0.0);
assert_eq!(zero64, 0.0);
assert!(1.0 - EPSILON32 <= one32 && one32 < 1.0);
assert!(1.0 - EPSILON64 <= one64 && one64 < 1.0);
Copy link
Member

Choose a reason for hiding this comment

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

As noted above, I believe this should be exactly 1 - ε / 2. Exact equality tests here would be preferable (excepting the Open01 case near 0, where there are many near representable values).

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Actually, one32 should be 1.0 - f32::EPSILON, not 1.0 - f32::EPSILON / 2, since the generator generates a multiple of ɛ in gen::<f32>. If we are to use the exact equality tests, they should be (zero32, one32) == (0.0, 1.0 - f32::EPSILON), (closed_zero32, closed_one32) == (0.0, 1.0) and (open_zero32, open_one32) == (f32::EPSILON / 2, 1.0 - f32::EPSILON / 2). Should I amend the commit to do that?


let Closed01(closed_zero32) = zeros.gen::<Closed01<f32>>();
let Closed01(closed_zero64) = zeros.gen::<Closed01<f64>>();
let Closed01(closed_one32) = ones.gen::<Closed01<f32>>();
let Closed01(closed_one64) = ones.gen::<Closed01<f64>>();
assert_eq!(closed_zero32, 0.0);
assert_eq!(closed_zero64, 0.0);
assert_eq!(closed_one32, 1.0);
assert_eq!(closed_one64, 1.0);

let Open01(open_zero32) = zeros.gen::<Open01<f32>>();
let Open01(open_zero64) = zeros.gen::<Open01<f64>>();
let Open01(open_one32) = ones.gen::<Open01<f32>>();
let Open01(open_one64) = ones.gen::<Open01<f64>>();
assert!(0.0 < open_zero32 && open_zero32 <= EPSILON32);
assert!(0.0 < open_zero64 && open_zero64 <= EPSILON64);
assert!(1.0 - EPSILON32 <= open_one32 && open_one32 < 1.0);
assert!(1.0 - EPSILON64 <= open_one64 && open_one64 < 1.0);
}

#[test]
Expand Down