Add a test against MPFR using random inputs

This commit is contained in:
Trevor Gross
2024-10-21 17:41:24 -05:00
parent c09e58be46
commit 64131ec5cd
7 changed files with 228 additions and 70 deletions

View File

@@ -7,7 +7,7 @@ use rand::{Rng, SeedableRng};
use rand_chacha::ChaCha8Rng;
use super::CachedInput;
use crate::GenerateInput;
use crate::{CheckCtx, GenerateInput};
const SEED: [u8; 32] = *b"3.141592653589793238462643383279";
@@ -40,9 +40,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") {
@@ -105,11 +106,10 @@ 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)
let inputs = if ctx.fname == "jn" || ctx.fname == "jnf" { &TEST_CASES_JN } else { &TEST_CASES };
inputs.get_cases()
}

View File

@@ -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,27 @@ 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,
"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 {

View File

@@ -248,27 +248,6 @@ macro_rules! impl_op_for_ty {
}
}
pub mod [<nextafter $suffix>] {
use super::*;
pub struct Operation(MpFloat, MpFloat);
impl MpOp for Operation {
type Input = ($fty, $fty);
type Output = $fty;
fn new() -> Self {
Self(new_mpfloat::<$fty>(), new_mpfloat::<$fty>())
}
fn run(&mut self, input: Self::Input) -> Self::Output {
self.0.assign(input.0);
self.1.assign(input.1);
self.0.next_toward(&self.1);
prep_retval::<Self::Output>(&mut self.0, Ordering::Equal)
}
}
}
pub mod [<pow $suffix>] {
use super::*;
pub struct Operation(MpFloat, MpFloat);

View File

@@ -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,40 @@ 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
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 +128,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
@@ -142,6 +178,11 @@ fn maybe_check_nan_bits<F: Float>(actual: F, expected: F, ctx: &CheckCtx) -> Opt
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;
@@ -158,9 +199,10 @@ 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 {
fn check_float<F: Float>(
input: (f64, f64),
@@ -169,47 +211,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 {
bessel_prec_dropoff(input, ulp, ctx)
}
}
}
}
}

View File

@@ -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.