psychopath/sub_crates/rmath/src/utils.rs

89 lines
2.9 KiB
Rust

const TOP_BIT: u32 = 1 << 31;
/// Compute how different two floats are in ulps.
///
/// Notes:
/// - Treats 0.0 and -0.0 as zero ulps apart, and extends the
/// implications of that to the rest of the numbers. E.g. the numbers
/// just above/below 0.0/-0.0 are only two ulps apart, not three.
/// - Infinity is one ulp past float max, and converse for -infinity.
/// - This doesn't handle NaNs in any really useful way.
#[inline(always)]
pub fn ulp_diff(a: f32, b: f32) -> u32 {
let a = a.to_bits();
let b = b.to_bits();
let a_sign = a & TOP_BIT;
let b_sign = b & TOP_BIT;
let a_abs = a & !TOP_BIT;
let b_abs = b & !TOP_BIT;
if a_sign == b_sign {
a_abs.max(b_abs) - a_abs.min(b_abs)
} else {
a_abs + b_abs
}
}
/// Checks if two floats are approximately equal, within `max_ulps`.
#[inline(always)]
pub fn ulps_eq(a: f32, b: f32, max_ulps: u32) -> bool {
!a.is_nan() && !b.is_nan() && (ulp_diff(a, b) <= max_ulps)
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn ulp_diff_test() {
assert_eq!(ulp_diff(1.0, 1.0), 0);
assert_eq!(ulp_diff(0.0, 0.0), 0);
assert_eq!(ulp_diff(0.0, -0.0), 0);
assert_eq!(ulp_diff(1.0, 2.0), 1 << 23);
assert_eq!(ulp_diff(2.0, 4.0), 1 << 23);
assert_eq!(ulp_diff(-1.0, -2.0), 1 << 23);
assert_eq!(ulp_diff(-2.0, -4.0), 1 << 23);
assert_eq!(ulp_diff(-1.0, 1.0), 0x7f000000);
assert_eq!(ulp_diff(0.0, 1.0), 0x3f800000);
assert_eq!(ulp_diff(-0.0, 1.0), 0x3f800000);
assert_eq!(ulp_diff(0.0, -1.0), 0x3f800000);
assert_eq!(ulp_diff(-0.0, -1.0), 0x3f800000);
assert_eq!(
ulp_diff(std::f32::INFINITY, -std::f32::INFINITY),
0xff000000
);
assert_eq!(ulp_diff(0.0, f32::from_bits(0.0f32.to_bits() + 1)), 1);
assert_eq!(ulp_diff(-0.0, f32::from_bits(0.0f32.to_bits() + 1)), 1);
}
#[test]
fn ulps_eq_test() {
assert!(ulps_eq(1.0, 1.0, 0));
assert!(ulps_eq(1.0, 1.0, 1));
assert!(ulps_eq(0.0, 0.0, 0));
assert!(ulps_eq(0.0, -0.0, 0));
assert!(ulps_eq(1.0, 2.0, 1 << 23));
assert!(!ulps_eq(1.0, 2.0, (1 << 23) - 1));
assert!(ulps_eq(0.0, f32::from_bits(0.0f32.to_bits() + 1), 1));
assert!(!ulps_eq(0.0, f32::from_bits(0.0f32.to_bits() + 1), 0));
assert!(ulps_eq(-0.0, f32::from_bits(0.0f32.to_bits() + 1), 1));
assert!(!ulps_eq(-0.0, f32::from_bits(0.0f32.to_bits() + 1), 0));
assert!(ulps_eq(std::f32::INFINITY, -std::f32::INFINITY, 0xff000000));
assert!(!ulps_eq(
std::f32::INFINITY,
-std::f32::INFINITY,
0xff000000 - 1
));
assert!(!ulps_eq(std::f32::NAN, std::f32::NAN, 0));
assert!(!ulps_eq(-std::f32::NAN, -std::f32::NAN, 0));
assert!(!ulps_eq(std::f32::NAN, std::f32::INFINITY, 1 << 31));
assert!(!ulps_eq(std::f32::INFINITY, std::f32::NAN, 1 << 31));
}
}