From ea658cb6a3a0619095134cca10c5bedb03ed04a3 Mon Sep 17 00:00:00 2001 From: Joshua Simmons Date: Sun, 16 Oct 2022 17:34:01 +0200 Subject: [PATCH] Fix errors in `sin/cos/tan` exhaustive checking Add fractional ulps error calculation. Fix where the multiplication by pi happens in checking code. Replace rug with small hand-rolled mpfr wrapper. --- Cargo.lock | 19 +- narcissus-maths/Cargo.toml | 2 +- narcissus-maths/src/sin_cos_pi.rs | 3 + narcissus-maths/src/tan_pi.rs | 2 + narcissus-maths/tests/exhaustive_f32.rs | 404 +++++++++++++++++------- narcissus-maths/tests/next_after_f32.rs | 3 +- 6 files changed, 300 insertions(+), 133 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 2351191..810d4b5 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -2,12 +2,6 @@ # It is not intended for manual editing. version = 3 -[[package]] -name = "az" -version = "1.2.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7b7e4c2464d97fe331d41de9d5db0def0a96f4d823b8b32a2efd503578988973" - [[package]] name = "fast-float" version = "0.2.0" @@ -69,7 +63,7 @@ dependencies = [ name = "narcissus-maths" version = "0.1.0" dependencies = [ - "rug", + "gmp-mpfr-sys", ] [[package]] @@ -83,17 +77,6 @@ dependencies = [ name = "renderdoc-sys" version = "0.1.0" -[[package]] -name = "rug" -version = "1.17.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "203180f444c95eac53586ed04793ecf6454c5d28f9eca8eead815fc19e136c47" -dependencies = [ - "az", - "gmp-mpfr-sys", - "libc", -] - [[package]] name = "sdl2-sys" version = "0.1.0" diff --git a/narcissus-maths/Cargo.toml b/narcissus-maths/Cargo.toml index 5065f8e..58e9fa7 100644 --- a/narcissus-maths/Cargo.toml +++ b/narcissus-maths/Cargo.toml @@ -10,4 +10,4 @@ edition = "2021" [dependencies] [dev-dependencies] -rug = "1.17.0" \ No newline at end of file +gmp-mpfr-sys = "1.4.10" \ No newline at end of file diff --git a/narcissus-maths/src/sin_cos_pi.rs b/narcissus-maths/src/sin_cos_pi.rs index ce6de9f..1d213c1 100644 --- a/narcissus-maths/src/sin_cos_pi.rs +++ b/narcissus-maths/src/sin_cos_pi.rs @@ -28,6 +28,9 @@ fn mulsign_f32(x: f32, s: u32) -> f32 { /// Simultaneously computes the sine and cosine of `a` expressed in multiples of *pi* radians, or half-turns. /// +/// Sin error <= 0.96563 ulp. +/// Cos error <= 0.96677 ulp. +/// /// Returns `(sin(a * pi), cos(a * pi))` /// /// # Examples diff --git a/narcissus-maths/src/tan_pi.rs b/narcissus-maths/src/tan_pi.rs index e220f99..3601d38 100644 --- a/narcissus-maths/src/tan_pi.rs +++ b/narcissus-maths/src/tan_pi.rs @@ -20,6 +20,8 @@ const F32_TAN_PI_15_K: [f32; 7] = unsafe { /// /// Returns `tan(a * pi)` /// +/// Error <= 1.60536 ulp. +/// /// # Examples /// /// ``` diff --git a/narcissus-maths/tests/exhaustive_f32.rs b/narcissus-maths/tests/exhaustive_f32.rs index 02295f6..ce723a9 100644 --- a/narcissus-maths/tests/exhaustive_f32.rs +++ b/narcissus-maths/tests/exhaustive_f32.rs @@ -1,65 +1,141 @@ /// Rough translation of Paul Zimmermann's binary32_exhaustive to Rust. /// https://gitlab.inria.fr/zimmerma/math_accuracy/-/blob/master/binary32_exhaustive/check_exhaustive.c use std::{ - ops::MulAssign, + mem::MaybeUninit, sync::atomic::{AtomicUsize, Ordering}, }; use narcissus_maths::{next_after_f32, sin_cos_pi_f32, tan_pi_f32}; -use rug::{ - ops::{AssignRound, DivAssignRound}, - Assign, -}; -fn ulp_error_f32(our_val: f32, ref_val: f32, x: f32, reference: fn(&mut rug::Float)) -> u32 { - if our_val == ref_val { - return 0; - } +use gmp_mpfr_sys::mpfr; + +// Less than 280 cause precision errors in the 2 * pi calculation, introducing extra failures. +const PREC: u32 = 280; +const COUNT: u32 = 2_139_095_040_u32; - let mut our_val = our_val; - let mut ref_val = ref_val; +pub enum Round { + /// Round to nearest + TiesToEven, // MPFR_RNDN, + /// Round toward minus infinity + TowardNegative, // MPFR_RNDD + /// Round toward plus infinity + TowardPositive, // MPFR_RNDU + /// Round toward zero + TowardZero, // MPFR_RNDZ + /// Away from zero. + AwayZero, // MPFR_RNDA +} - if ref_val.is_infinite() { - if our_val.is_infinite() { - // Then ours and reference have different signs. - assert!(our_val * ref_val < 0.0); - return u32::MAX; +impl Round { + fn to_raw(self) -> mpfr::rnd_t { + match self { + Round::TiesToEven => mpfr::rnd_t::RNDN, + Round::TowardNegative => mpfr::rnd_t::RNDD, + Round::TowardPositive => mpfr::rnd_t::RNDU, + Round::TowardZero => mpfr::rnd_t::RNDZ, + Round::AwayZero => mpfr::rnd_t::RNDA, } + } +} - our_val /= 2.0; - let mut ref_val_2 = rug::Float::new(24); - ref_val_2.assign_round(x, rug::float::Round::Down); - reference(&mut ref_val_2); - ref_val_2.div_assign_round(2, rug::float::Round::Down); - ref_val = ref_val_2.to_f32_round(rug::float::Round::Down); - if ref_val.is_infinite() { - ref_val = if ref_val > 0.0 { - f32::from_bits(0x7f00_0000) - } else { - f32::from_bits(0xff00_0000) - }; +pub struct Float(mpfr::mpfr_t); + +impl Float { + pub fn new(precision: u32) -> Self { + unsafe { + let mut x = MaybeUninit::uninit(); + mpfr::init2(x.as_mut_ptr(), precision as i64); + Self(x.assume_init()) } } - if our_val.is_infinite() { - assert!(!ref_val.is_infinite()); - // If ours gives +/-Inf but the correct rounding is in the binary32 range, assume we gave - // +/-2^128. - ref_val /= 2.0; - our_val = if our_val > 0.0 { - f32::from_bits(0x7f00_0000) - } else { - f32::from_bits(0xff00_0000) + pub fn with_value_f32(precision: u32, value: f32, round: Round) -> Self { + unsafe { + let mut x = MaybeUninit::uninit(); + mpfr::init2(x.as_mut_ptr(), precision as i64); + mpfr::set_flt(x.as_mut_ptr(), value, round.to_raw()); + Self(x.assume_init()) } } - let err = our_val - ref_val; - let ulp = next_after_f32(ref_val, our_val) - ref_val; - let err = (err / ulp).abs(); - if err >= u32::MAX as f32 { - u32::MAX - } else { - err as u32 + pub fn to_f32(&self, round: Round) -> f32 { + unsafe { mpfr::get_flt(&self.0, round.to_raw()) } + } + + pub fn to_f64(&self, round: Round) -> f64 { + unsafe { mpfr::get_d(&self.0, round.to_raw()) } + } + + pub fn exp(&self) -> i64 { + unsafe { mpfr::get_exp(&self.0) } + } + + pub fn set_const_pi(&mut self, round: Round) { + let _ = unsafe { mpfr::const_pi(&mut self.0, round.to_raw()) }; + } + + pub fn set_2_exp_u64(&mut self, x: u64, e: i64, round: Round) { + let _ = unsafe { mpfr::set_ui_2exp(&mut self.0, x, e, round.to_raw()) }; + } + + pub fn set_f32(&mut self, x: f32, round: Round) { + let _ = unsafe { mpfr::set_flt(&mut self.0, x, round.to_raw()) }; + } + + pub fn set_precision(&mut self, precision: u32, round: Round) { + unsafe { mpfr::prec_round(&mut self.0, precision as i64, round.to_raw()) }; + } + + pub fn add_assign(&mut self, rhs: &Float, round: Round) { + let _ = unsafe { mpfr::add(&mut self.0, &self.0, &rhs.0, round.to_raw()) }; + } + + pub fn sub_assign(&mut self, rhs: &Float, round: Round) { + let _ = unsafe { mpfr::sub(&mut self.0, &self.0, &rhs.0, round.to_raw()) }; + } + + pub fn mul_assign(&mut self, rhs: &Float, round: Round) { + let _ = unsafe { mpfr::mul(&mut self.0, &self.0, &rhs.0, round.to_raw()) }; + } + + pub fn div_assign(&mut self, rhs: &Float, round: Round) { + let _ = unsafe { mpfr::div(&mut self.0, &self.0, &rhs.0, round.to_raw()) }; + } + + /// multiplies self by 2 raised to `rhs`. + pub fn mul_2_assign_i64(&mut self, rhs: i64, round: Round) { + let _ = unsafe { mpfr::mul_2si(&mut self.0, &self.0, rhs, round.to_raw()) }; + } + + /// divides self by 2 raised to `rhs`. + pub fn div_2_assign_u64(&mut self, rhs: u64, round: Round) { + let _ = unsafe { mpfr::div_2ui(&mut self.0, &self.0, rhs, round.to_raw()) }; + } + + pub fn sin_mut(&mut self, round: Round) -> i32 { + unsafe { mpfr::sin(&mut self.0, &self.0, round.to_raw()) } + } + + pub fn cos_mut(&mut self, round: Round) -> i32 { + unsafe { mpfr::cos(&mut self.0, &self.0, round.to_raw()) } + } + + pub fn tan_mut(&mut self, round: Round) -> i32 { + unsafe { mpfr::tan(&mut self.0, &self.0, round.to_raw()) } + } + + pub fn abs_mut(&mut self, round: Round) -> i32 { + unsafe { mpfr::abs(&mut self.0, &self.0, round.to_raw()) } + } + + pub fn subnormalize(&mut self, t: i32, round: Round) { + let _ = unsafe { mpfr::subnormalize(&mut self.0, t, round.to_raw()) }; + } +} + +impl Drop for Float { + fn drop(&mut self) { + unsafe { mpfr::clear(&mut self.0) } } } @@ -69,66 +145,188 @@ struct TestRange { end: u32, } -#[derive(Clone, Copy, Debug, Default)] +#[derive(Clone, Copy, Default)] struct FloatErrors { num_errors: u32, + num_errors_2: u32, + max_error_f64: f64, max_error_ulp: u32, max_error_val: u32, } -const PREC: u32 = 280; -const COUNT: u32 = 2_139_095_040_u32; +impl std::fmt::Debug for FloatErrors { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + let errors = format!( + "{} ({:.2}%)", + self.num_errors, + (self.num_errors as f64 / (COUNT * 2) as f64) * 100.0 + ); + let errors_2 = format!( + "{} ({:.2}%)", + self.num_errors_2, + (self.num_errors_2 as f64 / (COUNT * 2) as f64) * 100.0 + ); + + f.debug_struct("FloatErrors") + .field("num_errors", &errors) + .field("num_errors_2", &errors_2) + .field("max_error_f64", &self.max_error_f64) + .field("max_error_ulp", &self.max_error_ulp) + .field("max_error_val", &self.max_error_val) + .finish() + } +} struct FloatCheck { - reference: fn(&mut rug::Float), - ours: fn(f32) -> f32, + ref_fn: fn(&mut Float, &Float) -> i32, + our_fn: fn(f32) -> f32, tan_mode: bool, - pi: rug::Float, - tmp: rug::Float, + pi: Float, + tmp: Float, errors: FloatErrors, } impl FloatCheck { - fn new(reference: fn(&mut rug::Float), ours: fn(f32) -> f32, tan_mode: bool) -> Self { + fn new(ref_fn: fn(&mut Float, &Float) -> i32, our_fn: fn(f32) -> f32, tan_mode: bool) -> Self { + unsafe { + mpfr::set_emin(-148); + mpfr::set_emax(128); + } + + let mut pi = Float::new(PREC); + pi.set_const_pi(Round::TiesToEven); + Self { - reference, - ours, + ref_fn, + our_fn, tan_mode, - pi: rug::Float::with_val(PREC, rug::float::Constant::Pi), - tmp: rug::Float::new(PREC), + pi, + tmp: Float::new(PREC), errors: Default::default(), } } + fn ulp_error_f64(&mut self, our_value: f32, x: f32) -> f64 { + let (emin, emax) = unsafe { (mpfr::get_emin(), mpfr::get_emax()) }; + unsafe { + mpfr::set_emin(mpfr::get_emin_min()); + mpfr::set_emax(mpfr::get_emax_max()); + } + + let mut y = Float::new(24); + if our_value.is_infinite() { + y.set_2_exp_u64(1, 128, Round::TowardNegative) + } else { + y.set_f32(our_value, Round::TowardNegative) + } + self.tmp.set_f32(x, Round::TowardNegative); + (self.ref_fn)(&mut self.tmp, &self.pi); + let e = self.tmp.exp(); + self.tmp.sub_assign(&y, Round::AwayZero); + self.tmp.abs_mut(Round::TowardNegative); + y.set_2_exp_u64(1, e - PREC as i64 - 1, Round::TowardNegative); + self.tmp.add_assign(&y, Round::AwayZero); + let e = if e - 24 < -149 { -149 } else { e - 24 }; + self.tmp.mul_2_assign_i64(-e, Round::TowardNegative); + let err = self.tmp.to_f64(Round::AwayZero); + + unsafe { + mpfr::set_emin(emin); + mpfr::set_emax(emax); + } + + err + } + + fn ulp_error(&mut self, our_value: f32, ref_value: f32, x: f32) -> u32 { + if our_value == ref_value { + return 0; + } + + let mut our_value = our_value; + let mut ref_value = ref_value; + + if ref_value.is_infinite() { + if our_value.is_infinite() { + // Then ours and reference have different signs. + assert!(our_value * ref_value < 0.0); + return u32::MAX; + } + + our_value /= 2.0; + self.tmp.set_f32(x, Round::TowardNegative); + (self.ref_fn)(&mut self.tmp, &self.pi); + self.tmp.div_2_assign_u64(1, Round::TowardNegative); + ref_value = self.tmp.to_f32(Round::TowardNegative); + + if ref_value.is_infinite() { + ref_value = if ref_value > 0.0 { + f32::from_bits(0x7f00_0000) + } else { + f32::from_bits(0xff00_0000) + }; + } + } + + if our_value.is_infinite() { + assert!(!ref_value.is_infinite()); + // If ours gives +/-Inf but the correct rounding is in the binary32 range, assume we gave + // +/-2^128. + ref_value /= 2.0; + our_value = if our_value > 0.0 { + f32::from_bits(0x7f00_0000) + } else { + f32::from_bits(0xff00_0000) + } + } + + let err = our_value - ref_value; + let ulp = next_after_f32(ref_value, our_value) - ref_value; + let err = (err / ulp).abs(); + if err >= u32::MAX as f32 { + u32::MAX + } else { + err as u32 + } + } + #[inline(always)] fn check(&mut self, u: u32) { let x = f32::from_bits(u); assert!(x.is_finite()); - let our_val = (self.ours)(x); + let our_value = (self.our_fn)(x); if self.tan_mode { let fract = x.fract(); if fract == 0.5 || fract == -0.5 { - assert!(our_val.is_infinite()); + assert!(our_value.is_infinite()); return; } } - self.tmp.assign(x); - self.tmp.mul_assign(&self.pi); - (self.reference)(&mut self.tmp); - self.tmp.subnormalize_ieee(); - let ref_val = self.tmp.to_f32_round(rug::float::Round::Nearest); + self.tmp.set_f32(x, Round::TiesToEven); + let inex = (self.ref_fn)(&mut self.tmp, &self.pi); + self.tmp.subnormalize(inex, Round::TiesToEven); + let ref_value = self.tmp.to_f32(Round::TiesToEven); - if our_val != ref_val && !(our_val.is_nan() && ref_val.is_nan()) { + if our_value != ref_value && !(our_value.is_nan() && ref_value.is_nan()) { self.errors.num_errors += 1; - let e = ulp_error_f32(our_val, ref_val, x, self.reference); - if e > self.errors.max_error_ulp { - self.errors.max_error_ulp = e; + + let err = self.ulp_error(our_value, ref_value, x); + if err > 1 { + self.errors.num_errors_2 += 1; + } + let err_f64 = self.ulp_error_f64(our_value, x); + + if err > self.errors.max_error_ulp + || (err == self.errors.max_error_ulp && err_f64 > self.errors.max_error_f64) + { + self.errors.max_error_ulp = err; + self.errors.max_error_f64 = err_f64; self.errors.max_error_val = u; } } @@ -136,8 +334,8 @@ impl FloatCheck { } fn check_exhaustive_f32( - reference: fn(&mut rug::Float), - ours: fn(f32) -> f32, + ref_fn: fn(&mut Float, &Float) -> i32, + our_fn: fn(f32) -> f32, tan_mode: bool, ) -> FloatErrors { const SPLIT: u32 = 512; @@ -155,11 +353,10 @@ fn check_exhaustive_f32( }) .collect::>(); - // Try and start with big numbers to avoid reallocs? + // The larger numbers towards the end of the test range are more costly to evaluate. So improve scheduling by + // running those long jobs first. work.reverse(); - let count = AtomicUsize::new(0); - let mut errors = FloatErrors::default(); // Spawn threads. @@ -169,25 +366,16 @@ fn check_exhaustive_f32( let threads = (0..num_threads) .map(|_| { s.spawn(|| { - let mut float_check = FloatCheck::new(reference, ours, tan_mode); + let mut float_check = FloatCheck::new(ref_fn, our_fn, tan_mode); loop { let index = work_index.fetch_add(1, Ordering::SeqCst); if let Some(range) = work.get(index) { - let start = std::time::Instant::now(); for i in range.base..range.end { float_check.check(i); } for i in range.base..range.end { float_check.check(0x8000_0000 + i); } - - let i = count.fetch_add(1, Ordering::SeqCst); - println!( - "{:.1}% chunk {index} took {:?}", - ((i + 1) as f32 * (1.0 / SPLIT as f32)) * 100.0, - std::time::Instant::now() - start - ); - continue; } break; @@ -200,8 +388,13 @@ fn check_exhaustive_f32( for thread in threads { let thread_errors = thread.join().unwrap(); errors.num_errors += thread_errors.num_errors; - if thread_errors.max_error_ulp > errors.max_error_ulp { + errors.num_errors_2 += thread_errors.num_errors_2; + + if thread_errors.max_error_ulp > errors.max_error_ulp + || thread_errors.max_error_f64 > errors.max_error_f64 + { errors.max_error_ulp = thread_errors.max_error_ulp; + errors.max_error_f64 = thread_errors.max_error_f64; errors.max_error_val = thread_errors.max_error_val; } } @@ -210,45 +403,36 @@ fn check_exhaustive_f32( errors } -fn ref_sin_pi_f32(x: &mut rug::Float) { - x.sin_mut(); +fn ref_sin_pi_f32(x: &mut Float, pi: &Float) -> i32 { + x.mul_assign(pi, Round::TiesToEven); + x.sin_mut(Round::TiesToEven) } -fn ref_cos_pi_f32(x: &mut rug::Float) { - x.cos_mut(); +fn ref_cos_pi_f32(x: &mut Float, pi: &Float) -> i32 { + x.mul_assign(pi, Round::TiesToEven); + x.cos_mut(Round::TiesToEven) } -fn ref_tan_pi_f32(x: &mut rug::Float) { - x.tan_mut(); +fn ref_tan_pi_f32(x: &mut Float, pi: &Float) -> i32 { + x.mul_assign(pi, Round::TiesToEven); + x.tan_mut(Round::TiesToEven) } #[test] #[ignore] pub fn exhaustive_sin_pi() { let errors = check_exhaustive_f32(ref_sin_pi_f32, |a| sin_cos_pi_f32(a).0, false); - println!( - "errors: {} ({:.2}%)", - errors.num_errors, - (errors.num_errors as f64 / (COUNT * 2) as f64) * 100.0 - ); - println!("max error ulps: {}", errors.max_error_ulp); - println!("max error error at: {:#08x}", errors.max_error_val); + println!("SIN: {:?}", errors); assert_eq!(errors.max_error_ulp, 1); - assert_eq!(errors.num_errors, 715_786_430); + assert_eq!(errors.num_errors, 744_496_052); } #[test] #[ignore] pub fn exhaustive_cos_pi() { let errors = check_exhaustive_f32(ref_cos_pi_f32, |a| sin_cos_pi_f32(a).1, false); - println!( - "errors: {} ({:.2}%)", - errors.num_errors, - (errors.num_errors as f64 / (COUNT * 2) as f64) * 100.0 - ); - println!("max error ulps: {}", errors.max_error_ulp); - println!("max error error at: {:#08x}", errors.max_error_val); - assert_eq!(errors.num_errors, 33_455_772); + println!("COS: {:?}", errors); + assert_eq!(errors.num_errors, 62_426_004); assert_eq!(errors.max_error_ulp, 1); } @@ -256,13 +440,7 @@ pub fn exhaustive_cos_pi() { #[ignore] pub fn exhaustive_tan_pi() { let errors = check_exhaustive_f32(ref_tan_pi_f32, tan_pi_f32, true); - println!( - "errors: {} ({:.2}%)", - errors.num_errors, - (errors.num_errors as f64 / (COUNT * 2) as f64) * 100.0 - ); - println!("max error ulps: {}", errors.max_error_ulp); - println!("max error error at: {:#08x}", errors.max_error_val); - assert_eq!(errors.num_errors, 68_630_324); + println!("TAN: {:?}", errors); + assert_eq!(errors.num_errors, 100_555_422); assert_eq!(errors.max_error_ulp, 2); } diff --git a/narcissus-maths/tests/next_after_f32.rs b/narcissus-maths/tests/next_after_f32.rs index e03e099..46b7603 100644 --- a/narcissus-maths/tests/next_after_f32.rs +++ b/narcissus-maths/tests/next_after_f32.rs @@ -1,7 +1,8 @@ use narcissus_maths::next_after_f32; mod libc { - use std::ffi::c_float; + use std::os::raw::c_float; + extern "C" { pub fn nextafterf(x: c_float, y: c_float) -> c_float; } -- 2.49.0