Add a test against musl libm

Check our functions against `musl-math-sys`. This is similar to the
existing musl tests that go through binary serialization, but works on
more platforms.

crates/libm-test/Cargo.toml

@@ -11,10 +11,15 @@ default = []
 # musl libc.
 test-musl-serialized = ["rand"]
 
+# Build our own musl for testing and benchmarks
+build-musl = ["dep:musl-math-sys"]
+
 [dependencies]
 anyhow = "1.0.90"
 libm = { path = "../.." }
 libm-macros = { path = "../libm-macros" }
+musl-math-sys = { path = "../musl-math-sys", optional = true }
+paste = "1.0.15"
 rand = "0.8.5"
 rand_chacha = "0.3.1"
 

crates/libm-test/src/lib.rs

@@ -14,6 +14,34 @@ pub type TestResult<T = (), E = anyhow::Error> = Result<T, E>;
 // List of all files present in libm's source
 include!(concat!(env!("OUT_DIR"), "/all_files.rs"));
 
+/// ULP allowed to differ from musl (note that musl itself may not be accurate).
+const MUSL_DEFAULT_ULP: u32 = 2;
+
+/// Certain functions have different allowed ULP (consider these xfail).
+///
+/// Note that these results were obtained using 400,000,000 rounds of random inputs, which
+/// is not a value used by default.
+pub fn musl_allowed_ulp(name: &str) -> u32 {
+    match name {
+        #[cfg(x86_no_sse)]
+        "asinh" | "asinhf" => 6,
+        "lgamma" | "lgamma_r" | "lgammaf" | "lgammaf_r" => 400,
+        "tanh" | "tanhf" => 4,
+        "tgamma" => 20,
+        "j0" | "j0f" | "j1" | "j1f" => {
+            // Results seem very target-dependent
+            if cfg!(target_arch = "x86_64") { 4000 } else { 800_000 }
+        }
+        "jn" | "jnf" => 1000,
+        "sincosf" => 500,
+        #[cfg(not(target_pointer_width = "64"))]
+        "exp10" => 4,
+        #[cfg(not(target_pointer_width = "64"))]
+        "exp10f" => 4,
+        _ => MUSL_DEFAULT_ULP,
+    }
+}
+
 /// Return the unsuffixed version of a function name; e.g. `abs` and `absf` both return `abs`,
 /// `lgamma_r` and `lgammaf_r` both return `lgamma_r`.
 pub fn canonical_name(name: &str) -> &str {

crates/libm-test/src/special_case.rs

@@ -1,7 +1,9 @@
 //! Configuration for skipping or changing the result for individual test cases (inputs) rather
 //! than ignoring entire tests.
 
-use crate::{CheckCtx, Float, Int, TestResult};
+use core::f32;
+
+use crate::{CheckBasis, CheckCtx, Float, Int, TestResult};
 
 /// Type implementing [`IgnoreCase`].
 pub struct SpecialCase;
@@ -49,47 +51,189 @@ pub trait MaybeOverride<Input> {
 
 impl MaybeOverride<(f32,)> for SpecialCase {
     fn check_float<F: Float>(
-        _input: (f32,),
+        input: (f32,),
         actual: F,
         expected: F,
         _ulp: &mut u32,
         ctx: &CheckCtx,
     ) -> Option<TestResult> {
+        if ctx.basis == CheckBasis::Musl {
+            if ctx.fname == "acoshf" && input.0 < -1.0 {
+                // acoshf is undefined for x <= 1.0, but we return a random result at lower
+                // values.
+                return XFAIL;
+            }
+
+            if ctx.fname == "sincosf" {
+                let factor_frac_pi_2 = input.0.abs() / f32::consts::FRAC_PI_2;
+                if (factor_frac_pi_2 - factor_frac_pi_2.round()).abs() < 1e-2 {
+                    // we have a bad approximation near multiples of pi/2
+                    return XFAIL;
+                }
+            }
+
+            if ctx.fname == "expm1f" && input.0 > 80.0 && actual.is_infinite() {
+                // we return infinity but the number is representable
+                return XFAIL;
+            }
+
+            if ctx.fname == "sinhf" && 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.fname == "lgammaf" || ctx.fname == "lgammaf_r" && input.0 < 0.0 {
+                // loggamma should not be defined for x < 0, yet we both return results
+                return XFAIL;
+            }
+        }
+
         maybe_check_nan_bits(actual, expected, ctx)
     }
 }
 
 impl MaybeOverride<(f64,)> for SpecialCase {
     fn check_float<F: Float>(
-        _input: (f64,),
+        input: (f64,),
         actual: F,
         expected: F,
         _ulp: &mut u32,
         ctx: &CheckCtx,
     ) -> Option<TestResult> {
+        if ctx.basis == CheckBasis::Musl {
+            if cfg!(target_arch = "x86") && ctx.fname == "acosh" && input.0 < 1.0 {
+                // The function is undefined, both implementations return random results
+                return SKIP;
+            }
+
+            if cfg!(x86_no_sse)
+                && ctx.fname == "ceil"
+                && input.0 < 0.0
+                && input.0 > -1.0
+                && expected == F::ZERO
+                && actual == F::ZERO
+            {
+                // musl returns -0.0, we return +0.0
+                return XFAIL;
+            }
+
+            if ctx.fname == "lgamma" || ctx.fname == "lgamma_r" && input.0 < 0.0 {
+                // loggamma should not be defined for x < 0, yet we both return results
+                return XFAIL;
+            }
+        }
+
         maybe_check_nan_bits(actual, expected, ctx)
     }
 }
 
-impl MaybeOverride<(f32, f32)> for SpecialCase {}
-impl MaybeOverride<(f64, f64)> for SpecialCase {}
-impl MaybeOverride<(f32, f32, f32)> for SpecialCase {}
-impl MaybeOverride<(f64, f64, f64)> for SpecialCase {}
-impl MaybeOverride<(i32, f32)> for SpecialCase {}
-impl MaybeOverride<(i32, f64)> for SpecialCase {}
-impl MaybeOverride<(f32, i32)> for SpecialCase {}
-impl MaybeOverride<(f64, i32)> for SpecialCase {}
-
 /// Check NaN bits if the function requires it
 fn maybe_check_nan_bits<F: Float>(actual: F, expected: F, ctx: &CheckCtx) -> Option<TestResult> {
-    if !(ctx.canonical_name == "abs" || ctx.canonical_name == "copysigh") {
+    if !(ctx.canonical_name == "fabs" || ctx.canonical_name == "copysign") {
         return None;
     }
 
+    // LLVM currently uses x87 instructions which quieten signalling NaNs to handle the i686
+    // `extern "C"` `f32`/`f64` return ABI.
+    // LLVM issue <https://github.com/llvm/llvm-project/issues/66803>
+    // Rust issue <https://github.com/rust-lang/rust/issues/115567>
+    if cfg!(target_arch = "x86") && ctx.basis == CheckBasis::Musl {
+        return SKIP;
+    }
+
     // abs and copysign require signaling NaNs to be propagated, so verify bit equality.
     if actual.to_bits() == expected.to_bits() {
         return SKIP;
     } else {
         Some(Err(anyhow::anyhow!("NaNs have different bitpatterns")))
     }
 }
+
+impl MaybeOverride<(f32, f32)> for SpecialCase {
+    fn check_float<F: Float>(
+        input: (f32, f32),
+        _actual: F,
+        expected: F,
+        _ulp: &mut u32,
+        ctx: &CheckCtx,
+    ) -> Option<TestResult> {
+        maybe_skip_min_max_nan(input, expected, ctx)
+    }
+}
+impl MaybeOverride<(f64, f64)> for SpecialCase {
+    fn check_float<F: Float>(
+        input: (f64, f64),
+        _actual: F,
+        expected: F,
+        _ulp: &mut u32,
+        ctx: &CheckCtx,
+    ) -> Option<TestResult> {
+        maybe_skip_min_max_nan(input, expected, ctx)
+    }
+}
+
+/// Musl propagates NaNs if one is provided as the input, but we return the other input.
+// F1 and F2 are always the same type, this is just to please generics
+fn maybe_skip_min_max_nan<F1: Float, F2: Float>(
+    input: (F1, F1),
+    expected: F2,
+    ctx: &CheckCtx,
+) -> Option<TestResult> {
+    if (ctx.canonical_name == "fmax" || ctx.canonical_name == "fmin")
+        && (input.0.is_nan() || input.1.is_nan())
+        && expected.is_nan()
+    {
+        return XFAIL;
+    } else {
+        None
+    }
+}
+
+impl MaybeOverride<(i32, f32)> for SpecialCase {
+    fn check_float<F: Float>(
+        input: (i32, f32),
+        _actual: F,
+        _expected: F,
+        ulp: &mut u32,
+        ctx: &CheckCtx,
+    ) -> Option<TestResult> {
+        bessel_prec_dropoff(input, ulp, ctx)
+    }
+}
+impl MaybeOverride<(i32, f64)> for SpecialCase {
+    fn check_float<F: Float>(
+        input: (i32, f64),
+        _actual: F,
+        _expected: F,
+        ulp: &mut u32,
+        ctx: &CheckCtx,
+    ) -> Option<TestResult> {
+        bessel_prec_dropoff(input, ulp, ctx)
+    }
+}
+
+/// Our bessel functions blow up with large N values
+fn bessel_prec_dropoff<F: Float>(
+    input: (i32, F),
+    ulp: &mut u32,
+    ctx: &CheckCtx,
+) -> Option<TestResult> {
+    if ctx.canonical_name == "jn" {
+        if input.0 > 4000 {
+            return XFAIL;
+        } else if input.0 > 2000 {
+            // *ulp = 20_000;
+            *ulp = 20000;
+        } else if input.0 > 1000 {
+            *ulp = 4000;
+        }
+    }
+
+    None
+}
+
+impl MaybeOverride<(f32, f32, f32)> for SpecialCase {}
+impl MaybeOverride<(f64, f64, f64)> for SpecialCase {}
+impl MaybeOverride<(f32, i32)> for SpecialCase {}
+impl MaybeOverride<(f64, i32)> for SpecialCase {}

crates/libm-test/src/test_traits.rs

@@ -49,7 +49,10 @@ impl CheckCtx {
 
 /// Possible items to test against
 #[derive(Clone, Debug, PartialEq, Eq)]
-pub enum CheckBasis {}
+pub enum CheckBasis {
+    /// Check against Musl's math sources.
+    Musl,
+}
 
 /// A trait to implement on any output type so we can verify it in a generic way.
 pub trait CheckOutput<Input>: Sized {
@@ -160,8 +163,7 @@ where
 
             // Check when both are NaNs
             if self.is_nan() && expected.is_nan() {
-                ensure!(self.to_bits() == expected.to_bits(), "NaNs have different bitpatterns");
-                // Nothing else to check
+                // By default, NaNs have nothing special to check.
                 return Ok(());
             } else if self.is_nan() || expected.is_nan() {
                 // Check when only one is a NaN

crates/libm-test/tests/compare_built_musl.rs

@@ -0,0 +1,52 @@
+//! Compare our implementations with the result of musl functions, as provided by `musl-math-sys`.
+//!
+//! Currently this only tests randomized inputs. In the future this may be improved to test edge
+//! cases or run exhaustive tests.
+//!
+//! Note that musl functions do not always provide 0.5ULP rounding, so our functions can do better
+//! than these results.
+
+// There are some targets we can't build musl for
+#![cfg(feature = "build-musl")]
+
+use libm_test::gen::random;
+use libm_test::{CheckBasis, CheckCtx, CheckOutput, TupleCall, musl_allowed_ulp};
+use musl_math_sys as musl;
+
+macro_rules! musl_rand_tests {
+    (
+        fn_name: $fn_name:ident,
+        CFn: $CFn:ty,
+        CArgs: $CArgs:ty,
+        CRet: $CRet:ty,
+        RustFn: $RustFn:ty,
+        RustArgs: $RustArgs:ty,
+        RustRet: $RustRet:ty,
+        attrs: [$($meta:meta)*]
+    ) => { paste::paste! {
+        #[test]
+        $(#[$meta])*
+        fn [< musl_random_ $fn_name >]() {
+            let fname = stringify!($fn_name);
+            let ulp = musl_allowed_ulp(fname);
+            let cases = random::get_test_cases::<$RustArgs>(fname);
+            let ctx = CheckCtx::new(ulp, fname, CheckBasis::Musl);
+
+            for input in cases {
+                let musl_res = input.call(musl::$fn_name as $CFn);
+                let crate_res = input.call(libm::$fn_name as $RustFn);
+
+                crate_res.validate(musl_res, input, &ctx).unwrap();
+            }
+        }
+    } };
+}
+
+libm_macros::for_each_function! {
+    callback: musl_rand_tests,
+    skip: [],
+    attributes: [
+        #[cfg_attr(x86_no_sse, ignore)] // FIXME(correctness): wrong result on i586
+        [exp10, exp10f, exp2, exp2f, rint]
+    ],
+}