diff --git a/src/distributions/float.rs b/src/distributions/float.rs index 522c9f3270..7f205e5a30 100644 --- a/src/distributions/float.rs +++ b/src/distributions/float.rs @@ -77,19 +77,19 @@ pub struct HighPrecision where F: HPFloatHelper { low_as_int: F::SignedInt, high_as_int: F::SignedInt, exponent: u16, - distribution: F::SignedIntDistribution, + mantissa_distribution: F::SignedIntDistribution, } impl HighPrecision { /// Create a new HighPrecision distribution. Sampling from this - /// distribution will return values `>= low` and `< high`. + /// distribution will return values `>= low` and `< high`. pub fn new(low: F, high: F) -> Self { let parsed = F::parse_new(low, high); HighPrecision { low_as_int: parsed.0, high_as_int: parsed.1, exponent: parsed.2, - distribution: parsed.3, + mantissa_distribution: parsed.3, } } } @@ -239,7 +239,7 @@ macro_rules! float_impls { if bits >= ::core::mem::size_of::<$uty>() as $uty * 8 { (-1 as $ity) as $uty } else { (1 as $uty << bits) - 1 } } loop { - let signed_mant = self.distribution.sample(rng); + let signed_mant = self.mantissa_distribution.sample(rng); // Operate on the absolute value so that we can count bit-sizes // correctly let is_neg = signed_mant < 0; @@ -302,6 +302,10 @@ macro_rules! float_impls { return <$ty>::from_bits(res); } + // If not start over. We're avoiding reusing any of the previous + // computation in order to avoid introducing bias, and to keep + // things simple since this should be rare. + // Assert that we got here due to rounding #[cfg(debug_assertions)] { @@ -318,15 +322,9 @@ macro_rules! float_impls { assert!(mant_high & bitmask(exp_low - exp_high) != 0); } } - - // If not start over. We're avoiding reusing any of the previous - // computation in order to avoid introducing bias, and to keep - // things simple since this should be rare. } } } - - } } float_impls! { f32, u32, i32, 23, 8, 127 } @@ -399,115 +397,174 @@ mod tests { let mut r = ::test::rng(601); macro_rules! float_test { - ($ty:ty, $uty:ty, $ity:ty, $test_vals:expr) => {{ - let mut vals: Vec<$ty> = - $test_vals.iter().cloned() - .flat_map(|x| (-2 as $ity..3).map(move |y| x + y)) - .map(|x| <$ty>::from_bits(x as $uty)) - .flat_map(|x| vec![x, -x].into_iter()) - .filter(|x| x.is_finite()) - .collect(); - vals.sort_by(|a, b| a.partial_cmp(b).unwrap()); - vals.dedup(); - - for a in vals.iter().cloned() { - for b in vals.iter().cloned().filter(|&b| b > a) { - fn to_signed_bits(val: $ty) -> $ity { - if val >= 0.0 { - val.to_bits() as $ity - } else { - -((-val).to_bits() as $ity) + ($ty:ty, $uty:ty, $ity:ty, $extra:expr, $test_vals:expr) => { + // Create a closure to make loop labels local + (|| { + let mut vals: Vec<$ty> = + $test_vals.iter().cloned() + .flat_map(|x| $extra.iter().map(move |y| x + y)) + .map(|x| <$ty>::from_bits(x as $uty)) + .flat_map(|x| vec![x, -x].into_iter()) + .filter(|x| x.is_finite()) + .collect(); + vals.sort_by(|a, b| a.partial_cmp(b).unwrap()); + vals.dedup(); + + for a in vals.iter().cloned() { + for b in vals.iter().cloned().filter(|&b| b > a) { + fn to_signed_bits(val: $ty) -> $ity { + if val >= 0.0 { + val.to_bits() as $ity + } else { + -((-val).to_bits() as $ity) + } } - } - fn from_signed_bits(val: $ity) -> $ty { - if val >= 0 { - <$ty>::from_bits(val as $uty) - } else { - -<$ty>::from_bits(-val as $uty) + fn from_signed_bits(val: $ity) -> $ty { + if val >= 0 { + <$ty>::from_bits(val as $uty) + } else { + -<$ty>::from_bits(-val as $uty) + } } - } - let hp = HighPrecision::new(a, b); - let a_bits = to_signed_bits(a); - let b_bits = to_signed_bits(b); - - // If a and b are "close enough", we can verify the full distribution - if (b_bits.wrapping_sub(a_bits) as $uty) < 100 { - let mut counts = Vec::::with_capacity((b_bits - a_bits) as usize); - counts.resize((b_bits - a_bits) as usize, 0); - for _ in 0..1000 { - let res = hp.sample(&mut r); - counts[(to_signed_bits(res) - a_bits) as usize] += 1; - } - for (count, i) in counts.iter().zip(0 as $ity..) { - let expected = 1000.0 as $ty * - (from_signed_bits(a_bits + i + 1) - - from_signed_bits(a_bits + i)) / (b - a); - let err = (*count as $ty - expected) / expected; - assert!(err.abs() <= 0.2); - } - } else { - // Rough estimate that the distribution is correct - let step = if (b - a).is_finite() { - (b - a) / 10.0 + let hp = HighPrecision::new(a, b); + let a_bits = to_signed_bits(a); + let b_bits = to_signed_bits(b); + + const N_RUNS: usize = 10; + const N_REPS_PER_RUN: usize = 1000; + + if (b_bits.wrapping_sub(a_bits) as $uty) < 100 { + // If a and b are "close enough", we can verify the full distribution + let mut counts = Vec::::with_capacity((b_bits - a_bits) as usize); + counts.resize((b_bits - a_bits) as usize, 0); + 'test_loop_exact: for test_run in 1..(N_RUNS+1) { + for _ in 0..N_REPS_PER_RUN { + let res = hp.sample(&mut r); + counts[(to_signed_bits(res) - a_bits) as usize] += 1; + } + for (count, i) in counts.iter().zip(0 as $ity..) { + let expected = (test_run * N_REPS_PER_RUN) as $ty * + ((from_signed_bits(a_bits + i + 1) - + from_signed_bits(a_bits + i)) / (b - a)); + let err = (*count as $ty - expected) / expected; + if err.abs() > 0.2 { + if test_run < N_RUNS { + continue 'test_loop_exact; + } + panic!(format!("Failed {}-bit exact test: a: 0x{:x}, b: 0x{:x}, err: {:.2}", + ::core::mem::size_of::<$ty>() * 8, + a.to_bits(), + b.to_bits(), + err.abs())); + } + } + } } else { - b / 10.0 - a / 10.0 - }; - assert!(step.is_finite()); - let mut counts = Vec::::with_capacity(10); - counts.resize(10, 0); - for _ in 0..3000 { - let res = hp.sample(&mut r); - assert!(a <= res && res < b); - let index = if (res - a).is_finite() { - (res - a) / step + // Otherwise divide range into 10 sections + let step = if (b - a).is_finite() { + (b - a) / 10.0 } else { - res / step - a / step - } as usize; - counts[::core::cmp::min(index, 9)] += 1; - } - for count in &counts { - let expected = 3000.0 as $ty / 10.0; - let err = (*count as $ty - expected) / expected; - assert!(err.abs() <= 0.25); + b / 10.0 - a / 10.0 + }; + assert!(step.is_finite()); + let mut counts = Vec::::with_capacity(10); + counts.resize(10, 0); + + 'test_loop_rough: for test_run in 1..(N_RUNS+1) { + for _ in 0..N_REPS_PER_RUN { + let res = hp.sample(&mut r); + assert!(a <= res && res < b); + let index = (res / step - a / step) as usize; + counts[::core::cmp::min(index, 9)] += 1; + } + for count in &counts { + let expected = (test_run * N_REPS_PER_RUN) as $ty / 10.0; + let err = (*count as $ty - expected) / expected; + if err.abs() > 0.2 { + if test_run < N_RUNS { + continue 'test_loop_rough; + } + panic!(format!("Failed {}-bit rough test: a: 0x{:x}, b: 0x{:x}, err: {:.2}", + ::core::mem::size_of::<$ty>() * 8, + a.to_bits(), + b.to_bits(), + err.abs())); + } + } + } } } } - } - }} + })() + } } - float_test!(f64, u64, i64, - [0i64, - 0x0000_0f00_0000_0000, - 0x0001_0000_0000_0000, - 0x0004_0000_0000_0000, - 0x0008_0000_0000_0000, - 0x0010_0000_0000_0000, - 0x0020_0000_0000_0000, - 0x0040_0000_0000_0000, - 0x0100_0000_0000_0000, - 0x00cd_ef12_3456_789a, - 0x0100_ffff_ffff_ffff, - 0x010f_ffff_ffff_ffff, - 0x0400_1234_5678_abcd, - 0x7fef_ffff_ffff_ffff, - ]); - float_test!(f32, u32, i32, - [0i32, - 0x000f_0000, - 0x0008_0000, - 0x0020_0000, - 0x0040_0000, - 0x0080_0000, - 0x0100_0000, - 0x0200_0000, - 0x0800_0000, - 0x5678_abcd, - 0x0807_ffff, - 0x087f_ffff, - 0x4012_3456, - 0x7f7f_ffff, - ]); + const SLOW_TESTS: bool = false; + if SLOW_TESTS { + // These test cases are commented out since they + // take too long to run. + float_test!(f64, u64, i64, + [-5, -1, 0, 1, 7], + [0i64, + 0x0000_0f00_0000_0000, + 0x0001_0000_0000_0000, + 0x0004_0000_0000_0000, + 0x0008_0000_0000_0000, + 0x0010_0000_0000_0000, + 0x0020_0000_0000_0000, + 0x0040_0000_0000_0000, + 0x0100_0000_0000_0000, + 0x00cd_ef12_3456_789a, + 0x0100_ffff_ffff_ffff, + 0x010f_ffff_ffff_ffff, + 0x0400_1234_5678_abcd, + 0x7fef_ffff_ffff_ffff, + ]); + float_test!(f32, u32, i32, + [-5, -1, 0, 1, 7], + [0i32, + 0x000f_0000, + 0x0008_0000, + 0x0020_0000, + 0x0040_0000, + 0x0080_0000, + 0x0100_0000, + 0x0200_0000, + 0x0800_0000, + 0x5678_abcd, + 0x0807_ffff, + 0x087f_ffff, + 0x4012_3456, + 0x7f7f_ffff, + ]); + } else { + float_test!(f64, u64, i64, + [0], + [0i64, + 1, + 0x0000_0f00_0000_0000, + 0x0000_0f00_0000_0005, + 0x000f_ffff_ffff_fffd, + 0x0010_0000_0000_0000, + 0x0040_0000_0000_0000, + 0x0100_ffff_ffff_ffff, + 0x0101_0000_0000_0004, + 0x7fef_ffff_ffff_ffff, + ]); + float_test!(f32, u32, i32, + [0], + [0i32, + 1, + 0x000f_0000, + 0x000f_0005, + 0x007f_fffd, + 0x0080_0000, + 0x0200_0000, + 0x0807_ffff, + 0x0808_0004, + 0x7f7f_ffff, + ]); + } } }