From e80dad0eb512ad5c350d758ec108769bac3fdb55 Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Sun, 29 Dec 2024 07:37:01 +0000 Subject: [PATCH 1/9] Always enable `unstable-float` in CI Since these add new API but do not affect runtime, we can enable it for all tests that run with nightly. --- ci/run.sh | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/ci/run.sh b/ci/run.sh index d89c8bdf0..7e514a1cd 100755 --- a/ci/run.sh +++ b/ci/run.sh @@ -62,22 +62,26 @@ esac cargo check --no-default-features cargo check --features "force-soft-floats" +# Always enable `unstable-float` since it expands available API but does not +# change any implementations. +extra_flags="$extra_flags --features unstable-float" + if [ "${BUILD_ONLY:-}" = "1" ]; then cmd="cargo build --target $target --package libm" $cmd - $cmd --features "unstable-intrinsics" + $cmd --features unstable-intrinsics echo "can't run tests on $target; skipping" else cmd="cargo test --all --target $target $extra_flags" - # stable by default + # Test without intrinsics $cmd $cmd --release - # unstable with a feature - $cmd --features "unstable-intrinsics" - $cmd --release --features "unstable-intrinsics" + # Test with intrinsic use + $cmd --features unstable-intrinsics + $cmd --release --features unstable-intrinsics # Make sure benchmarks have correct results $cmd --benches From b1b30caa2111131b06b9d773a7ba2fb79edf4a7f Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Thu, 19 Dec 2024 11:47:04 +0000 Subject: [PATCH 2/9] Update and slightly refactor some of the `Float` trait Add a constant for negative pi and provide a standalone const `from_bits`, which can be combined with what we already had in `hex_float`. Also provide another default method to reduce what needs to be provided by the macro. --- src/math/support/float_traits.rs | 47 ++++++++++++++++++++++---------- src/math/support/hex_float.rs | 12 ++------ src/math/support/mod.rs | 1 + 3 files changed, 36 insertions(+), 24 deletions(-) diff --git a/src/math/support/float_traits.rs b/src/math/support/float_traits.rs index 7b3f6904b..68ba60030 100644 --- a/src/math/support/float_traits.rs +++ b/src/math/support/float_traits.rs @@ -38,6 +38,7 @@ pub trait Float: const MAX: Self; const MIN: Self; const PI: Self; + const NEG_PI: Self; const FRAC_PI_2: Self; /// The bitwidth of the float type @@ -71,7 +72,9 @@ pub trait Float: fn to_bits(self) -> Self::Int; /// Returns `self` transmuted to `Self::SignedInt` - fn to_bits_signed(self) -> Self::SignedInt; + fn to_bits_signed(self) -> Self::SignedInt { + self.to_bits().signed() + } /// Checks if two floats have the same bit representation. *Except* for NaNs! NaN can be /// represented in multiple different ways. This method returns `true` if two NaNs are @@ -158,7 +161,15 @@ pub trait Float: pub type IntTy = ::Int; macro_rules! float_impl { - ($ty:ident, $ity:ident, $sity:ident, $expty:ident, $bits:expr, $significand_bits:expr) => { + ( + $ty:ident, + $ity:ident, + $sity:ident, + $expty:ident, + $bits:expr, + $significand_bits:expr, + $from_bits:path + ) => { impl Float for $ty { type Int = $ity; type SignedInt = $sity; @@ -173,13 +184,10 @@ macro_rules! float_impl { const NAN: Self = Self::NAN; const MAX: Self = -Self::MIN; // Sign bit set, saturated mantissa, saturated exponent with last bit zeroed - // FIXME(msrv): just use `from_bits` when available - // SAFETY: POD cast with no preconditions - const MIN: Self = unsafe { - mem::transmute::(Self::Int::MAX & !(1 << Self::SIG_BITS)) - }; + const MIN: Self = $from_bits(Self::Int::MAX & !(1 << Self::SIG_BITS)); const PI: Self = core::$ty::consts::PI; + const NEG_PI: Self = -Self::PI; const FRAC_PI_2: Self = core::$ty::consts::FRAC_PI_2; const BITS: u32 = $bits; @@ -193,9 +201,6 @@ macro_rules! float_impl { fn to_bits(self) -> Self::Int { self.to_bits() } - fn to_bits_signed(self) -> Self::SignedInt { - self.to_bits() as Self::SignedInt - } fn is_nan(self) -> bool { self.is_nan() } @@ -220,8 +225,22 @@ macro_rules! float_impl { } #[cfg(f16_enabled)] -float_impl!(f16, u16, i16, i8, 16, 10); -float_impl!(f32, u32, i32, i16, 32, 23); -float_impl!(f64, u64, i64, i16, 64, 52); +float_impl!(f16, u16, i16, i8, 16, 10, f16::from_bits); +float_impl!(f32, u32, i32, i16, 32, 23, f32_from_bits); +float_impl!(f64, u64, i64, i16, 64, 52, f64_from_bits); #[cfg(f128_enabled)] -float_impl!(f128, u128, i128, i16, 128, 112); +float_impl!(f128, u128, i128, i16, 128, 112, f128::from_bits); + +/* FIXME(msrv): vendor some things that are not const stable at our MSRV */ + +/// `f32::from_bits` +pub const fn f32_from_bits(bits: u32) -> f32 { + // SAFETY: POD cast with no preconditions + unsafe { mem::transmute::(bits) } +} + +/// `f64::from_bits` +pub const fn f64_from_bits(bits: u64) -> f64 { + // SAFETY: POD cast with no preconditions + unsafe { mem::transmute::(bits) } +} diff --git a/src/math/support/hex_float.rs b/src/math/support/hex_float.rs index 80434a5ec..1666c6153 100644 --- a/src/math/support/hex_float.rs +++ b/src/math/support/hex_float.rs @@ -2,6 +2,8 @@ #![allow(dead_code)] // FIXME: remove once this gets used +use super::{f32_from_bits, f64_from_bits}; + /// Construct a 32-bit float from hex float representation (C-style) pub const fn hf32(s: &str) -> f32 { f32_from_bits(parse_any(s, 32, 23) as u32) @@ -159,16 +161,6 @@ const fn hex_digit(c: u8) -> u8 { /* FIXME(msrv): vendor some things that are not const stable at our MSRV */ -/// `f32::from_bits` -const fn f32_from_bits(v: u32) -> f32 { - unsafe { core::mem::transmute(v) } -} - -/// `f64::from_bits` -const fn f64_from_bits(v: u64) -> f64 { - unsafe { core::mem::transmute(v) } -} - /// `u128::ilog2` const fn u128_ilog2(v: u128) -> u32 { assert!(v != 0); diff --git a/src/math/support/mod.rs b/src/math/support/mod.rs index 25681c307..e2f4e0e98 100644 --- a/src/math/support/mod.rs +++ b/src/math/support/mod.rs @@ -6,6 +6,7 @@ mod int_traits; #[allow(unused_imports)] pub use float_traits::{Float, IntTy}; +pub(crate) use float_traits::{f32_from_bits, f64_from_bits}; #[allow(unused_imports)] pub use hex_float::{hf32, hf64}; pub use int_traits::{CastFrom, CastInto, DInt, HInt, Int, MinInt}; From f0817728e58d4b0d6e8e03d06cab82ce6ec2b289 Mon Sep 17 00:00:00 2001 From: beetrees Date: Thu, 19 Dec 2024 11:59:09 +0000 Subject: [PATCH 3/9] Remove an `is_nan` workaround that is no longer needed --- src/math/support/float_traits.rs | 12 +----------- 1 file changed, 1 insertion(+), 11 deletions(-) diff --git a/src/math/support/float_traits.rs b/src/math/support/float_traits.rs index 68ba60030..e64640a0d 100644 --- a/src/math/support/float_traits.rs +++ b/src/math/support/float_traits.rs @@ -80,17 +80,7 @@ pub trait Float: /// represented in multiple different ways. This method returns `true` if two NaNs are /// compared. fn eq_repr(self, rhs: Self) -> bool { - let is_nan = |x: Self| -> bool { - // } - // fn is_nan(x: Self) -> bool { - // When using mangled-names, the "real" compiler-builtins might not have the - // necessary builtin (__unordtf2) to test whether `f128` is NaN. - // FIXME(f16_f128): Remove once the nightly toolchain has the __unordtf2 builtin - // x is NaN if all the bits of the exponent are set and the significand is non-0 - x.to_bits() & Self::EXP_MASK == Self::EXP_MASK - && x.to_bits() & Self::SIG_MASK != Self::Int::ZERO - }; - if is_nan(self) && is_nan(rhs) { true } else { self.to_bits() == rhs.to_bits() } + if self.is_nan() && rhs.is_nan() { true } else { self.to_bits() == rhs.to_bits() } } /// Returns true if the value is NaN. From 69cf64f8f5f5275802ca2d0c2ac3d7ec27ac3f6f Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Thu, 19 Dec 2024 11:18:33 +0000 Subject: [PATCH 4/9] Add an 8-bit float type for testing purposes Introduce `f8`, which is an 8-bit float compliant with IEEE-754. This type is useful for testing since it is easily possible to enumerate all values. --- crates/libm-test/src/f8_impl.rs | 487 ++++++++++++++++++++++++++++++++ crates/libm-test/src/lib.rs | 4 + 2 files changed, 491 insertions(+) create mode 100644 crates/libm-test/src/f8_impl.rs diff --git a/crates/libm-test/src/f8_impl.rs b/crates/libm-test/src/f8_impl.rs new file mode 100644 index 000000000..babcc6357 --- /dev/null +++ b/crates/libm-test/src/f8_impl.rs @@ -0,0 +1,487 @@ +//! An IEEE-compliant 8-bit float type for testing purposes. + +use std::cmp::{self, Ordering}; +use std::{fmt, ops}; + +use crate::Float; + +/// Sometimes verifying float logic is easiest when all values can quickly be checked exhaustively +/// or by hand. +/// +/// IEEE-754 compliant type that includes a 1 bit sign, 4 bit exponent, and 3 bit significand. +/// Bias is -7. +/// +/// Based on . +#[derive(Clone, Copy)] +#[repr(transparent)] +#[allow(non_camel_case_types)] +pub struct f8(u8); + +impl Float for f8 { + type Int = u8; + type SignedInt = i8; + type ExpInt = i8; + + const ZERO: Self = Self(0b0_0000_000); + const NEG_ZERO: Self = Self(0b1_0000_000); + const ONE: Self = Self(0b0_0111_000); + const NEG_ONE: Self = Self(0b1_0111_000); + const MAX: Self = Self(0b0_1110_111); + const MIN: Self = Self(0b1_1110_111); + const INFINITY: Self = Self(0b0_1111_000); + const NEG_INFINITY: Self = Self(0b1_1111_000); + const NAN: Self = Self(0b0_1111_100); + const PI: Self = Self::ZERO; + const NEG_PI: Self = Self::ZERO; + const FRAC_PI_2: Self = Self::ZERO; + + const BITS: u32 = 8; + const SIG_BITS: u32 = 3; + const SIGN_MASK: Self::Int = 0b1_0000_000; + const SIG_MASK: Self::Int = 0b0_0000_111; + const EXP_MASK: Self::Int = 0b0_1111_000; + const IMPLICIT_BIT: Self::Int = 0b0_0001_000; + + fn to_bits(self) -> Self::Int { + self.0 + } + + fn to_bits_signed(self) -> Self::SignedInt { + self.0 as i8 + } + + fn is_nan(self) -> bool { + self.0 & Self::EXP_MASK == Self::EXP_MASK && self.0 & Self::SIG_MASK != 0 + } + + fn is_infinite(self) -> bool { + self.0 & Self::EXP_MASK == Self::EXP_MASK && self.0 & Self::SIG_MASK == 0 + } + + fn is_sign_negative(self) -> bool { + self.0 & Self::SIGN_MASK != 0 + } + + fn exp(self) -> Self::ExpInt { + unimplemented!() + } + + fn from_bits(a: Self::Int) -> Self { + Self(a) + } + + fn normalize(_significand: Self::Int) -> (i32, Self::Int) { + unimplemented!() + } +} + +impl f8 { + pub const ALL_LEN: usize = 240; + + /// All non-infinite non-NaN values of `f8` + pub const ALL: [Self; Self::ALL_LEN] = [ + // -m*2^7 + Self(0b1_1110_111), // -240 + Self(0b1_1110_110), + Self(0b1_1110_101), + Self(0b1_1110_100), + Self(0b1_1110_011), + Self(0b1_1110_010), + Self(0b1_1110_001), + Self(0b1_1110_000), // -128 + // -m*2^6 + Self(0b1_1101_111), // -120 + Self(0b1_1101_110), + Self(0b1_1101_101), + Self(0b1_1101_100), + Self(0b1_1101_011), + Self(0b1_1101_010), + Self(0b1_1101_001), + Self(0b1_1101_000), // -64 + // -m*2^5 + Self(0b1_1100_111), // -60 + Self(0b1_1100_110), + Self(0b1_1100_101), + Self(0b1_1100_100), + Self(0b1_1100_011), + Self(0b1_1100_010), + Self(0b1_1100_001), + Self(0b1_1100_000), // -32 + // -m*2^4 + Self(0b1_1011_111), // -30 + Self(0b1_1011_110), + Self(0b1_1011_101), + Self(0b1_1011_100), + Self(0b1_1011_011), + Self(0b1_1011_010), + Self(0b1_1011_001), + Self(0b1_1011_000), // -16 + // -m*2^3 + Self(0b1_1010_111), // -15 + Self(0b1_1010_110), + Self(0b1_1010_101), + Self(0b1_1010_100), + Self(0b1_1010_011), + Self(0b1_1010_010), + Self(0b1_1010_001), + Self(0b1_1010_000), // -8 + // -m*2^2 + Self(0b1_1001_111), // -7.5 + Self(0b1_1001_110), + Self(0b1_1001_101), + Self(0b1_1001_100), + Self(0b1_1001_011), + Self(0b1_1001_010), + Self(0b1_1001_001), + Self(0b1_1001_000), // -4 + // -m*2^1 + Self(0b1_1000_111), // -3.75 + Self(0b1_1000_110), + Self(0b1_1000_101), + Self(0b1_1000_100), + Self(0b1_1000_011), + Self(0b1_1000_010), + Self(0b1_1000_001), + Self(0b1_1000_000), // -2 + // -m*2^0 + Self(0b1_0111_111), // -1.875 + Self(0b1_0111_110), + Self(0b1_0111_101), + Self(0b1_0111_100), + Self(0b1_0111_011), + Self(0b1_0111_010), + Self(0b1_0111_001), + Self(0b1_0111_000), // -1 + // -m*2^-1 + Self(0b1_0110_111), // −0.9375 + Self(0b1_0110_110), + Self(0b1_0110_101), + Self(0b1_0110_100), + Self(0b1_0110_011), + Self(0b1_0110_010), + Self(0b1_0110_001), + Self(0b1_0110_000), // -0.5 + // -m*2^-2 + Self(0b1_0101_111), // −0.46875 + Self(0b1_0101_110), + Self(0b1_0101_101), + Self(0b1_0101_100), + Self(0b1_0101_011), + Self(0b1_0101_010), + Self(0b1_0101_001), + Self(0b1_0101_000), // -0.25 + // -m*2^-3 + Self(0b1_0100_111), // −0.234375 + Self(0b1_0100_110), + Self(0b1_0100_101), + Self(0b1_0100_100), + Self(0b1_0100_011), + Self(0b1_0100_010), + Self(0b1_0100_001), + Self(0b1_0100_000), // -0.125 + // -m*2^-4 + Self(0b1_0011_111), // −0.1171875 + Self(0b1_0011_110), + Self(0b1_0011_101), + Self(0b1_0011_100), + Self(0b1_0011_011), + Self(0b1_0011_010), + Self(0b1_0011_001), + Self(0b1_0011_000), // −0.0625 + // -m*2^-5 + Self(0b1_0010_111), // −0.05859375 + Self(0b1_0010_110), + Self(0b1_0010_101), + Self(0b1_0010_100), + Self(0b1_0010_011), + Self(0b1_0010_010), + Self(0b1_0010_001), + Self(0b1_0010_000), // −0.03125 + // -m*2^-6 + Self(0b1_0001_111), // −0.029296875 + Self(0b1_0001_110), + Self(0b1_0001_101), + Self(0b1_0001_100), + Self(0b1_0001_011), + Self(0b1_0001_010), + Self(0b1_0001_001), + Self(0b1_0001_000), // −0.015625 + // -m*2^-7 subnormal numbers + Self(0b1_0000_111), // −0.013671875 + Self(0b1_0000_110), + Self(0b1_0000_101), + Self(0b1_0000_100), + Self(0b1_0000_011), + Self(0b1_0000_010), + Self(0b1_0000_001), // −0.001953125 + // Zeroes + Self(0b1_0000_000), // -0.0 + Self(0b0_0000_000), // 0.0 + // m*2^-7 // subnormal numbers + Self(0b0_0000_001), + Self(0b0_0000_010), + Self(0b0_0000_011), + Self(0b0_0000_100), + Self(0b0_0000_101), + Self(0b0_0000_110), + Self(0b0_0000_111), // 0.013671875 + // m*2^-6 + Self(0b0_0001_000), // 0.015625 + Self(0b0_0001_001), + Self(0b0_0001_010), + Self(0b0_0001_011), + Self(0b0_0001_100), + Self(0b0_0001_101), + Self(0b0_0001_110), + Self(0b0_0001_111), // 0.029296875 + // m*2^-5 + Self(0b0_0010_000), // 0.03125 + Self(0b0_0010_001), + Self(0b0_0010_010), + Self(0b0_0010_011), + Self(0b0_0010_100), + Self(0b0_0010_101), + Self(0b0_0010_110), + Self(0b0_0010_111), // 0.05859375 + // m*2^-4 + Self(0b0_0011_000), // 0.0625 + Self(0b0_0011_001), + Self(0b0_0011_010), + Self(0b0_0011_011), + Self(0b0_0011_100), + Self(0b0_0011_101), + Self(0b0_0011_110), + Self(0b0_0011_111), // 0.1171875 + // m*2^-3 + Self(0b0_0100_000), // 0.125 + Self(0b0_0100_001), + Self(0b0_0100_010), + Self(0b0_0100_011), + Self(0b0_0100_100), + Self(0b0_0100_101), + Self(0b0_0100_110), + Self(0b0_0100_111), // 0.234375 + // m*2^-2 + Self(0b0_0101_000), // 0.25 + Self(0b0_0101_001), + Self(0b0_0101_010), + Self(0b0_0101_011), + Self(0b0_0101_100), + Self(0b0_0101_101), + Self(0b0_0101_110), + Self(0b0_0101_111), // 0.46875 + // m*2^-1 + Self(0b0_0110_000), // 0.5 + Self(0b0_0110_001), + Self(0b0_0110_010), + Self(0b0_0110_011), + Self(0b0_0110_100), + Self(0b0_0110_101), + Self(0b0_0110_110), + Self(0b0_0110_111), // 0.9375 + // m*2^0 + Self(0b0_0111_000), // 1 + Self(0b0_0111_001), + Self(0b0_0111_010), + Self(0b0_0111_011), + Self(0b0_0111_100), + Self(0b0_0111_101), + Self(0b0_0111_110), + Self(0b0_0111_111), // 1.875 + // m*2^1 + Self(0b0_1000_000), // 2 + Self(0b0_1000_001), + Self(0b0_1000_010), + Self(0b0_1000_011), + Self(0b0_1000_100), + Self(0b0_1000_101), + Self(0b0_1000_110), + Self(0b0_1000_111), // 3.75 + // m*2^2 + Self(0b0_1001_000), // 4 + Self(0b0_1001_001), + Self(0b0_1001_010), + Self(0b0_1001_011), + Self(0b0_1001_100), + Self(0b0_1001_101), + Self(0b0_1001_110), + Self(0b0_1001_111), // 7.5 + // m*2^3 + Self(0b0_1010_000), // 8 + Self(0b0_1010_001), + Self(0b0_1010_010), + Self(0b0_1010_011), + Self(0b0_1010_100), + Self(0b0_1010_101), + Self(0b0_1010_110), + Self(0b0_1010_111), // 15 + // m*2^4 + Self(0b0_1011_000), // 16 + Self(0b0_1011_001), + Self(0b0_1011_010), + Self(0b0_1011_011), + Self(0b0_1011_100), + Self(0b0_1011_101), + Self(0b0_1011_110), + Self(0b0_1011_111), // 30 + // m*2^5 + Self(0b0_1100_000), // 32 + Self(0b0_1100_001), + Self(0b0_1100_010), + Self(0b0_1100_011), + Self(0b0_1100_100), + Self(0b0_1100_101), + Self(0b0_1100_110), + Self(0b0_1100_111), // 60 + // m*2^6 + Self(0b0_1101_000), // 64 + Self(0b0_1101_001), + Self(0b0_1101_010), + Self(0b0_1101_011), + Self(0b0_1101_100), + Self(0b0_1101_101), + Self(0b0_1101_110), + Self(0b0_1101_111), // 120 + // m*2^7 + Self(0b0_1110_000), // 128 + Self(0b0_1110_001), + Self(0b0_1110_010), + Self(0b0_1110_011), + Self(0b0_1110_100), + Self(0b0_1110_101), + Self(0b0_1110_110), + Self(0b0_1110_111), // 240 + ]; +} + +impl ops::Add for f8 { + type Output = Self; + fn add(self, _rhs: Self) -> Self::Output { + unimplemented!() + } +} + +impl ops::Sub for f8 { + type Output = Self; + fn sub(self, _rhs: Self) -> Self::Output { + unimplemented!() + } +} +impl ops::Mul for f8 { + type Output = Self; + fn mul(self, _rhs: Self) -> Self::Output { + unimplemented!() + } +} +impl ops::Div for f8 { + type Output = Self; + fn div(self, _rhs: Self) -> Self::Output { + unimplemented!() + } +} + +impl ops::Neg for f8 { + type Output = Self; + fn neg(self) -> Self::Output { + Self(self.0 ^ Self::SIGN_MASK) + } +} + +impl ops::Rem for f8 { + type Output = Self; + fn rem(self, _rhs: Self) -> Self::Output { + unimplemented!() + } +} + +impl ops::AddAssign for f8 { + fn add_assign(&mut self, _rhs: Self) { + unimplemented!() + } +} + +impl ops::SubAssign for f8 { + fn sub_assign(&mut self, _rhs: Self) { + unimplemented!() + } +} + +impl ops::MulAssign for f8 { + fn mul_assign(&mut self, _rhs: Self) { + unimplemented!() + } +} + +impl cmp::PartialEq for f8 { + fn eq(&self, other: &Self) -> bool { + if self.is_nan() || other.is_nan() { + false + } else if self.abs().to_bits() | other.abs().to_bits() == 0 { + true + } else { + self.0 == other.0 + } + } +} +impl cmp::PartialOrd for f8 { + fn partial_cmp(&self, other: &Self) -> Option { + let inf_rep = f8::EXP_MASK; + + let a_abs = self.abs().to_bits(); + let b_abs = other.abs().to_bits(); + + // If either a or b is NaN, they are unordered. + if a_abs > inf_rep || b_abs > inf_rep { + return None; + } + + // If a and b are both zeros, they are equal. + if a_abs | b_abs == 0 { + return Some(Ordering::Equal); + } + + let a_srep = self.to_bits_signed(); + let b_srep = other.to_bits_signed(); + let res = a_srep.cmp(&b_srep); + + if a_srep & b_srep >= 0 { + // If at least one of a and b is positive, we get the same result comparing + // a and b as signed integers as we would with a fp_ting-point compare. + Some(res) + } else { + // Otherwise, both are negative, so we need to flip the sense of the + // comparison to get the correct result. + Some(res.reverse()) + } + } +} +impl fmt::Display for f8 { + fn fmt(&self, _f: &mut fmt::Formatter<'_>) -> fmt::Result { + unimplemented!() + } +} + +impl fmt::Debug for f8 { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + fmt::Binary::fmt(self, f) + } +} + +impl fmt::Binary for f8 { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + let v = self.0; + write!( + f, + "0b{:b}_{:04b}_{:03b}", + v >> 7, + (v & Self::EXP_MASK) >> Self::SIG_BITS, + v & Self::SIG_MASK + ) + } +} + +impl fmt::LowerHex for f8 { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + self.0.fmt(f) + } +} diff --git a/crates/libm-test/src/lib.rs b/crates/libm-test/src/lib.rs index 17a06b3be..ed7131713 100644 --- a/crates/libm-test/src/lib.rs +++ b/crates/libm-test/src/lib.rs @@ -1,3 +1,6 @@ +#![allow(clippy::unusual_byte_groupings)] // sometimes we group by sign_exp_sig + +mod f8_impl; pub mod gen; #[cfg(feature = "test-multiprecision")] pub mod mpfloat; @@ -5,6 +8,7 @@ pub mod op; mod precision; mod test_traits; +pub use f8_impl::f8; pub use libm::support::{Float, Int, IntTy}; pub use op::{BaseName, Identifier, MathOp, OpCFn, OpFTy, OpRustFn, OpRustRet}; pub use precision::{MaybeOverride, SpecialCase, default_ulp}; From 842b08b4bd78ad04329c17298703203219654f3f Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Thu, 19 Dec 2024 15:10:35 +0000 Subject: [PATCH 5/9] Introduce a float extension trait and some numerical routines --- crates/libm-test/src/lib.rs | 4 +- crates/libm-test/src/num.rs | 458 ++++++++++++++++++++++++++++++++++++ 2 files changed, 461 insertions(+), 1 deletion(-) create mode 100644 crates/libm-test/src/num.rs diff --git a/crates/libm-test/src/lib.rs b/crates/libm-test/src/lib.rs index ed7131713..48b382d20 100644 --- a/crates/libm-test/src/lib.rs +++ b/crates/libm-test/src/lib.rs @@ -4,12 +4,14 @@ mod f8_impl; pub mod gen; #[cfg(feature = "test-multiprecision")] pub mod mpfloat; +mod num; pub mod op; mod precision; mod test_traits; pub use f8_impl::f8; -pub use libm::support::{Float, Int, IntTy}; +pub use libm::support::{Float, Int, IntTy, MinInt}; +pub use num::{FloatExt, logspace}; pub use op::{BaseName, Identifier, MathOp, OpCFn, OpFTy, OpRustFn, OpRustRet}; pub use precision::{MaybeOverride, SpecialCase, default_ulp}; pub use test_traits::{CheckBasis, CheckCtx, CheckOutput, GenerateInput, Hex, TupleCall}; diff --git a/crates/libm-test/src/num.rs b/crates/libm-test/src/num.rs new file mode 100644 index 000000000..4aa7f61b0 --- /dev/null +++ b/crates/libm-test/src/num.rs @@ -0,0 +1,458 @@ +//! Helpful numeric operations. + +use std::cmp::min; + +use libm::support::{CastInto, Float}; + +use crate::{Int, MinInt}; + +/// Extension to `libm`'s `Float` trait with methods that are useful for tests but not +/// needed in `libm` itself. +pub trait FloatExt: Float { + /// The minimum subnormal number. + const TINY_BITS: Self::Int = Self::Int::ONE; + + /// Retrieve additional constants for this float type. + fn consts() -> Consts { + Consts::new() + } + + /// Increment by one ULP, saturating at infinity. + fn next_up(self) -> Self { + let bits = self.to_bits(); + if self.is_nan() || bits == Self::INFINITY.to_bits() { + return self; + } + + let abs = self.abs().to_bits(); + let next_bits = if abs == Self::Int::ZERO { + // Next up from 0 is the smallest subnormal + Self::TINY_BITS + } else if bits == abs { + // Positive: counting up is more positive + bits + Self::Int::ONE + } else { + // Negative: counting down is more positive + bits - Self::Int::ONE + }; + Self::from_bits(next_bits) + } + + /// A faster way to effectively call `next_up` `n` times. + fn n_up(self, n: Self::Int) -> Self { + let bits = self.to_bits(); + if self.is_nan() || bits == Self::INFINITY.to_bits() || n == Self::Int::ZERO { + return self; + } + + let abs = self.abs().to_bits(); + let is_positive = bits == abs; + let crosses_zero = !is_positive && n > abs; + let inf_bits = Self::INFINITY.to_bits(); + + let next_bits = if abs == Self::Int::ZERO { + min(n, inf_bits) + } else if crosses_zero { + min(n - abs, inf_bits) + } else if is_positive { + // Positive, counting up is more positive but this may overflow + match bits.checked_add(n) { + Some(v) if v >= inf_bits => inf_bits, + Some(v) => v, + None => inf_bits, + } + } else { + // Negative, counting down is more positive + bits - n + }; + Self::from_bits(next_bits) + } + + /// Decrement by one ULP, saturating at negative infinity. + fn next_down(self) -> Self { + let bits = self.to_bits(); + if self.is_nan() || bits == Self::NEG_INFINITY.to_bits() { + return self; + } + + let abs = self.abs().to_bits(); + let next_bits = if abs == Self::Int::ZERO { + // Next up from 0 is the smallest negative subnormal + Self::TINY_BITS | Self::SIGN_MASK + } else if bits == abs { + // Positive: counting down is more negative + bits - Self::Int::ONE + } else { + // Negative: counting up is more negative + bits + Self::Int::ONE + }; + Self::from_bits(next_bits) + } + + /// A faster way to effectively call `next_down` `n` times. + fn n_down(self, n: Self::Int) -> Self { + let bits = self.to_bits(); + if self.is_nan() || bits == Self::NEG_INFINITY.to_bits() || n == Self::Int::ZERO { + return self; + } + + let abs = self.abs().to_bits(); + let is_positive = bits == abs; + let crosses_zero = is_positive && n > abs; + let inf_bits = Self::INFINITY.to_bits(); + let ninf_bits = Self::NEG_INFINITY.to_bits(); + + let next_bits = if abs == Self::Int::ZERO { + min(n, inf_bits) | Self::SIGN_MASK + } else if crosses_zero { + min(n - abs, inf_bits) | Self::SIGN_MASK + } else if is_positive { + // Positive, counting down is more negative + bits - n + } else { + // Negative, counting up is more negative but this may overflow + match bits.checked_add(n) { + Some(v) if v > ninf_bits => ninf_bits, + Some(v) => v, + None => ninf_bits, + } + }; + Self::from_bits(next_bits) + } +} + +impl FloatExt for F where F: Float {} + +/// Extra constants that are useful for tests. +#[derive(Debug, Clone, Copy)] +pub struct Consts { + /// The default quiet NaN, which is also the minimum quiet NaN. + pub pos_nan: F, + /// The default quiet NaN with negative sign. + pub neg_nan: F, + /// NaN with maximum (unsigned) significand to be a quiet NaN. The significand is saturated. + pub max_qnan: F, + /// NaN with minimum (unsigned) significand to be a signaling NaN. + pub min_snan: F, + /// NaN with maximum (unsigned) significand to be a signaling NaN. + pub max_snan: F, + pub neg_max_qnan: F, + pub neg_min_snan: F, + pub neg_max_snan: F, +} + +impl Consts { + fn new() -> Self { + let top_sigbit_mask = F::Int::ONE << (F::SIG_BITS - 1); + let pos_nan = F::EXP_MASK | top_sigbit_mask; + let max_qnan = F::EXP_MASK | F::SIG_MASK; + let min_snan = F::EXP_MASK | F::Int::ONE; + let max_snan = (F::EXP_MASK | F::SIG_MASK) ^ top_sigbit_mask; + + let neg_nan = pos_nan | F::SIGN_MASK; + let neg_max_qnan = max_qnan | F::SIGN_MASK; + let neg_min_snan = min_snan | F::SIGN_MASK; + let neg_max_snan = max_snan | F::SIGN_MASK; + + Self { + pos_nan: F::from_bits(pos_nan), + neg_nan: F::from_bits(neg_nan), + max_qnan: F::from_bits(max_qnan), + min_snan: F::from_bits(min_snan), + max_snan: F::from_bits(max_snan), + neg_max_qnan: F::from_bits(neg_max_qnan), + neg_min_snan: F::from_bits(neg_min_snan), + neg_max_snan: F::from_bits(neg_max_snan), + } + } + + pub fn iter(self) -> impl Iterator { + // Destructure so we get unused warnings if we forget a list entry. + let Self { + pos_nan, + neg_nan, + max_qnan, + min_snan, + max_snan, + neg_max_qnan, + neg_min_snan, + neg_max_snan, + } = self; + + [pos_nan, neg_nan, max_qnan, min_snan, max_snan, neg_max_qnan, neg_min_snan, neg_max_snan] + .into_iter() + } +} + +/// Return the number of steps between two floats, returning `None` if either input is NaN. +/// +/// This is the number of steps needed for `n_up` or `n_down` to go between values. Infinities +/// are treated the same as those functions (will return the nearest finite value), and only one +/// of `-0` or `+0` is counted. It does not matter which value is greater. +pub fn ulp_between(x: F, y: F) -> Option { + let a = as_ulp_steps(x)?; + let b = as_ulp_steps(y)?; + Some(a.abs_diff(b)) +} + +/// Return the (signed) number of steps from zero to `x`. +fn as_ulp_steps(x: F) -> Option { + let s = x.to_bits_signed(); + let val = if s >= F::SignedInt::ZERO { + // each increment from `s = 0` is one step up from `x = 0.0` + s + } else { + // each increment from `s = F::SignedInt::MIN` is one step down from `x = -0.0` + F::SignedInt::MIN - s + }; + + // If `x` is NaN, return `None` + (!x.is_nan()).then_some(val) +} + +/// An iterator that returns floats with linearly spaced integer representations, which translates +/// to logarithmic spacing of their values. +/// +/// Note that this tends to skip negative zero, so that needs to be checked explicitly. +pub fn logspace(start: F, end: F, steps: F::Int) -> impl Iterator { + assert!(!start.is_nan()); + assert!(!end.is_nan()); + assert!(end >= start); + + let mut steps = steps.checked_sub(F::Int::ONE).expect("`steps` must be at least 2"); + let between = ulp_between(start, end).expect("`start` or `end` is NaN"); + let spacing = (between / steps).max(F::Int::ONE); + steps = steps.min(between); // At maximum, one step per ULP + + let mut x = start; + (0..=steps.cast()).map(move |_| { + let ret = x; + x = x.n_up(spacing); + ret + }) +} + +#[cfg(test)] +mod tests { + use std::cmp::max; + + use super::*; + use crate::f8; + + #[test] + fn test_next_up_down() { + for (i, v) in f8::ALL.into_iter().enumerate() { + let down = v.next_down().to_bits(); + let up = v.next_up().to_bits(); + + if i == 0 { + assert_eq!(down, f8::NEG_INFINITY.to_bits(), "{i} next_down({v:#010b})"); + } else { + let expected = + if v == f8::ZERO { 1 | f8::SIGN_MASK } else { f8::ALL[i - 1].to_bits() }; + assert_eq!(down, expected, "{i} next_down({v:#010b})"); + } + + if i == f8::ALL_LEN - 1 { + assert_eq!(up, f8::INFINITY.to_bits(), "{i} next_up({v:#010b})"); + } else { + let expected = if v == f8::NEG_ZERO { 1 } else { f8::ALL[i + 1].to_bits() }; + assert_eq!(up, expected, "{i} next_up({v:#010b})"); + } + } + } + + #[test] + fn test_next_up_down_inf_nan() { + assert_eq!(f8::NEG_INFINITY.next_up().to_bits(), f8::ALL[0].to_bits(),); + assert_eq!(f8::NEG_INFINITY.next_down().to_bits(), f8::NEG_INFINITY.to_bits(),); + assert_eq!(f8::INFINITY.next_down().to_bits(), f8::ALL[f8::ALL_LEN - 1].to_bits(),); + assert_eq!(f8::INFINITY.next_up().to_bits(), f8::INFINITY.to_bits(),); + assert_eq!(f8::NAN.next_up().to_bits(), f8::NAN.to_bits(),); + assert_eq!(f8::NAN.next_down().to_bits(), f8::NAN.to_bits(),); + } + + #[test] + fn test_n_up_down_quick() { + assert_eq!(f8::ALL[0].n_up(4).to_bits(), f8::ALL[4].to_bits(),); + assert_eq!( + f8::ALL[f8::ALL_LEN - 1].n_down(4).to_bits(), + f8::ALL[f8::ALL_LEN - 5].to_bits(), + ); + + // Check around zero + assert_eq!(f8::from_bits(0b0).n_up(7).to_bits(), 0b0_0000_111); + assert_eq!(f8::from_bits(0b0).n_down(7).to_bits(), 0b1_0000_111); + + // Check across zero + assert_eq!(f8::from_bits(0b1_0000_111).n_up(8).to_bits(), 0b0_0000_001); + assert_eq!(f8::from_bits(0b0_0000_111).n_down(8).to_bits(), 0b1_0000_001); + } + + #[test] + fn test_n_up_down_one() { + // Verify that `n_up(1)` and `n_down(1)` are the same as `next_up()` and next_down()`.` + for i in 0..u8::MAX { + let v = f8::from_bits(i); + assert_eq!(v.next_up().to_bits(), v.n_up(1).to_bits()); + assert_eq!(v.next_down().to_bits(), v.n_down(1).to_bits()); + } + } + + #[test] + fn test_n_up_down_inf_nan_zero() { + assert_eq!(f8::NEG_INFINITY.n_up(1).to_bits(), f8::ALL[0].to_bits()); + assert_eq!(f8::NEG_INFINITY.n_up(239).to_bits(), f8::ALL[f8::ALL_LEN - 1].to_bits()); + assert_eq!(f8::NEG_INFINITY.n_up(240).to_bits(), f8::INFINITY.to_bits()); + assert_eq!(f8::NEG_INFINITY.n_down(u8::MAX).to_bits(), f8::NEG_INFINITY.to_bits()); + + assert_eq!(f8::INFINITY.n_down(1).to_bits(), f8::ALL[f8::ALL_LEN - 1].to_bits()); + assert_eq!(f8::INFINITY.n_down(239).to_bits(), f8::ALL[0].to_bits()); + assert_eq!(f8::INFINITY.n_down(240).to_bits(), f8::NEG_INFINITY.to_bits()); + assert_eq!(f8::INFINITY.n_up(u8::MAX).to_bits(), f8::INFINITY.to_bits()); + + assert_eq!(f8::NAN.n_up(u8::MAX).to_bits(), f8::NAN.to_bits()); + assert_eq!(f8::NAN.n_down(u8::MAX).to_bits(), f8::NAN.to_bits()); + + assert_eq!(f8::ZERO.n_down(1).to_bits(), f8::TINY_BITS | f8::SIGN_MASK); + assert_eq!(f8::NEG_ZERO.n_up(1).to_bits(), f8::TINY_BITS); + } + + /// True if the specified range of `f8::ALL` includes both +0 and -0 + fn crossed_zero(start: usize, end: usize) -> bool { + let crossed = &f8::ALL[start..=end]; + crossed.iter().any(|f| f8::eq_repr(*f, f8::ZERO)) + && crossed.iter().any(|f| f8::eq_repr(*f, f8::NEG_ZERO)) + } + + #[test] + fn test_n_up_down() { + for (i, v) in f8::ALL.into_iter().enumerate() { + for n in 0..f8::ALL_LEN { + let down = v.n_down(n as u8).to_bits(); + let up = v.n_up(n as u8).to_bits(); + + if let Some(down_exp_idx) = i.checked_sub(n) { + // No overflow + let mut expected = f8::ALL[down_exp_idx].to_bits(); + if n >= 1 && crossed_zero(down_exp_idx, i) { + // If both -0 and +0 are included, we need to adjust our expected value + match down_exp_idx.checked_sub(1) { + Some(v) => expected = f8::ALL[v].to_bits(), + // Saturate to -inf if we are out of values + None => expected = f8::NEG_INFINITY.to_bits(), + } + } + assert_eq!(down, expected, "{i} {n} n_down({v:#010b})"); + } else { + // Overflow to -inf + assert_eq!(down, f8::NEG_INFINITY.to_bits(), "{i} {n} n_down({v:#010b})"); + } + + let mut up_exp_idx = i + n; + if up_exp_idx < f8::ALL_LEN { + // No overflow + if n >= 1 && up_exp_idx < f8::ALL_LEN && crossed_zero(i, up_exp_idx) { + // If both -0 and +0 are included, we need to adjust our expected value + up_exp_idx += 1; + } + + let expected = if up_exp_idx >= f8::ALL_LEN { + f8::INFINITY.to_bits() + } else { + f8::ALL[up_exp_idx].to_bits() + }; + + assert_eq!(up, expected, "{i} {n} n_up({v:#010b})"); + } else { + // Overflow to +inf + assert_eq!(up, f8::INFINITY.to_bits(), "{i} {n} n_up({v:#010b})"); + } + } + } + } + + #[test] + fn test_ulp_between() { + for (i, x) in f8::ALL.into_iter().enumerate() { + for (j, y) in f8::ALL.into_iter().enumerate() { + let ulp = ulp_between(x, y).unwrap(); + let make_msg = || format!("i: {i} j: {j} x: {x:b} y: {y:b} ulp {ulp}"); + + let i_low = min(i, j); + let i_hi = max(i, j); + let mut expected = u8::try_from(i_hi - i_low).unwrap(); + if crossed_zero(i_low, i_hi) { + expected -= 1; + } + + assert_eq!(ulp, expected, "{}", make_msg()); + + // Skip if either are zero since `next_{up,down}` will count over it + let either_zero = x == f8::ZERO || y == f8::ZERO; + if x < y && !either_zero { + assert_eq!(x.n_up(ulp).to_bits(), y.to_bits(), "{}", make_msg()); + assert_eq!(y.n_down(ulp).to_bits(), x.to_bits(), "{}", make_msg()); + } else if !either_zero { + assert_eq!(y.n_up(ulp).to_bits(), x.to_bits(), "{}", make_msg()); + assert_eq!(x.n_down(ulp).to_bits(), y.to_bits(), "{}", make_msg()); + } + } + } + } + + #[test] + fn test_ulp_between_inf_nan_zero() { + assert_eq!(ulp_between(f8::NEG_INFINITY, f8::INFINITY).unwrap(), f8::ALL_LEN as u8); + assert_eq!(ulp_between(f8::INFINITY, f8::NEG_INFINITY).unwrap(), f8::ALL_LEN as u8); + assert_eq!( + ulp_between(f8::NEG_INFINITY, f8::ALL[f8::ALL_LEN - 1]).unwrap(), + f8::ALL_LEN as u8 - 1 + ); + assert_eq!(ulp_between(f8::INFINITY, f8::ALL[0]).unwrap(), f8::ALL_LEN as u8 - 1); + + assert_eq!(ulp_between(f8::ZERO, f8::NEG_ZERO).unwrap(), 0); + assert_eq!(ulp_between(f8::NAN, f8::ZERO), None); + assert_eq!(ulp_between(f8::ZERO, f8::NAN), None); + } + + #[test] + fn test_logspace() { + let ls: Vec<_> = logspace(f8::from_bits(0x0), f8::from_bits(0x4), 2).collect(); + let exp = [f8::from_bits(0x0), f8::from_bits(0x4)]; + assert_eq!(ls, exp); + + let ls: Vec<_> = logspace(f8::from_bits(0x0), f8::from_bits(0x4), 3).collect(); + let exp = [f8::from_bits(0x0), f8::from_bits(0x2), f8::from_bits(0x4)]; + assert_eq!(ls, exp); + + // Check that we include all values with no repeats if `steps` exceeds the maximum number + // of steps. + let ls: Vec<_> = logspace(f8::from_bits(0x0), f8::from_bits(0x3), 10).collect(); + let exp = [f8::from_bits(0x0), f8::from_bits(0x1), f8::from_bits(0x2), f8::from_bits(0x3)]; + assert_eq!(ls, exp); + } + + #[test] + fn test_consts() { + let Consts { + pos_nan, + neg_nan, + max_qnan, + min_snan, + max_snan, + neg_max_qnan, + neg_min_snan, + neg_max_snan, + } = f8::consts(); + + assert_eq!(pos_nan.to_bits(), 0b0_1111_100); + assert_eq!(neg_nan.to_bits(), 0b1_1111_100); + assert_eq!(max_qnan.to_bits(), 0b0_1111_111); + assert_eq!(min_snan.to_bits(), 0b0_1111_001); + assert_eq!(max_snan.to_bits(), 0b0_1111_011); + assert_eq!(neg_max_qnan.to_bits(), 0b1_1111_111); + assert_eq!(neg_min_snan.to_bits(), 0b1_1111_001); + assert_eq!(neg_max_snan.to_bits(), 0b1_1111_011); + } +} From a114878d53f470b5a76f8905318392e18f617c01 Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Thu, 19 Dec 2024 11:19:01 +0000 Subject: [PATCH 6/9] Add interfaces and tests based on function domains Create a type representing a function's domain and a test that does a logarithmic sweep of points within the domain. --- crates/libm-test/src/domain.rs | 186 ++++++++++++++++++++ crates/libm-test/src/gen.rs | 1 + crates/libm-test/src/gen/domain_logspace.rs | 43 +++++ crates/libm-test/src/lib.rs | 1 + crates/libm-test/tests/multiprecision.rs | 101 ++++++++++- 5 files changed, 327 insertions(+), 5 deletions(-) create mode 100644 crates/libm-test/src/domain.rs create mode 100644 crates/libm-test/src/gen/domain_logspace.rs diff --git a/crates/libm-test/src/domain.rs b/crates/libm-test/src/domain.rs new file mode 100644 index 000000000..43ba21974 --- /dev/null +++ b/crates/libm-test/src/domain.rs @@ -0,0 +1,186 @@ +//! Traits and operations related to bounds of a function. + +use std::fmt; +use std::ops::{self, Bound}; + +use crate::Float; + +/// Representation of a function's domain. +#[derive(Clone, Debug)] +pub struct Domain { + /// Start of the region for which a function is defined (ignoring poles). + pub start: Bound, + /// Endof the region for which a function is defined (ignoring poles). + pub end: Bound, + /// Additional points to check closer around. These can be e.g. undefined asymptotes or + /// inflection points. + pub check_points: Option BoxIter>, +} + +type BoxIter = Box>; + +impl Domain { + /// The start of this domain, saturating at negative infinity. + pub fn range_start(&self) -> F { + match self.start { + Bound::Included(v) => v, + Bound::Excluded(v) => v.next_up(), + Bound::Unbounded => F::NEG_INFINITY, + } + } + + /// The end of this domain, saturating at infinity. + pub fn range_end(&self) -> F { + match self.end { + Bound::Included(v) => v, + Bound::Excluded(v) => v.next_down(), + Bound::Unbounded => F::INFINITY, + } + } +} + +impl Domain { + /// x ∈ ℝ + pub const UNBOUNDED: Self = + Self { start: Bound::Unbounded, end: Bound::Unbounded, check_points: None }; + + /// x ∈ ℝ >= 0 + pub const POSITIVE: Self = + Self { start: Bound::Included(F::ZERO), end: Bound::Unbounded, check_points: None }; + + /// x ∈ ℝ > 0 + pub const STRICTLY_POSITIVE: Self = + Self { start: Bound::Excluded(F::ZERO), end: Bound::Unbounded, check_points: None }; + + /// Used for versions of `asin` and `acos`. + pub const INVERSE_TRIG_PERIODIC: Self = Self { + start: Bound::Included(F::NEG_ONE), + end: Bound::Included(F::ONE), + check_points: None, + }; + + /// Domain for `acosh` + pub const ACOSH: Self = + Self { start: Bound::Included(F::ONE), end: Bound::Unbounded, check_points: None }; + + /// Domain for `atanh` + pub const ATANH: Self = Self { + start: Bound::Excluded(F::NEG_ONE), + end: Bound::Excluded(F::ONE), + check_points: None, + }; + + /// Domain for `sin`, `cos`, and `tan` + pub const TRIG: Self = Self { + // TODO + check_points: Some(|| Box::new([-F::PI, -F::FRAC_PI_2, F::FRAC_PI_2, F::PI].into_iter())), + ..Self::UNBOUNDED + }; + + /// Domain for `log` in various bases + pub const LOG: Self = Self::STRICTLY_POSITIVE; + + /// Domain for `log1p` i.e. `log(1 + x)` + pub const LOG1P: Self = + Self { start: Bound::Excluded(F::NEG_ONE), end: Bound::Unbounded, check_points: None }; + + /// Domain for `sqrt` + pub const SQRT: Self = Self::POSITIVE; + + /// Domain for `gamma` + pub const GAMMA: Self = Self { + check_points: Some(|| { + // Negative integers are asymptotes + Box::new((0..u8::MAX).map(|scale| { + let mut base = F::ZERO; + for _ in 0..scale { + base = base - F::ONE; + } + base + })) + }), + // Whether or not gamma is defined for negative numbers is implementation dependent + ..Self::UNBOUNDED + }; + + /// Domain for `loggamma` + pub const LGAMMA: Self = Self::STRICTLY_POSITIVE; +} + +/// Implement on `op::*` types to indicate how they are bounded. +pub trait HasDomain +where + T: Copy + fmt::Debug + ops::Add + ops::Sub + PartialOrd + 'static, +{ + const DOMAIN: Domain; +} + +/// Implement [`HasDomain`] for both the `f32` and `f64` variants of a function. +macro_rules! impl_has_domain { + ($($fn_name:ident => $domain:expr;)*) => { + paste::paste! { + $( + // Implement for f64 functions + impl HasDomain for $crate::op::$fn_name::Routine { + const DOMAIN: Domain = Domain::::$domain; + } + + // Implement for f32 functions + impl HasDomain for $crate::op::[< $fn_name f >]::Routine { + const DOMAIN: Domain = Domain::::$domain; + } + )* + } + }; +} + +// Tie functions together with their domains. +impl_has_domain! { + acos => INVERSE_TRIG_PERIODIC; + acosh => ACOSH; + asin => INVERSE_TRIG_PERIODIC; + asinh => UNBOUNDED; + atan => UNBOUNDED; + atanh => ATANH; + cbrt => UNBOUNDED; + ceil => UNBOUNDED; + cos => TRIG; + cosh => UNBOUNDED; + erf => UNBOUNDED; + exp => UNBOUNDED; + exp10 => UNBOUNDED; + exp2 => UNBOUNDED; + expm1 => UNBOUNDED; + fabs => UNBOUNDED; + floor => UNBOUNDED; + frexp => UNBOUNDED; + ilogb => UNBOUNDED; + j0 => UNBOUNDED; + j1 => UNBOUNDED; + lgamma => LGAMMA; + log => LOG; + log10 => LOG; + log1p => LOG1P; + log2 => LOG; + modf => UNBOUNDED; + rint => UNBOUNDED; + round => UNBOUNDED; + sin => TRIG; + sincos => TRIG; + sinh => UNBOUNDED; + sqrt => SQRT; + tan => TRIG; + tanh => UNBOUNDED; + tgamma => GAMMA; + trunc => UNBOUNDED; +} + +/* Manual implementations, these functions don't follow `foo`->`foof` naming */ + +impl HasDomain for crate::op::lgammaf_r::Routine { + const DOMAIN: Domain = Domain::::LGAMMA; +} + +impl HasDomain for crate::op::lgamma_r::Routine { + const DOMAIN: Domain = Domain::::LGAMMA; +} diff --git a/crates/libm-test/src/gen.rs b/crates/libm-test/src/gen.rs index 3e9eca37a..e3c88c44a 100644 --- a/crates/libm-test/src/gen.rs +++ b/crates/libm-test/src/gen.rs @@ -1,6 +1,7 @@ //! Different generators that can create random or systematic bit patterns. use crate::GenerateInput; +pub mod domain_logspace; pub mod random; /// Helper type to turn any reusable input into a generator. diff --git a/crates/libm-test/src/gen/domain_logspace.rs b/crates/libm-test/src/gen/domain_logspace.rs new file mode 100644 index 000000000..e8cdb9d2b --- /dev/null +++ b/crates/libm-test/src/gen/domain_logspace.rs @@ -0,0 +1,43 @@ +//! A generator that produces logarithmically spaced values within domain bounds. + +use libm::support::{IntTy, MinInt}; + +use crate::domain::HasDomain; +use crate::op::OpITy; +use crate::{MathOp, logspace}; + +/// Number of tests to run. +// FIXME(ntests): replace this with a more logical algorithm +const NTESTS: usize = { + if cfg!(optimizations_enabled) { + if crate::emulated() + || !cfg!(target_pointer_width = "64") + || cfg!(all(target_arch = "x86_64", target_vendor = "apple")) + { + // Tests are pretty slow on non-64-bit targets, x86 MacOS, and targets that run + // in QEMU. + 100_000 + } else { + 5_000_000 + } + } else { + // Without optimizations just run a quick check + 800 + } +}; + +/// Create a range of logarithmically spaced inputs within a function's domain. +/// +/// This allows us to get reasonably thorough coverage without wasting time on values that are +/// NaN or out of range. Random tests will still cover values that are excluded here. +pub fn get_test_cases() -> impl Iterator +where + Op: MathOp + HasDomain, + IntTy: TryFrom, +{ + let domain = Op::DOMAIN; + let start = domain.range_start(); + let end = domain.range_end(); + let steps = OpITy::::try_from(NTESTS).unwrap_or(OpITy::::MAX); + logspace(start, end, steps).map(|v| (v,)) +} diff --git a/crates/libm-test/src/lib.rs b/crates/libm-test/src/lib.rs index 48b382d20..622b2dec9 100644 --- a/crates/libm-test/src/lib.rs +++ b/crates/libm-test/src/lib.rs @@ -1,5 +1,6 @@ #![allow(clippy::unusual_byte_groupings)] // sometimes we group by sign_exp_sig +pub mod domain; mod f8_impl; pub mod gen; #[cfg(feature = "test-multiprecision")] diff --git a/crates/libm-test/tests/multiprecision.rs b/crates/libm-test/tests/multiprecision.rs index 0b41fba82..e643f3c9c 100644 --- a/crates/libm-test/tests/multiprecision.rs +++ b/crates/libm-test/tests/multiprecision.rs @@ -2,11 +2,14 @@ #![cfg(feature = "test-multiprecision")] -use libm_test::gen::{CachedInput, random}; +use libm_test::domain::HasDomain; +use libm_test::gen::{CachedInput, domain_logspace, random}; use libm_test::mpfloat::MpOp; -use libm_test::{CheckBasis, CheckCtx, CheckOutput, GenerateInput, MathOp, TupleCall}; +use libm_test::{ + CheckBasis, CheckCtx, CheckOutput, GenerateInput, MathOp, OpFTy, OpRustFn, OpRustRet, TupleCall, +}; -/// Implement a test against MPFR with random inputs. +/// Test against MPFR with random inputs. macro_rules! mp_rand_tests { ( fn_name: $fn_name:ident, @@ -16,13 +19,14 @@ macro_rules! mp_rand_tests { #[test] $(#[$meta])* fn [< mp_random_ $fn_name >]() { - test_one::(); + test_one_random::(); } } }; } -fn test_one() +/// Test a single routine with random inputs +fn test_one_random() where Op: MathOp + MpOp, CachedInput: GenerateInput, @@ -67,3 +71,90 @@ libm_macros::for_each_function! { nextafterf, ], } + +/// Test against MPFR with generators from a domain. +macro_rules! mp_domain_tests { + ( + fn_name: $fn_name:ident, + attrs: [$($meta:meta)*] + ) => { + paste::paste! { + #[test] + $(#[$meta])* + fn [< mp_logspace_ $fn_name >]() { + type Op = libm_test::op::$fn_name::Routine; + domain_test_runner::(domain_logspace::get_test_cases::()); + } + } + }; +} + +/// Test a single routine against domaine-aware inputs. +fn domain_test_runner(cases: impl Iterator) +where + // Complicated generics... + // The operation must take a single float argument (unary only) + Op: MathOp::FTy,)>, + // It must also support multiprecision operations + Op: MpOp, + // And it must have a domain specified + Op: HasDomain, + // The single float argument tuple must be able to call the `RustFn` and return `RustRet` + (OpFTy,): TupleCall, Output = OpRustRet>, +{ + let mut mp_vals = Op::new_mp(); + let ctx = CheckCtx::new(Op::IDENTIFIER, CheckBasis::Mpfr); + + for input in cases { + let mp_res = Op::run(&mut mp_vals, input); + let crate_res = input.call(Op::ROUTINE); + + crate_res.validate(mp_res, input, &ctx).unwrap(); + } +} + +libm_macros::for_each_function! { + callback: mp_domain_tests, + attributes: [], + skip: [ + // Functions with multiple inputs + atan2, + atan2f, + copysign, + copysignf, + fdim, + fdimf, + fma, + fmaf, + fmax, + fmaxf, + fmin, + fminf, + fmod, + fmodf, + hypot, + hypotf, + jn, + jnf, + ldexp, + ldexpf, + nextafter, + nextafterf, + pow, + powf, + remainder, + remainderf, + remquo, + remquof, + scalbn, + scalbnf, + + // FIXME: MPFR tests needed + frexp, + frexpf, + ilogb, + ilogbf, + modf, + modff, + ], +} From e1935690fe724669209ce286a29a0e74bfd59a54 Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Thu, 19 Dec 2024 11:22:02 +0000 Subject: [PATCH 7/9] Add tests for edge cases Introduce a generator that will tests various points of interest including zeros, infinities, and NaNs. --- crates/libm-test/src/domain.rs | 4 +- crates/libm-test/src/gen.rs | 1 + crates/libm-test/src/gen/edge_cases.rs | 90 ++++++++++++++++++++++++ crates/libm-test/tests/multiprecision.rs | 9 ++- 4 files changed, 101 insertions(+), 3 deletions(-) create mode 100644 crates/libm-test/src/gen/edge_cases.rs diff --git a/crates/libm-test/src/domain.rs b/crates/libm-test/src/domain.rs index 43ba21974..9ee8a19b9 100644 --- a/crates/libm-test/src/domain.rs +++ b/crates/libm-test/src/domain.rs @@ -3,7 +3,7 @@ use std::fmt; use std::ops::{self, Bound}; -use crate::Float; +use crate::{Float, FloatExt}; /// Representation of a function's domain. #[derive(Clone, Debug)] @@ -19,7 +19,7 @@ pub struct Domain { type BoxIter = Box>; -impl Domain { +impl Domain { /// The start of this domain, saturating at negative infinity. pub fn range_start(&self) -> F { match self.start { diff --git a/crates/libm-test/src/gen.rs b/crates/libm-test/src/gen.rs index e3c88c44a..2d15915d9 100644 --- a/crates/libm-test/src/gen.rs +++ b/crates/libm-test/src/gen.rs @@ -2,6 +2,7 @@ use crate::GenerateInput; pub mod domain_logspace; +pub mod edge_cases; pub mod random; /// Helper type to turn any reusable input into a generator. diff --git a/crates/libm-test/src/gen/edge_cases.rs b/crates/libm-test/src/gen/edge_cases.rs new file mode 100644 index 000000000..625e18bc7 --- /dev/null +++ b/crates/libm-test/src/gen/edge_cases.rs @@ -0,0 +1,90 @@ +//! A generator that checks a handful of cases near infinities, zeros, asymptotes, and NaNs. + +use libm::support::Float; + +use crate::domain::HasDomain; +use crate::{FloatExt, MathOp}; + +/// Number of values near an interesting point to check. +// FIXME(ntests): replace this with a more logical algorithm +const AROUND: usize = 100; + +/// Functions have infinite asymptotes, limit how many we check. +// FIXME(ntests): replace this with a more logical algorithm +const MAX_CHECK_POINTS: usize = 10; + +/// Create a list of values around interesting points (infinities, zeroes, NaNs). +pub fn get_test_cases() -> impl Iterator +where + Op: MathOp + HasDomain, + F: Float, +{ + let mut ret = Vec::new(); + let values = &mut ret; + let domain = Op::DOMAIN; + let domain_start = domain.range_start(); + let domain_end = domain.range_end(); + + // Check near some notable constants + count_up(F::ONE, values); + count_up(F::ZERO, values); + count_up(F::NEG_ONE, values); + count_down(F::ONE, values); + count_down(F::ZERO, values); + count_down(F::NEG_ONE, values); + values.push(F::NEG_ZERO); + + // Check values near the extremes + count_up(F::NEG_INFINITY, values); + count_down(F::INFINITY, values); + count_down(domain_end, values); + count_up(domain_start, values); + count_down(domain_start, values); + count_up(domain_end, values); + count_down(domain_end, values); + + // Check some special values that aren't included in the above ranges + values.push(F::NAN); + values.extend(F::consts().iter()); + + // Check around asymptotes + if let Some(f) = domain.check_points { + let iter = f(); + for x in iter.take(MAX_CHECK_POINTS) { + count_up(x, values); + count_down(x, values); + } + } + + // Some results may overlap so deduplicate the vector to save test cycles. + values.sort_by_key(|x| x.to_bits()); + values.dedup_by_key(|x| x.to_bits()); + + ret.into_iter().map(|v| (v,)) +} + +/// Add `AROUND` values starting at and including `x` and counting up. Uses the smallest possible +/// increments (1 ULP). +fn count_up(mut x: F, values: &mut Vec) { + assert!(!x.is_nan()); + + let mut count = 0; + while x < F::INFINITY && count < AROUND { + values.push(x); + x = x.next_up(); + count += 1; + } +} + +/// Add `AROUND` values starting at and including `x` and counting down. Uses the smallest possible +/// increments (1 ULP). +fn count_down(mut x: F, values: &mut Vec) { + assert!(!x.is_nan()); + + let mut count = 0; + while x > F::NEG_INFINITY && count < AROUND { + values.push(x); + x = x.next_down(); + count += 1; + } +} diff --git a/crates/libm-test/tests/multiprecision.rs b/crates/libm-test/tests/multiprecision.rs index e643f3c9c..5255dc1cf 100644 --- a/crates/libm-test/tests/multiprecision.rs +++ b/crates/libm-test/tests/multiprecision.rs @@ -3,7 +3,7 @@ #![cfg(feature = "test-multiprecision")] use libm_test::domain::HasDomain; -use libm_test::gen::{CachedInput, domain_logspace, random}; +use libm_test::gen::{CachedInput, domain_logspace, edge_cases, random}; use libm_test::mpfloat::MpOp; use libm_test::{ CheckBasis, CheckCtx, CheckOutput, GenerateInput, MathOp, OpFTy, OpRustFn, OpRustRet, TupleCall, @@ -79,6 +79,13 @@ macro_rules! mp_domain_tests { attrs: [$($meta:meta)*] ) => { paste::paste! { + #[test] + $(#[$meta])* + fn [< mp_edge_case_ $fn_name >]() { + type Op = libm_test::op::$fn_name::Routine; + domain_test_runner::(edge_cases::get_test_cases::()); + } + #[test] $(#[$meta])* fn [< mp_logspace_ $fn_name >]() { From 5b26d7901e65cd5ef109fe117305010663f40bc6 Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Thu, 19 Dec 2024 14:23:05 +0000 Subject: [PATCH 8/9] Update allowed precision to account for new tests --- crates/libm-test/src/precision.rs | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/crates/libm-test/src/precision.rs b/crates/libm-test/src/precision.rs index c7f9d9e30..b878212fa 100644 --- a/crates/libm-test/src/precision.rs +++ b/crates/libm-test/src/precision.rs @@ -41,10 +41,11 @@ pub fn default_ulp(ctx: &CheckCtx) -> u32 { (Musl, Id::Tgamma) => 20, // Overrides for MPFR + (Mpfr, Id::Acosh) => 4, (Mpfr, Id::Acoshf) => 4, (Mpfr, Id::Asinh | Id::Asinhf) => 2, (Mpfr, Id::Atanh | Id::Atanhf) => 2, - (Mpfr, Id::Exp10 | Id::Exp10f) => 3, + (Mpfr, Id::Exp10 | Id::Exp10f) => 6, (Mpfr, Id::Lgamma | Id::LgammaR | Id::Lgammaf | Id::LgammafR) => 16, (Mpfr, Id::Sinh | Id::Sinhf) => 2, (Mpfr, Id::Tanh | Id::Tanhf) => 2, @@ -105,17 +106,14 @@ impl MaybeOverride<(f32,)> for SpecialCase { _ulp: &mut u32, ctx: &CheckCtx, ) -> Option { - if ctx.basis == CheckBasis::Musl { - if ctx.base_name == BaseName::Expm1 && input.0 > 80.0 && actual.is_infinite() { - // we return infinity but the number is representable - return XFAIL; - } + if ctx.base_name == BaseName::Expm1 && input.0 > 80.0 && actual.is_infinite() { + // we return infinity but the number is representable + return XFAIL; + } - if ctx.base_name == BaseName::Sinh && input.0.abs() > 80.0 && actual.is_nan() { - // we return some NaN that should be real values or infinite - // doesn't seem to happen on x86 - return XFAIL; - } + if ctx.base_name == BaseName::Sinh && input.0.abs() > 80.0 && actual.is_nan() { + // we return some NaN that should be real values or infinite + return XFAIL; } if ctx.base_name == BaseName::Acosh && input.0 < -1.0 { From f4d97cd3f64d7a6682a2cbde6ed6233cc20c2e0f Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Thu, 19 Dec 2024 11:22:22 +0000 Subject: [PATCH 9/9] Add a way to plot the output from generators For visualization, add a simple script for generating scatter plots and a binary (via examples) to plot the inputs given various domains. --- crates/libm-test/examples/plot_domains.rs | 105 +++++++++++++++ crates/libm-test/examples/plot_file.jl | 157 ++++++++++++++++++++++ 2 files changed, 262 insertions(+) create mode 100644 crates/libm-test/examples/plot_domains.rs create mode 100644 crates/libm-test/examples/plot_file.jl diff --git a/crates/libm-test/examples/plot_domains.rs b/crates/libm-test/examples/plot_domains.rs new file mode 100644 index 000000000..630a0c233 --- /dev/null +++ b/crates/libm-test/examples/plot_domains.rs @@ -0,0 +1,105 @@ +//! Program to write all inputs from a generator to a file, then invoke a Julia script to plot +//! them. Output is in `target/plots`. +//! +//! Requires Julia with the `CairoMakie` dependency. +//! +//! Note that running in release mode by default generates a _lot_ more datapoints, which +//! causes plotting to be extremely slow (some simplification to be done in the script). + +use std::fmt::Write as _; +use std::io::{BufWriter, Write}; +use std::path::Path; +use std::process::Command; +use std::{env, fs}; + +use libm_test::domain::HasDomain; +use libm_test::gen::{domain_logspace, edge_cases}; +use libm_test::{MathOp, op}; + +const JL_PLOT: &str = "examples/plot_file.jl"; + +fn main() { + let manifest_env = env::var("CARGO_MANIFEST_DIR").unwrap(); + let manifest_dir = Path::new(&manifest_env); + let out_dir = manifest_dir.join("../../target/plots"); + if !out_dir.exists() { + fs::create_dir(&out_dir).unwrap(); + } + + let jl_script = manifest_dir.join(JL_PLOT); + let mut config = format!(r#"out_dir = "{}""#, out_dir.display()); + config.write_str("\n\n").unwrap(); + + // Plot a few domains with some functions that use them. + plot_one_operator::(&out_dir, &mut config); + plot_one_operator::(&out_dir, &mut config); + plot_one_operator::(&out_dir, &mut config); + + let config_path = out_dir.join("config.toml"); + fs::write(&config_path, config).unwrap(); + + // The script expects a path to `config.toml` to be passed as its only argument + let mut cmd = Command::new("julia"); + if cfg!(optimizations_enabled) { + cmd.arg("-O3"); + } + cmd.arg(jl_script).arg(config_path); + + println!("launching script... {cmd:?}"); + cmd.status().unwrap(); +} + +/// Run multiple generators for a single operator. +fn plot_one_operator(out_dir: &Path, config: &mut String) +where + Op: MathOp + HasDomain, +{ + plot_one_generator( + out_dir, + Op::BASE_NAME.as_str(), + "logspace", + config, + domain_logspace::get_test_cases::(), + ); + plot_one_generator( + out_dir, + Op::BASE_NAME.as_str(), + "edge_cases", + config, + edge_cases::get_test_cases::(), + ); +} + +/// Plot the output of a single generator. +fn plot_one_generator( + out_dir: &Path, + fn_name: &str, + gen_name: &str, + config: &mut String, + gen: impl Iterator, +) { + let text_file = out_dir.join(format!("input-{fn_name}-{gen_name}.txt")); + + let f = fs::File::create(&text_file).unwrap(); + let mut w = BufWriter::new(f); + let mut count = 0u64; + + for input in gen { + writeln!(w, "{:e}", input.0).unwrap(); + count += 1; + } + + w.flush().unwrap(); + println!("generated {count} inputs for {fn_name}-{gen_name}"); + + writeln!( + config, + r#"[[input]] +function = "{fn_name}" +generator = "{gen_name}" +input_file = "{}" +"#, + text_file.to_str().unwrap() + ) + .unwrap() +} diff --git a/crates/libm-test/examples/plot_file.jl b/crates/libm-test/examples/plot_file.jl new file mode 100644 index 000000000..14a128303 --- /dev/null +++ b/crates/libm-test/examples/plot_file.jl @@ -0,0 +1,157 @@ +"A quick script for plotting a list of floats. + +Takes a path to a TOML file (Julia has builtin TOML support but not JSON) which +specifies a list of source files to plot. Plots are done with both a linear and +a log scale. + +Requires [Makie] (specifically CairoMakie) for plotting. + +[Makie]: https://docs.makie.org/stable/ +" + +using CairoMakie +using TOML + +function main()::Nothing + CairoMakie.activate!(px_per_unit=10) + config_path = ARGS[1] + + cfg = Dict() + open(config_path, "r") do f + cfg = TOML.parse(f) + end + + out_dir = cfg["out_dir"] + for input in cfg["input"] + fn_name = input["function"] + gen_name = input["generator"] + input_file = input["input_file"] + + plot_one(input_file, out_dir, fn_name, gen_name) + end +end + +"Read inputs from a file, create both linear and log plots for one function" +function plot_one( + input_file::String, + out_dir::String, + fn_name::String, + gen_name::String, +)::Nothing + fig = Figure() + + lin_out_file = joinpath(out_dir, "plot-$fn_name-$gen_name.png") + log_out_file = joinpath(out_dir, "plot-$fn_name-$gen_name-log.png") + + # Map string function names to callable functions + if fn_name == "cos" + orig_func = cos + xlims = (-6.0, 6.0) + xlims_log = (-pi * 10, pi * 10) + elseif fn_name == "cbrt" + orig_func = cbrt + xlims = (-2.0, 2.0) + xlims_log = (-1000.0, 1000.0) + elseif fn_name == "sqrt" + orig_func = sqrt + xlims = (-1.1, 6.0) + xlims_log = (-1.1, 5000.0) + else + println("unrecognized function name `$fn_name`; update plot_file.jl") + exit(1) + end + + # Edge cases don't do much beyond +/-1, except for infinity. + if gen_name == "edge_cases" + xlims = (-1.1, 1.1) + xlims_log = (-1.1, 1.1) + end + + # Turn domain errors into NaN + func(x) = map_or(x, orig_func, NaN) + + # Parse a series of X values produced by the generator + inputs = readlines(input_file) + gen_x = map((v) -> parse(Float32, v), inputs) + + do_plot( + fig, gen_x, func, xlims[1], xlims[2], + "$fn_name $gen_name (linear scale)", + lin_out_file, false, + ) + + do_plot( + fig, gen_x, func, xlims_log[1], xlims_log[2], + "$fn_name $gen_name (log scale)", + log_out_file, true, + ) +end + +"Create a single plot" +function do_plot( + fig::Figure, + gen_x::Vector{F}, + func::Function, + xmin::AbstractFloat, + xmax::AbstractFloat, + title::String, + out_file::String, + logscale::Bool, +)::Nothing where F<:AbstractFloat + println("plotting $title") + + # `gen_x` is the values the generator produces. `actual_x` is for plotting a + # continuous function. + input_min = xmin - 1.0 + input_max = xmax + 1.0 + gen_x = filter((v) -> v >= input_min && v <= input_max, gen_x) + markersize = length(gen_x) < 10_000 ? 6.0 : 4.0 + + steps = 10_000 + if logscale + r = LinRange(symlog10(input_min), symlog10(input_max), steps) + actual_x = sympow10.(r) + xscale = Makie.pseudolog10 + else + actual_x = LinRange(input_min, input_max, steps) + xscale = identity + end + + gen_y = @. func(gen_x) + actual_y = @. func(actual_x) + + ax = Axis(fig[1, 1], xscale=xscale, title=title) + + lines!( + ax, actual_x, actual_y, color=(:lightblue, 0.6), + linewidth=6.0, label="true function", + ) + scatter!( + ax, gen_x, gen_y, color=(:darkblue, 0.9), + markersize=markersize, label="checked inputs", + ) + axislegend(ax, position=:rb, framevisible=false) + + save(out_file, fig) + delete!(ax) +end + +"Apply a function, returning the default if there is a domain error" +function map_or( + input::AbstractFloat, + f::Function, + default::Any +)::Union{AbstractFloat,Any} + try + return f(input) + catch + return default + end +end + +# Operations for logarithms that are symmetric about 0 +C = 10 +symlog10(x::Number) = sign(x) * (log10(1 + abs(x)/(10^C))) +sympow10(x::Number) = (10^C) * (10^x - 1) + +main()