Add a test against MPFR using random inputs

tgross35 · tgross35 · commit f36264ac8f71 · 2024-10-28T01:03:50.000-05:00
diff --git a/crates/libm-test/Cargo.toml b/crates/libm-test/Cargo.toml
@@ -18,6 +18,7 @@ build-musl = ["dep:musl-math-sys"]
 [dependencies]
 anyhow = "1.0.90"
 az = { version = "1.2.1", optional = true }
+either = "1.13.0"
 libm = { path = "../.." }
 libm-macros = { path = "../libm-macros" }
 musl-math-sys = { path = "../musl-math-sys", optional = true }
diff --git a/crates/libm-test/src/gen/random.rs b/crates/libm-test/src/gen/random.rs
@@ -3,11 +3,12 @@
 
 use std::sync::LazyLock;
 
+use either::Either;
 use rand::{Rng, SeedableRng};
 use rand_chacha::ChaCha8Rng;
 
 use super::CachedInput;
-use crate::GenerateInput;
+use crate::{CheckBasis, CheckCtx, GenerateInput};
 
 const SEED: [u8; 32] = *b"3.141592653589793238462643383279";
 
@@ -50,9 +51,10 @@ static TEST_CASES_JN: LazyLock<CachedInput> = LazyLock::new(|| {
     let mut cases = (&*TEST_CASES).clone();
 
     // These functions are extremely slow, limit them
-    cases.inputs_i32.truncate((NTESTS / 1000).max(80));
-    cases.inputs_f32.truncate((NTESTS / 1000).max(80));
-    cases.inputs_f64.truncate((NTESTS / 1000).max(80));
+    let ntests_jn = (NTESTS / 1000).max(80);
+    cases.inputs_i32.truncate(ntests_jn);
+    cases.inputs_f32.truncate(ntests_jn);
+    cases.inputs_f64.truncate(ntests_jn);
 
     // It is easy to overflow the stack with these in debug mode
     let max_iterations = if cfg!(optimizations_enabled) && cfg!(target_pointer_width = "64") {
@@ -115,11 +117,17 @@ fn make_test_cases(ntests: usize) -> CachedInput {
 }
 
 /// Create a test case iterator.
-pub fn get_test_cases<RustArgs>(fname: &str) -> impl Iterator<Item = RustArgs>
+pub fn get_test_cases<RustArgs>(ctx: &CheckCtx) -> impl Iterator<Item = RustArgs>
 where
     CachedInput: GenerateInput<RustArgs>,
 {
-    let inputs = if fname == "jn" || fname == "jnf" { &TEST_CASES_JN } else { &TEST_CASES };
-
-    CachedInput::get_cases(inputs)
+    // jn is extremely slow on MPFR, we need to reduce it further
+    let ntests_jn_multi = (NTESTS / 5000).max(20);
+
+    match (&ctx.basis, ctx.fname) {
+        (CheckBasis::Musl, "jnf") => Either::Left(TEST_CASES_JN.get_cases()),
+        (_, "jn") => Either::Left(TEST_CASES_JN.get_cases()),
+        (CheckBasis::Mpfr, "jnf") => Either::Right(TEST_CASES_JN.get_cases().take(ntests_jn_multi)),
+        _ => Either::Left(TEST_CASES.get_cases()),
+    }
 }
diff --git a/crates/libm-test/src/lib.rs b/crates/libm-test/src/lib.rs
@@ -16,14 +16,18 @@ 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).
+/// Default 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).
+/// Default ULP allowed to differ from multiprecision (i.e. infinite) results.
+const MULTIPREC_DEFAULT_ULP: u32 = 1;
+
+/// ULP allowed to differ from muls results.
 ///
 /// 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 {
+    // Consider overrides xfail
     match name {
         #[cfg(x86_no_sse)]
         "asinh" | "asinhf" => 6,
@@ -44,6 +48,28 @@ pub fn musl_allowed_ulp(name: &str) -> u32 {
     }
 }
 
+/// ULP allowed to differ from multiprecision results.
+pub fn multiprec_allowed_ulp(name: &str) -> u32 {
+    // Consider overrides xfail
+    match name {
+        "asinh" | "asinhf" => 2,
+        "acoshf" => 4,
+        "atanh" | "atanhf" => 2,
+        "exp10" | "exp10f" => 3,
+        "j0" | "j0f" | "j1" | "j1f" => {
+            // Results seem very target-dependent
+            if cfg!(target_arch = "x86_64") { 4000 } else { 800_000 }
+        }
+        "jn" | "jnf" => 1000,
+        "lgamma" | "lgammaf" | "lgamma_r" | "lgammaf_r" => 16,
+        "sinh" | "sinhf" => 2,
+        "sincosf" => 500,
+        "tanh" | "tanhf" => 2,
+        "tgamma" => 20,
+        _ => MULTIPREC_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 {
diff --git a/crates/libm-test/src/special_case.rs b/crates/libm-test/src/special_case.rs
@@ -58,20 +58,6 @@ impl MaybeOverride<(f32,)> for SpecialCase {
         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;
@@ -82,15 +68,48 @@ impl MaybeOverride<(f32,)> for SpecialCase {
                 // 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
+        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 == "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 == "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)
     }
+
+    fn check_int<I: Int>(
+        input: (f32,),
+        actual: I,
+        expected: I,
+        ctx: &CheckCtx,
+    ) -> Option<anyhow::Result<()>> {
+        // On MPFR for lgammaf_r, we set -1 as the integer result for negative infinity but MPFR
+        // sets +1
+        if ctx.basis == CheckBasis::Mpfr
+            && ctx.fname == "lgammaf_r"
+            && input.0 == f32::NEG_INFINITY
+            && actual.abs() == expected.abs()
+        {
+            XFAIL
+        } else {
+            None
+        }
+    }
 }
 
 impl MaybeOverride<(f64,)> for SpecialCase {
@@ -117,15 +136,40 @@ impl MaybeOverride<(f64,)> for SpecialCase {
                 // 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;
-            }
+        if ctx.fname == "acosh" && input.0 < 1.0 {
+            // The function is undefined for the inputs, musl and our libm both return
+            // random results.
+            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)
     }
+
+    fn check_int<I: Int>(
+        input: (f64,),
+        actual: I,
+        expected: I,
+        ctx: &CheckCtx,
+    ) -> Option<anyhow::Result<()>> {
+        // On MPFR for lgamma_r, we set -1 as the integer result for negative infinity but MPFR
+        // sets +1
+        if ctx.basis == CheckBasis::Mpfr
+            && ctx.fname == "lgamma_r"
+            && input.0 == f64::NEG_INFINITY
+            && actual.abs() == expected.abs()
+        {
+            XFAIL
+        } else {
+            None
+        }
+    }
 }
 
 /// Check NaN bits if the function requires it
@@ -138,6 +182,12 @@ fn maybe_check_nan_bits<F: Float>(actual: F, expected: F, ctx: &CheckCtx) -> Opt
     if cfg!(target_arch = "x86") && ctx.basis == CheckBasis::Musl && ctx.canonical_name == "fabs" {
         return SKIP;
     }
+
+    // MPFR only has one NaN bitpattern; allow the default `.is_nan()` checks to validate.
+    if ctx.basis == CheckBasis::Mpfr {
+        return SKIP;
+    }
+
     // abs and copysign require signaling NaNs to be propagated, so verify bit equality.
     if actual.to_bits() == expected.to_bits() {
         return SKIP;
@@ -154,7 +204,7 @@ impl MaybeOverride<(f32, f32)> for SpecialCase {
         _ulp: &mut u32,
         ctx: &CheckCtx,
     ) -> Option<TestResult> {
-        maybe_skip_min_max_nan(input, expected, ctx)
+        maybe_skip_binop_nan(input, expected, ctx)
     }
 }
 impl MaybeOverride<(f64, f64)> for SpecialCase {
@@ -165,47 +215,86 @@ impl MaybeOverride<(f64, f64)> for SpecialCase {
         _ulp: &mut u32,
         ctx: &CheckCtx,
     ) -> Option<TestResult> {
-        maybe_skip_min_max_nan(input, expected, ctx)
+        maybe_skip_binop_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>(
+fn maybe_skip_binop_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
+    match ctx.basis {
+        CheckBasis::Musl => {
+            if (ctx.canonical_name == "fmax" || ctx.canonical_name == "fmin")
+                && (input.0.is_nan() || input.1.is_nan())
+                && expected.is_nan()
+            {
+                XFAIL
+            } else {
+                None
+            }
+        }
+        CheckBasis::Mpfr => {
+            if ctx.canonical_name == "copysign" && input.1.is_nan() {
+                SKIP
+            } else {
+                None
+            }
+        }
     }
 }
 
 impl MaybeOverride<(i32, f32)> for SpecialCase {
     fn check_float<F: Float>(
         input: (i32, f32),
-        _actual: F,
-        _expected: F,
+        actual: F,
+        expected: F,
         ulp: &mut u32,
         ctx: &CheckCtx,
     ) -> Option<TestResult> {
-        bessel_prec_dropoff(input, ulp, ctx)
+        match ctx.basis {
+            CheckBasis::Musl => bessel_prec_dropoff(input, ulp, ctx),
+            CheckBasis::Mpfr => {
+                // We return +0.0, MPFR returns -0.0
+                if ctx.fname == "jnf"
+                    && input.1 == f32::NEG_INFINITY
+                    && actual == F::ZERO
+                    && expected == F::ZERO
+                {
+                    XFAIL
+                } else {
+                    None
+                }
+            }
+        }
     }
 }
 impl MaybeOverride<(i32, f64)> for SpecialCase {
     fn check_float<F: Float>(
         input: (i32, f64),
-        _actual: F,
-        _expected: F,
+        actual: F,
+        expected: F,
         ulp: &mut u32,
         ctx: &CheckCtx,
     ) -> Option<TestResult> {
-        bessel_prec_dropoff(input, ulp, ctx)
+        match ctx.basis {
+            CheckBasis::Musl => bessel_prec_dropoff(input, ulp, ctx),
+            CheckBasis::Mpfr => {
+                // We return +0.0, MPFR returns -0.0
+                if ctx.fname == "jn"
+                    && input.1 == f64::NEG_INFINITY
+                    && actual == F::ZERO
+                    && expected == F::ZERO
+                {
+                    XFAIL
+                } else {
+                    None
+                }
+            }
+        }
     }
 }
 
diff --git a/crates/libm-test/src/test_traits.rs b/crates/libm-test/src/test_traits.rs
@@ -52,6 +52,8 @@ impl CheckCtx {
 pub enum CheckBasis {
     /// Check against Musl's math sources.
     Musl,
+    /// Check against infinite precision (MPFR).
+    Mpfr,
 }
 
 /// A trait to implement on any output type so we can verify it in a generic way.
diff --git a/crates/libm-test/tests/compare_built_musl.rs b/crates/libm-test/tests/compare_built_musl.rs
@@ -29,8 +29,8 @@ macro_rules! musl_rand_tests {
         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);
+            let cases = random::get_test_cases::<$RustArgs>(&ctx);
 
             for input in cases {
                 let musl_res = input.call(musl::$fn_name as $CFn);
diff --git a/crates/libm-test/tests/multiprecision.rs b/crates/libm-test/tests/multiprecision.rs

Original file line number	Diff line number	Diff line change
`@@ -52,6 +52,8 @@ impl CheckCtx {`
`52`	`52`	`pub enum CheckBasis {`
`53`	`53`	`/// Check against Musl's math sources.`
`54`	`54`	`Musl,`
	`55`	`+ /// Check against infinite precision (MPFR).`
	`56`	`+ Mpfr,`
`55`	`57`	`}`
`56`	`58`
`57`	`59`	`/// A trait to implement on any output type so we can verify it in a generic way.`