From cec3240673b691b2b1472ed3108162c662f53dad Mon Sep 17 00:00:00 2001 From: Srayan Jana Date: Fri, 27 Mar 2026 11:36:34 -0700 Subject: [PATCH 01/31] add sqrt from miri + tests --- src/lib.rs | 82 +++++++++++++++++++++++++++++++++++++++++++++++++++ tests/ieee.rs | 18 +++++++++++ 2 files changed, 100 insertions(+) diff --git a/src/lib.rs b/src/lib.rs index 87de8a0..17ff340 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -484,6 +484,88 @@ pub trait Float: } } + fn sqrt(self) -> Self { + match self.category() { + // preserve zero sign + Category::Zero => self, + // propagate NaN + Category::NaN => self, + // sqrt of negative number is NaN + _ if self.is_negative() => Self::NAN, + // sqrt(∞) = ∞ + Category::Infinity => Self::INFINITY, + Category::Normal => { + // Floating point precision, excluding the integer bit + let prec = i32::try_from(Self::PRECISION).unwrap() - 1; + + // x = 2^(exp - prec) * mant + // where mant is an integer with prec+1 bits + // mant is a u128, which should be large enough for the largest prec (112 for f128) + let mut exp = self.ilogb(); + let mut mant = self.scalbn(prec - exp).to_u128(128).value; + + if exp % 2 != 0 { + // Make exponent even, so it can be divided by 2 + exp -= 1; + mant <<= 1; + } + + // Bit-by-bit (base-2 digit-by-digit) sqrt of mant. + // mant is treated here as a fixed point number with prec fractional bits. + // mant will be shifted left by one bit to have an extra fractional bit, which + // will be used to determine the rounding direction. + + // res is the truncated sqrt of mant, where one bit is added at each iteration. + let mut res = 0u128; + // rem is the remainder with the current res + // rem_i = 2^i * ((mant<<1) - res_i^2) + // starting with res = 0, rem = mant<<1 + let mut rem = mant << 1; + // s_i = 2*res_i + let mut s = 0u128; + // d is used to iterate over bits, from high to low (d_i = 2^(-i)) + let mut d = 1u128 << (prec + 1); + + // For iteration j=i+1, we need to find largest b_j = 0 or 1 such that + // (res_i + b_j * 2^(-j))^2 <= mant<<1 + // Expanding (a + b)^2 = a^2 + b^2 + 2*a*b: + // res_i^2 + (b_j * 2^(-j))^2 + 2 * res_i * b_j * 2^(-j) <= mant<<1 + // And rearranging the terms: + // b_j^2 * 2^(-j) + 2 * res_i * b_j <= 2^j * (mant<<1 - res_i^2) + // b_j^2 * 2^(-j) + 2 * res_i * b_j <= rem_i + + while d != 0 { + // Probe b_j^2 * 2^(-j) + 2 * res_i * b_j <= rem_i with b_j = 1: + // t = 2*res_i + 2^(-j) + let t = s + d; + if rem >= t { + // b_j should be 1, so make res_j = res_i + 2^(-j) and adjust rem + res += d; + s += d + d; + rem -= t; + } + // Adjust rem for next iteration + rem <<= 1; + // Shift iterator + d >>= 1; + } + + // Remove extra fractional bit from result, rounding to nearest. + // If the last bit is 0, then the nearest neighbor is definitely the lower one. + // If the last bit is 1, it sounds like this may either be a tie (if there's + // infinitely many 0s after this 1), or the nearest neighbor is the upper one. + // However, since square roots are either exact or irrational, and an exact root + // would lead to the last "extra" bit being 0, we can exclude a tie in this case. + // We therefore always round up if the last bit is 1. When the last bit is 0, + // adding 1 will not do anything since the shift will discard it. + res = (res + 1) >> 1; + + // Build resulting value with res as mantissa and exp/2 as exponent + Self::from_u128(res).value.scalbn(exp / 2 - prec) + } + } + } + /// IEEE-754R isSignMinus: Returns true if and only if the current value is /// negative. /// diff --git a/tests/ieee.rs b/tests/ieee.rs index 0c356d9..f646ff4 100644 --- a/tests/ieee.rs +++ b/tests/ieee.rs @@ -795,6 +795,24 @@ fn maximum() { assert!(nan.maximum(f1).to_f64().is_nan()); } +#[test] +fn sqrt() { + assert_eq!(64_f32.sqrt(), 8_f32); + assert_eq!(64_f64.sqrt(), 8_f64); + assert_eq!(f32::INFINITY.sqrt(), f32::INFINITY); + assert_eq!(f64::INFINITY.sqrt(), f64::INFINITY); + assert_eq!(0.0_f32.sqrt().total_cmp(&0.0), std::cmp::Ordering::Equal); + assert_eq!(0.0_f64.sqrt().total_cmp(&0.0), std::cmp::Ordering::Equal); + assert_eq!((-0.0_f32).sqrt().total_cmp(&-0.0), std::cmp::Ordering::Equal); + assert_eq!((-0.0_f64).sqrt().total_cmp(&-0.0), std::cmp::Ordering::Equal); + assert!((-5.0_f32).sqrt().is_nan()); + assert!((-5.0_f64).sqrt().is_nan()); + assert!(f32::NEG_INFINITY.sqrt().is_nan()); + assert!(f64::NEG_INFINITY.sqrt().is_nan()); + assert!(f32::NAN.sqrt().is_nan()); + assert!(f64::NAN.sqrt().is_nan()); +} + #[test] fn denormal() { // Test single precision From 9108da8127cfc472dd9dd81ab324068514b0dd02 Mon Sep 17 00:00:00 2001 From: Srayan Jana Date: Fri, 27 Mar 2026 12:40:09 -0700 Subject: [PATCH 02/31] actually test double --- tests/ieee.rs | 27 +++++++++++++-------------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/tests/ieee.rs b/tests/ieee.rs index f646ff4..3d84f0f 100644 --- a/tests/ieee.rs +++ b/tests/ieee.rs @@ -2,6 +2,7 @@ extern crate rustc_apfloat; use core::cmp::Ordering; +use std::ops::Neg; use rustc_apfloat::ieee::{BFloat, Double, Float8E4M3FN, Float8E5M2, Half, Quad, Single, X87DoubleExtended}; use rustc_apfloat::{Category, ExpInt, IEK_INF, IEK_NAN, IEK_ZERO}; use rustc_apfloat::{Float, FloatConvert, Round, Status}; @@ -797,20 +798,18 @@ fn maximum() { #[test] fn sqrt() { - assert_eq!(64_f32.sqrt(), 8_f32); - assert_eq!(64_f64.sqrt(), 8_f64); - assert_eq!(f32::INFINITY.sqrt(), f32::INFINITY); - assert_eq!(f64::INFINITY.sqrt(), f64::INFINITY); - assert_eq!(0.0_f32.sqrt().total_cmp(&0.0), std::cmp::Ordering::Equal); - assert_eq!(0.0_f64.sqrt().total_cmp(&0.0), std::cmp::Ordering::Equal); - assert_eq!((-0.0_f32).sqrt().total_cmp(&-0.0), std::cmp::Ordering::Equal); - assert_eq!((-0.0_f64).sqrt().total_cmp(&-0.0), std::cmp::Ordering::Equal); - assert!((-5.0_f32).sqrt().is_nan()); - assert!((-5.0_f64).sqrt().is_nan()); - assert!(f32::NEG_INFINITY.sqrt().is_nan()); - assert!(f64::NEG_INFINITY.sqrt().is_nan()); - assert!(f32::NAN.sqrt().is_nan()); - assert!(f64::NAN.sqrt().is_nan()); + let f1 = Double::from_f64(64.); + let f2 = Double::from_f64(8.); + let infinity = Double::INFINITY; + let nan = Double::NAN; + let negative_infinity = Double::INFINITY.neg(); + assert_eq!(f1.sqrt().to_f64(), f2.to_f64()); + assert_eq!(infinity.sqrt().to_f64(), infinity.to_f64()); + assert_eq!(Double::ZERO.sqrt().to_f64().total_cmp(&0.0), std::cmp::Ordering::Equal); + assert_eq!((-Double::ZERO).sqrt().to_f64().total_cmp(&-0.0), std::cmp::Ordering::Equal); + assert!((-Double::from_f64(5.0)).sqrt().is_nan()); + assert!(negative_infinity.sqrt().is_nan()); + assert!(nan.sqrt().is_nan()); } #[test] From 0a204b9f8c7950bb9675884ae9f3300a4d760620 Mon Sep 17 00:00:00 2001 From: Srayan Jana Date: Fri, 27 Mar 2026 12:41:30 -0700 Subject: [PATCH 03/31] formatting --- tests/ieee.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/ieee.rs b/tests/ieee.rs index 3d84f0f..c1e4734 100644 --- a/tests/ieee.rs +++ b/tests/ieee.rs @@ -2,10 +2,10 @@ extern crate rustc_apfloat; use core::cmp::Ordering; -use std::ops::Neg; use rustc_apfloat::ieee::{BFloat, Double, Float8E4M3FN, Float8E5M2, Half, Quad, Single, X87DoubleExtended}; use rustc_apfloat::{Category, ExpInt, IEK_INF, IEK_NAN, IEK_ZERO}; use rustc_apfloat::{Float, FloatConvert, Round, Status}; +use std::ops::Neg; // FIXME(eddyb) maybe include this in `rustc_apfloat` itself? macro_rules! define_for_each_float_type { @@ -807,7 +807,7 @@ fn sqrt() { assert_eq!(infinity.sqrt().to_f64(), infinity.to_f64()); assert_eq!(Double::ZERO.sqrt().to_f64().total_cmp(&0.0), std::cmp::Ordering::Equal); assert_eq!((-Double::ZERO).sqrt().to_f64().total_cmp(&-0.0), std::cmp::Ordering::Equal); - assert!((-Double::from_f64(5.0)).sqrt().is_nan()); + assert!(Double::from_f64(-5.0).sqrt().is_nan()); assert!(negative_infinity.sqrt().is_nan()); assert!(nan.sqrt().is_nan()); } From b16d17fcf2a2c4a2984d57de31a93eaf993d88e4 Mon Sep 17 00:00:00 2001 From: Srayan Jana Date: Fri, 27 Mar 2026 12:43:13 -0700 Subject: [PATCH 04/31] remove unnecessary variables --- tests/ieee.rs | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/tests/ieee.rs b/tests/ieee.rs index c1e4734..62da4ee 100644 --- a/tests/ieee.rs +++ b/tests/ieee.rs @@ -800,16 +800,13 @@ fn maximum() { fn sqrt() { let f1 = Double::from_f64(64.); let f2 = Double::from_f64(8.); - let infinity = Double::INFINITY; - let nan = Double::NAN; - let negative_infinity = Double::INFINITY.neg(); assert_eq!(f1.sqrt().to_f64(), f2.to_f64()); - assert_eq!(infinity.sqrt().to_f64(), infinity.to_f64()); + assert_eq!(Double::INFINITY.sqrt().to_f64(), f64::INFINITY); assert_eq!(Double::ZERO.sqrt().to_f64().total_cmp(&0.0), std::cmp::Ordering::Equal); assert_eq!((-Double::ZERO).sqrt().to_f64().total_cmp(&-0.0), std::cmp::Ordering::Equal); assert!(Double::from_f64(-5.0).sqrt().is_nan()); - assert!(negative_infinity.sqrt().is_nan()); - assert!(nan.sqrt().is_nan()); + assert!(Double::INFINITY.neg().sqrt().is_nan()); + assert!(Double::NAN.sqrt().is_nan()); } #[test] From 816b4ea32b8405d7bdb23bfa068c5da9c0f8fd18 Mon Sep 17 00:00:00 2001 From: Srayan Jana Date: Fri, 27 Mar 2026 12:53:52 -0700 Subject: [PATCH 05/31] add check against host cpu sqrt --- tests/ieee.rs | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/ieee.rs b/tests/ieee.rs index 62da4ee..02e3eae 100644 --- a/tests/ieee.rs +++ b/tests/ieee.rs @@ -801,6 +801,8 @@ fn sqrt() { let f1 = Double::from_f64(64.); let f2 = Double::from_f64(8.); assert_eq!(f1.sqrt().to_f64(), f2.to_f64()); + assert_eq!(f1.sqrt().to_f64(), 64_f64.sqrt()); + assert_eq!(f2.sqrt().to_f64(), 8_f64.sqrt()); assert_eq!(Double::INFINITY.sqrt().to_f64(), f64::INFINITY); assert_eq!(Double::ZERO.sqrt().to_f64().total_cmp(&0.0), std::cmp::Ordering::Equal); assert_eq!((-Double::ZERO).sqrt().to_f64().total_cmp(&-0.0), std::cmp::Ordering::Equal); From b2fe989bab69a86da3bca9e45851be59c26df2d1 Mon Sep 17 00:00:00 2001 From: Srayan Jana Date: Fri, 27 Mar 2026 13:10:59 -0700 Subject: [PATCH 06/31] check all float versions --- tests/ieee.rs | 33 +++++++++++++++++++++------------ 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/tests/ieee.rs b/tests/ieee.rs index 02e3eae..cfb40e3 100644 --- a/tests/ieee.rs +++ b/tests/ieee.rs @@ -5,7 +5,6 @@ use core::cmp::Ordering; use rustc_apfloat::ieee::{BFloat, Double, Float8E4M3FN, Float8E5M2, Half, Quad, Single, X87DoubleExtended}; use rustc_apfloat::{Category, ExpInt, IEK_INF, IEK_NAN, IEK_ZERO}; use rustc_apfloat::{Float, FloatConvert, Round, Status}; -use std::ops::Neg; // FIXME(eddyb) maybe include this in `rustc_apfloat` itself? macro_rules! define_for_each_float_type { @@ -798,17 +797,27 @@ fn maximum() { #[test] fn sqrt() { - let f1 = Double::from_f64(64.); - let f2 = Double::from_f64(8.); - assert_eq!(f1.sqrt().to_f64(), f2.to_f64()); - assert_eq!(f1.sqrt().to_f64(), 64_f64.sqrt()); - assert_eq!(f2.sqrt().to_f64(), 8_f64.sqrt()); - assert_eq!(Double::INFINITY.sqrt().to_f64(), f64::INFINITY); - assert_eq!(Double::ZERO.sqrt().to_f64().total_cmp(&0.0), std::cmp::Ordering::Equal); - assert_eq!((-Double::ZERO).sqrt().to_f64().total_cmp(&-0.0), std::cmp::Ordering::Equal); - assert!(Double::from_f64(-5.0).sqrt().is_nan()); - assert!(Double::INFINITY.neg().sqrt().is_nan()); - assert!(Double::NAN.sqrt().is_nan()); + for_each_float_type!(for test::()); + fn test() { + assert!(F::ZERO.sqrt().bitwise_eq(F::ZERO)); + assert!((-F::ZERO).sqrt().bitwise_eq(-F::ZERO)); + assert!(F::INFINITY.sqrt().bitwise_eq(F::INFINITY)); + assert!(F::NAN.sqrt().is_nan()); + assert!((-F::INFINITY).sqrt().is_nan()); + assert!((-F::from_u128(5).value).sqrt().is_nan()); + let one = F::from_u128(1).value; + assert!(one.sqrt().bitwise_eq(one)); + let f1 = F::from_u128(64).value; + let f2 = F::from_u128(8).value; + assert!(f1.sqrt().bitwise_eq(f2)); + } + + // Cross-check against host sqrt + // TODO: This is only f32 and f64 for now since they are stable, update to other values later. + for x in 0..1000 { + assert_eq!(Single::from_u128(x).value.sqrt().to_f32(), f32::sqrt(x as f32)); + assert_eq!(Double::from_u128(x).value.sqrt().to_f64(), f64::sqrt(x as f64)); + } } #[test] From 3fc058aedcc64b8f10e248d4c36e0dd24ee8dd52 Mon Sep 17 00:00:00 2001 From: Srayan Jana Date: Fri, 27 Mar 2026 13:12:02 -0700 Subject: [PATCH 07/31] minor formatting --- tests/ieee.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/ieee.rs b/tests/ieee.rs index cfb40e3..ab1efe9 100644 --- a/tests/ieee.rs +++ b/tests/ieee.rs @@ -814,7 +814,7 @@ fn sqrt() { // Cross-check against host sqrt // TODO: This is only f32 and f64 for now since they are stable, update to other values later. - for x in 0..1000 { + for x in 0..1000 { assert_eq!(Single::from_u128(x).value.sqrt().to_f32(), f32::sqrt(x as f32)); assert_eq!(Double::from_u128(x).value.sqrt().to_f64(), f64::sqrt(x as f64)); } From f08e8bea811cbbb2a80d8902a6e7167781d252c5 Mon Sep 17 00:00:00 2001 From: Srayan Jana Date: Wed, 1 Apr 2026 16:14:29 -0700 Subject: [PATCH 08/31] added comment documenting sqrt behavior --- src/lib.rs | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/lib.rs b/src/lib.rs index 17ff340..ae64147 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -484,6 +484,8 @@ pub trait Float: } } + /// Technically this should be StatusAnd, but since sqrt is exact for all supported formats, + /// we can just round to the nearest, ignore exceptions, and return Self. fn sqrt(self) -> Self { match self.category() { // preserve zero sign From 059fc87f3b97f16d6c13f2c7a219922315f31bd2 Mon Sep 17 00:00:00 2001 From: Srayan Jana Date: Wed, 1 Apr 2026 16:34:15 -0700 Subject: [PATCH 09/31] push and hope fuzz gets triggered by CI --- fuzz/ops.rs | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/fuzz/ops.rs b/fuzz/ops.rs index fb1a45d..6522324 100644 --- a/fuzz/ops.rs +++ b/fuzz/ops.rs @@ -16,6 +16,9 @@ enum OpKind { Binary(char), Ternary(Rust<&'static str>, Cxx<&'static str>), + // method-call ops, this one has the arity of 1 (so only 1 argument) + RustMethodArity1(Rust<&'static str>), + // HACK(eddyb) all other ops have floating-point inputs *and* outputs, so // the easiest way to fuzz conversions from/to other types, even if it won't // cover *all possible* inputs, is to do a round-trip through the other type. @@ -41,7 +44,7 @@ impl Type { impl OpKind { fn inputs<'a, T>(&self, all_inputs: &'a [T; 3]) -> &'a [T] { match self { - Unary(_) | Roundtrip(_) => &all_inputs[..1], + Unary(_) | RustMethodArity1(_) | Roundtrip(_) => &all_inputs[..1], Binary(_) => &all_inputs[..2], Ternary(..) => &all_inputs[..3], } @@ -59,6 +62,9 @@ const OPS: &[(&str, OpKind)] = &[ ("Rem", Binary('%')), // Ternary (`(F, F) -> F`) ops. ("MulAdd", Ternary(Rust("mul_add"), Cxx("fusedMultiplyAdd"))), + // Method-call ops. + // For now, sqrt is Rust-only, there is no C++ `APFloat` equivalent + ("Sqrt", RustMethodArity1(Rust("sqrt"))), // Roundtrip (`F -> T -> F`) ops. ("FToI128ToF", Roundtrip(Type::SInt(128))), ("FToU128ToF", Roundtrip(Type::UInt(128))), @@ -154,6 +160,9 @@ impl FuzzOp Ternary(Rust(method), _) => { format!("{}.{method}({}, {})", inputs[0], inputs[1], inputs[2]) } + RustMethodArity1(Rust(method)) => { + format!("{}.{method}()", inputs[0]) + } Roundtrip(ty) => format!( "<{ty} as num_traits::AsPrimitive::>::as_( >::as_({}))", @@ -189,6 +198,9 @@ impl FuzzOp Ternary(Rust(method), _) => { format!("{}.{method}({}).value", inputs[0], inputs[1..].join(", ")) } + RustMethodArity1(Rust(method)) => { + format!("{}.{method}()", inputs[0]) + } Roundtrip(ty @ (Type::SInt(_) | Type::UInt(_))) => { let (w, i_or_u) = match ty { Type::SInt(w) => (w, "i"), @@ -266,6 +278,16 @@ struct FuzzOp { + &all_ops_map_concat(|_tag, name, kind| { let inputs = kind.inputs(&["a.to_apf()", "b.to_apf()", "c.to_apf()"]); let expr = match kind { + RustMethodArity1(method_name) => { + if method_name.0 == "sqrt" { + // For now, sqrt is the only Rust method-call op, and it has no C++ `APFloat` equivalent + // so don't generate any C++ code for it. + return String::new(); + } else { + unreachable!() + } + } + // HACK(eddyb) `APFloat` doesn't overload `operator%`, so we have // to go through the `mod` method instead. Binary('%') => format!("((r = {}), r.mod({}), r)", inputs[0], inputs[1]), From 1deb61a9c858013cee8dafa63988ce2315141451 Mon Sep 17 00:00:00 2001 From: Srayan Jana Date: Wed, 1 Apr 2026 16:41:00 -0700 Subject: [PATCH 10/31] formatting --- fuzz/ops.rs | 2 +- src/lib.rs | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/fuzz/ops.rs b/fuzz/ops.rs index 6522324..1d8f0b5 100644 --- a/fuzz/ops.rs +++ b/fuzz/ops.rs @@ -286,7 +286,7 @@ struct FuzzOp { } else { unreachable!() } - } + } // HACK(eddyb) `APFloat` doesn't overload `operator%`, so we have // to go through the `mod` method instead. diff --git a/src/lib.rs b/src/lib.rs index ae64147..d5afa3f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -484,7 +484,7 @@ pub trait Float: } } - /// Technically this should be StatusAnd, but since sqrt is exact for all supported formats, + /// Technically this should be StatusAnd, but since sqrt is exact for all supported formats, /// we can just round to the nearest, ignore exceptions, and return Self. fn sqrt(self) -> Self { match self.category() { From b6404fe4d046c5fe7f2e73057f26b09a596d76eb Mon Sep 17 00:00:00 2001 From: Srayan Jana Date: Wed, 1 Apr 2026 17:28:17 -0700 Subject: [PATCH 11/31] add comment --- src/lib.rs | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/lib.rs b/src/lib.rs index d5afa3f..31000c0 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -491,6 +491,13 @@ pub trait Float: // preserve zero sign Category::Zero => self, // propagate NaN + // If the input is a signalling NaN, then IEEE 754 requires the result to be converted to a quiet NaN. + // On most CPUs that means the most significant bit of the significand field is 0 for signalling NaNs and 1 for quiet NaNs. + // On most CPUs they quiet a NaN by setting that bit to a 1, RISC-V instead returns the canonical NaN with positive sign, + // the most significant significand bit set and all other significand bits cleared. + // However, Rust and LLVM allow input NaNs to be returned unmodified as well as a few other options -- see Rust's rules for NaNs. + // https://doc.rust-lang.org/std/primitive.f32.html#nan-bit-patterns + // (Thanks @programmerjake for the comment) Category::NaN => self, // sqrt of negative number is NaN _ if self.is_negative() => Self::NAN, From 64387f74b2c9173323ae17be5a7824f87502aa8b Mon Sep 17 00:00:00 2001 From: Srayan Jana Date: Wed, 1 Apr 2026 21:01:42 -0700 Subject: [PATCH 12/31] code cleanup, make sqrt function ieee only --- fuzz/ops.rs | 14 ++++---- src/ieee.rs | 89 +++++++++++++++++++++++++++++++++++++++++++++++++ src/lib.rs | 91 ++------------------------------------------------- src/ppc.rs | 4 +++ tests/ieee.rs | 39 +++++++++++++++++++++- 5 files changed, 141 insertions(+), 96 deletions(-) diff --git a/fuzz/ops.rs b/fuzz/ops.rs index 1d8f0b5..3befe42 100644 --- a/fuzz/ops.rs +++ b/fuzz/ops.rs @@ -16,8 +16,8 @@ enum OpKind { Binary(char), Ternary(Rust<&'static str>, Cxx<&'static str>), - // method-call ops, this one has the arity of 1 (so only 1 argument) - RustMethodArity1(Rust<&'static str>), + // A operation that's actually a method call in Rust, but only takes one argument (unary). + RustUnary(Rust<&'static str>), // HACK(eddyb) all other ops have floating-point inputs *and* outputs, so // the easiest way to fuzz conversions from/to other types, even if it won't @@ -44,7 +44,7 @@ impl Type { impl OpKind { fn inputs<'a, T>(&self, all_inputs: &'a [T; 3]) -> &'a [T] { match self { - Unary(_) | RustMethodArity1(_) | Roundtrip(_) => &all_inputs[..1], + Unary(_) | RustUnary(_) | Roundtrip(_) => &all_inputs[..1], Binary(_) => &all_inputs[..2], Ternary(..) => &all_inputs[..3], } @@ -64,7 +64,7 @@ const OPS: &[(&str, OpKind)] = &[ ("MulAdd", Ternary(Rust("mul_add"), Cxx("fusedMultiplyAdd"))), // Method-call ops. // For now, sqrt is Rust-only, there is no C++ `APFloat` equivalent - ("Sqrt", RustMethodArity1(Rust("sqrt"))), + ("Sqrt", RustUnary(Rust("sqrt"))), // Roundtrip (`F -> T -> F`) ops. ("FToI128ToF", Roundtrip(Type::SInt(128))), ("FToU128ToF", Roundtrip(Type::UInt(128))), @@ -160,7 +160,7 @@ impl FuzzOp Ternary(Rust(method), _) => { format!("{}.{method}({}, {})", inputs[0], inputs[1], inputs[2]) } - RustMethodArity1(Rust(method)) => { + RustUnary(Rust(method)) => { format!("{}.{method}()", inputs[0]) } Roundtrip(ty) => format!( @@ -198,7 +198,7 @@ impl FuzzOp Ternary(Rust(method), _) => { format!("{}.{method}({}).value", inputs[0], inputs[1..].join(", ")) } - RustMethodArity1(Rust(method)) => { + RustUnary(Rust(method)) => { format!("{}.{method}()", inputs[0]) } Roundtrip(ty @ (Type::SInt(_) | Type::UInt(_))) => { @@ -278,7 +278,7 @@ struct FuzzOp { + &all_ops_map_concat(|_tag, name, kind| { let inputs = kind.inputs(&["a.to_apf()", "b.to_apf()", "c.to_apf()"]); let expr = match kind { - RustMethodArity1(method_name) => { + RustUnary(method_name) => { if method_name.0 == "sqrt" { // For now, sqrt is the only Rust method-call op, and it has no C++ `APFloat` equivalent // so don't generate any C++ code for it. diff --git a/src/ieee.rs b/src/ieee.rs index 9b8c5a0..eeef2ec 100644 --- a/src/ieee.rs +++ b/src/ieee.rs @@ -1563,6 +1563,95 @@ impl Float for IeeeFloat { } } + fn sqrt(self) -> Self { + match self.category() { + // preserve zero sign + Category::Zero => return self, + // propagate NaN + // If the input is a signalling NaN, then IEEE 754 requires the result to be converted to a quiet NaN. + // On most CPUs that means the most significant bit of the significand field is 0 for signalling NaNs and 1 for quiet NaNs. + // On most CPUs they quiet a NaN by setting that bit to a 1, RISC-V instead returns the canonical NaN with positive sign, + // the most significant significand bit set and all other significand bits cleared. + // However, Rust and LLVM allow input NaNs to be returned unmodified as well as a few other options -- see Rust's rules for NaNs. + // https://doc.rust-lang.org/std/primitive.f32.html#nan-bit-patterns + // (Thanks @programmerjake for the comment) + Category::NaN => return IeeeDefaultExceptionHandling::result_from_nan(self).value, + // sqrt of negative number is NaN + _ if self.is_negative() => return Self::NAN, + // sqrt(inf) = inf + Category::Infinity => return Self::INFINITY, + Category::Normal => (), + } + + // Floating point precision, excluding the integer bit. + let prec = i32::try_from(Self::PRECISION).unwrap() - 1; + + // x = 2^(exp - prec) * mant + // where mant is an integer with prec+1 bits. + // mant is a u128, which is large enough for the largest prec (112 for f128). + let mut exp = self.ilogb(); + let mut mant = self.scalbn(prec - exp).to_u128(128).value; + + if exp % 2 != 0 { + // Make exponent even, so it can be divided by 2. + exp -= 1; + mant <<= 1; + } + + // Bit-by-bit (base-2 digit-by-digit) sqrt of mant. + // mant is treated here as a fixed point number with prec fractional bits. + // mant will be shifted left by one bit to have an extra fractional bit, which + // will be used to determine the rounding direction. + + // res is the truncated sqrt of mant, where one bit is added at each iteration. + let mut res = 0u128; + // rem is the remainder with the current res + // rem_i = 2^i * ((mant<<1) - res_i^2) + // starting with res = 0, rem = mant<<1 + let mut rem = mant << 1; + // s_i = 2*res_i + let mut s = 0u128; + // d is used to iterate over bits, from high to low (d_i = 2^(-i)) + let mut d = 1u128 << (prec + 1); + + // For iteration j=i+1, we need to find largest b_j = 0 or 1 such that + // (res_i + b_j * 2^(-j))^2 <= mant<<1 + // Expanding (a + b)^2 = a^2 + b^2 + 2*a*b: + // res_i^2 + (b_j * 2^(-j))^2 + 2 * res_i * b_j * 2^(-j) <= mant<<1 + // And rearranging the terms: + // b_j^2 * 2^(-j) + 2 * res_i * b_j <= 2^j * (mant<<1 - res_i^2) + // b_j^2 * 2^(-j) + 2 * res_i * b_j <= rem_i + + while d != 0 { + // Probe b_j^2 * 2^(-j) + 2 * res_i * b_j <= rem_i with b_j = 1: + // t = 2*res_i + 2^(-j) + let t = s + d; + if rem >= t { + // b_j should be 1, so make res_j = res_i + 2^(-j) and adjust rem + res += d; + s += d + d; + rem -= t; + } + // Adjust rem for next iteration + rem <<= 1; + // Shift iterator + d >>= 1; + } + + // Remove extra fractional bit from result, rounding to nearest. + // If the last bit is 0, then the nearest neighbor is definitely the lower one. + // If the last bit is 1, it sounds like this may either be a tie (if there's + // infinitely many 0s after this 1), or the nearest neighbor is the upper one. + // However, since square roots are either exact or irrational, and an exact root + // would lead to the last "extra" bit being 0, we can exclude a tie in this case. + // We therefore always round up if the last bit is 1. When the last bit is 0, + // adding 1 will not do anything since the shift will discard it. + res = (res + 1) >> 1; + + // Build resulting value with res as mantissa and exp/2 as exponent. + Self::from_u128(res).value.scalbn(exp / 2 - prec) + } + fn from_bits(input: u128) -> Self { // Dispatch to semantics. S::from_bits(input) diff --git a/src/lib.rs b/src/lib.rs index 31000c0..9a4c8c7 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -484,95 +484,10 @@ pub trait Float: } } - /// Technically this should be StatusAnd, but since sqrt is exact for all supported formats, - /// we can just round to the nearest, ignore exceptions, and return Self. + /// IEEE-754R sqrt: Returns the correctly rounded square root of the current value + /// Note: we currently don't support raising any exceptions from sqrt, so the result is always exact and the status is always OK. fn sqrt(self) -> Self { - match self.category() { - // preserve zero sign - Category::Zero => self, - // propagate NaN - // If the input is a signalling NaN, then IEEE 754 requires the result to be converted to a quiet NaN. - // On most CPUs that means the most significant bit of the significand field is 0 for signalling NaNs and 1 for quiet NaNs. - // On most CPUs they quiet a NaN by setting that bit to a 1, RISC-V instead returns the canonical NaN with positive sign, - // the most significant significand bit set and all other significand bits cleared. - // However, Rust and LLVM allow input NaNs to be returned unmodified as well as a few other options -- see Rust's rules for NaNs. - // https://doc.rust-lang.org/std/primitive.f32.html#nan-bit-patterns - // (Thanks @programmerjake for the comment) - Category::NaN => self, - // sqrt of negative number is NaN - _ if self.is_negative() => Self::NAN, - // sqrt(∞) = ∞ - Category::Infinity => Self::INFINITY, - Category::Normal => { - // Floating point precision, excluding the integer bit - let prec = i32::try_from(Self::PRECISION).unwrap() - 1; - - // x = 2^(exp - prec) * mant - // where mant is an integer with prec+1 bits - // mant is a u128, which should be large enough for the largest prec (112 for f128) - let mut exp = self.ilogb(); - let mut mant = self.scalbn(prec - exp).to_u128(128).value; - - if exp % 2 != 0 { - // Make exponent even, so it can be divided by 2 - exp -= 1; - mant <<= 1; - } - - // Bit-by-bit (base-2 digit-by-digit) sqrt of mant. - // mant is treated here as a fixed point number with prec fractional bits. - // mant will be shifted left by one bit to have an extra fractional bit, which - // will be used to determine the rounding direction. - - // res is the truncated sqrt of mant, where one bit is added at each iteration. - let mut res = 0u128; - // rem is the remainder with the current res - // rem_i = 2^i * ((mant<<1) - res_i^2) - // starting with res = 0, rem = mant<<1 - let mut rem = mant << 1; - // s_i = 2*res_i - let mut s = 0u128; - // d is used to iterate over bits, from high to low (d_i = 2^(-i)) - let mut d = 1u128 << (prec + 1); - - // For iteration j=i+1, we need to find largest b_j = 0 or 1 such that - // (res_i + b_j * 2^(-j))^2 <= mant<<1 - // Expanding (a + b)^2 = a^2 + b^2 + 2*a*b: - // res_i^2 + (b_j * 2^(-j))^2 + 2 * res_i * b_j * 2^(-j) <= mant<<1 - // And rearranging the terms: - // b_j^2 * 2^(-j) + 2 * res_i * b_j <= 2^j * (mant<<1 - res_i^2) - // b_j^2 * 2^(-j) + 2 * res_i * b_j <= rem_i - - while d != 0 { - // Probe b_j^2 * 2^(-j) + 2 * res_i * b_j <= rem_i with b_j = 1: - // t = 2*res_i + 2^(-j) - let t = s + d; - if rem >= t { - // b_j should be 1, so make res_j = res_i + 2^(-j) and adjust rem - res += d; - s += d + d; - rem -= t; - } - // Adjust rem for next iteration - rem <<= 1; - // Shift iterator - d >>= 1; - } - - // Remove extra fractional bit from result, rounding to nearest. - // If the last bit is 0, then the nearest neighbor is definitely the lower one. - // If the last bit is 1, it sounds like this may either be a tie (if there's - // infinitely many 0s after this 1), or the nearest neighbor is the upper one. - // However, since square roots are either exact or irrational, and an exact root - // would lead to the last "extra" bit being 0, we can exclude a tie in this case. - // We therefore always round up if the last bit is 1. When the last bit is 0, - // adding 1 will not do anything since the shift will discard it. - res = (res + 1) >> 1; - - // Build resulting value with res as mantissa and exp/2 as exponent - Self::from_u128(res).value.scalbn(exp / 2 - prec) - } - } + unimplemented!() } /// IEEE-754R isSignMinus: Returns true if and only if the current value is diff --git a/src/ppc.rs b/src/ppc.rs index b32cf75..483f225 100644 --- a/src/ppc.rs +++ b/src/ppc.rs @@ -365,6 +365,10 @@ where Fallback::from(self).next_up().map(Self::from) } + fn sqrt(self) -> Self { + unimplemented!() + } + fn from_bits(input: u128) -> Self { let (a, b) = (input, input >> F::BITS); DoubleFloat(F::from_bits(a & ((1 << F::BITS) - 1)), F::from_bits(b & ((1 << F::BITS) - 1))) diff --git a/tests/ieee.rs b/tests/ieee.rs index ab1efe9..2af6e88 100644 --- a/tests/ieee.rs +++ b/tests/ieee.rs @@ -35,6 +35,43 @@ define_for_each_float_type! { rustc_apfloat::ppc::DoubleDouble, } +macro_rules! for_each_ieee_float_type { + (for<$ty_var:ident: Float> $e:expr) => {{ + { + type $ty_var = Half; + $e; + } + { + type $ty_var = Single; + $e; + } + { + type $ty_var = Double; + $e; + } + { + type $ty_var = Quad; + $e; + } + { + type $ty_var = BFloat; + $e; + } + { + type $ty_var = Float8E5M2; + $e; + } + { + type $ty_var = Float8E4M3FN; + $e; + } + { + type $ty_var = X87DoubleExtended; + $e; + } + }}; +} + trait SingleExt { fn from_f32(input: f32) -> Self; fn to_f32(self) -> f32; @@ -797,7 +834,7 @@ fn maximum() { #[test] fn sqrt() { - for_each_float_type!(for test::()); + for_each_ieee_float_type!(for test::()); fn test() { assert!(F::ZERO.sqrt().bitwise_eq(F::ZERO)); assert!((-F::ZERO).sqrt().bitwise_eq(-F::ZERO)); From 562f7efb3cf941588dca86978245465d23b333b1 Mon Sep 17 00:00:00 2001 From: Srayan Jana Date: Thu, 2 Apr 2026 19:04:17 -0700 Subject: [PATCH 13/31] add rounding mode as parameter --- src/ieee.rs | 51 +++++++++++++++++++++++++++++++++++++++------------ src/lib.rs | 3 ++- src/ppc.rs | 3 ++- tests/ieee.rs | 36 ++++++++++++++++++------------------ 4 files changed, 61 insertions(+), 32 deletions(-) diff --git a/src/ieee.rs b/src/ieee.rs index eeef2ec..90b399e 100644 --- a/src/ieee.rs +++ b/src/ieee.rs @@ -1563,7 +1563,7 @@ impl Float for IeeeFloat { } } - fn sqrt(self) -> Self { + fn sqrt(self, round: Round) -> Self { match self.category() { // preserve zero sign Category::Zero => return self, @@ -1638,18 +1638,45 @@ impl Float for IeeeFloat { d >>= 1; } - // Remove extra fractional bit from result, rounding to nearest. - // If the last bit is 0, then the nearest neighbor is definitely the lower one. - // If the last bit is 1, it sounds like this may either be a tie (if there's - // infinitely many 0s after this 1), or the nearest neighbor is the upper one. - // However, since square roots are either exact or irrational, and an exact root - // would lead to the last "extra" bit being 0, we can exclude a tie in this case. - // We therefore always round up if the last bit is 1. When the last bit is 0, - // adding 1 will not do anything since the shift will discard it. - res = (res + 1) >> 1; + let mut status = Status::OK; + + // A nonzero remainder indicates that we could continue processing sqrt if we had + // more precision, potentially indefinitely. We don't because we have enough bits + // to fill our significand already, and only need the one extra bit to determine + // rounding. + if rem != 0 { + status = Status::INEXACT; + + match round { + // If the LSB is 0, we should round down and this 1 gets cut off. If the LSB + // is 1, it is either a tie (if all remaining bits would be 0) or something + // that should be rounded up. + // + // Square roots are either exact or irrational, so a `1` in the extra bit + // already implies an irrational result with more `1`s in the infinite + // precision tail that should be rounded up, which this does. We are in a + // `rem != 0` block but could technically add the `1` unconditionally, given + // that a 0 in the extra bit would imply an exact result to be rounded down + // (and the extra bit is just shifted out). + Round::NearestTiesToEven => res += 1, + // We know we have an inexact result that needs rounding up. If the round + // bit is 1, adding 1 is sufficient and adding 2 does nothing extra (the + // new LSB will get truncated). If the round bit is 0, we need to add + // two anyway to affect the significand. + Round::TowardPositive => res += 2, + // By default, shifting will round down. + Round::TowardNegative => (), + // Same as negative since the result of sqrt is positive. + Round::TowardZero => (), + Round::NearestTiesToAway => unimplemented!("unsupported rounding mode"), + }; + } + + // Remove the extra fractional bit. + res >>= 1; - // Build resulting value with res as mantissa and exp/2 as exponent. - Self::from_u128(res).value.scalbn(exp / 2 - prec) + // Build resulting value with res as mantissa and exp/2 as exponent + status.and(Self::from_u128(res).value.scalbn(exp / 2 - prec)).value } fn from_bits(input: u128) -> Self { diff --git a/src/lib.rs b/src/lib.rs index 9a4c8c7..0467bf0 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -486,7 +486,8 @@ pub trait Float: /// IEEE-754R sqrt: Returns the correctly rounded square root of the current value /// Note: we currently don't support raising any exceptions from sqrt, so the result is always exact and the status is always OK. - fn sqrt(self) -> Self { + #[allow(unused_variables)] + fn sqrt(self, round: Round) -> Self { unimplemented!() } diff --git a/src/ppc.rs b/src/ppc.rs index 483f225..eb6872b 100644 --- a/src/ppc.rs +++ b/src/ppc.rs @@ -365,7 +365,8 @@ where Fallback::from(self).next_up().map(Self::from) } - fn sqrt(self) -> Self { + #[allow(unused_variables)] + fn sqrt(self, round: Round) -> Self { unimplemented!() } diff --git a/tests/ieee.rs b/tests/ieee.rs index 2af6e88..fc8ae3b 100644 --- a/tests/ieee.rs +++ b/tests/ieee.rs @@ -836,24 +836,24 @@ fn maximum() { fn sqrt() { for_each_ieee_float_type!(for test::()); fn test() { - assert!(F::ZERO.sqrt().bitwise_eq(F::ZERO)); - assert!((-F::ZERO).sqrt().bitwise_eq(-F::ZERO)); - assert!(F::INFINITY.sqrt().bitwise_eq(F::INFINITY)); - assert!(F::NAN.sqrt().is_nan()); - assert!((-F::INFINITY).sqrt().is_nan()); - assert!((-F::from_u128(5).value).sqrt().is_nan()); - let one = F::from_u128(1).value; - assert!(one.sqrt().bitwise_eq(one)); - let f1 = F::from_u128(64).value; - let f2 = F::from_u128(8).value; - assert!(f1.sqrt().bitwise_eq(f2)); - } - - // Cross-check against host sqrt - // TODO: This is only f32 and f64 for now since they are stable, update to other values later. - for x in 0..1000 { - assert_eq!(Single::from_u128(x).value.sqrt().to_f32(), f32::sqrt(x as f32)); - assert_eq!(Double::from_u128(x).value.sqrt().to_f64(), f64::sqrt(x as f64)); + for round in [ + Round::NearestTiesToEven, + Round::TowardPositive, + Round::TowardNegative, + Round::TowardZero, + ] { + assert!(F::ZERO.sqrt(round).bitwise_eq(F::ZERO)); + assert!((-F::ZERO).sqrt(round).bitwise_eq(-F::ZERO)); + assert!(F::INFINITY.sqrt(round).bitwise_eq(F::INFINITY)); + assert!(F::NAN.sqrt(round).is_nan()); + assert!((-F::INFINITY).sqrt(round).is_nan()); + assert!((-F::from_u128(5).value).sqrt(round).is_nan()); + let one = F::from_u128(1).value; + assert!(one.sqrt(round).bitwise_eq(one)); + let f1 = F::from_u128(64).value; + let f2 = F::from_u128(8).value; + assert!(f1.sqrt(round).bitwise_eq(f2)); + } } } From c13103d24b6eae5d2a64498259daf0113d89f2f8 Mon Sep 17 00:00:00 2001 From: Srayan Jana Date: Thu, 2 Apr 2026 19:18:44 -0700 Subject: [PATCH 14/31] move code around to downstream.rs --- src/downstream.rs | 124 +++++++++++++++++++++++++++++++++++++++++++ src/ieee.rs | 126 +++----------------------------------------- src/lib.rs | 1 + tests/downstream.rs | 66 ++++++++++++++++++++++- tests/ieee.rs | 62 ---------------------- 5 files changed, 196 insertions(+), 183 deletions(-) create mode 100644 src/downstream.rs diff --git a/src/downstream.rs b/src/downstream.rs new file mode 100644 index 0000000..9728b3f --- /dev/null +++ b/src/downstream.rs @@ -0,0 +1,124 @@ +use crate::{ + ieee::{IeeeDefaultExceptionHandling, IeeeFloat, Semantics}, + Category, Float, Round, Status, +}; + +impl IeeeFloat { + /// This is a spec conformant implementation of the IEEE Float sqrt function + /// This is put in downstream.rs because this function hasn't been implemented in the upstream C++ version yet. + pub(crate) fn ieee_sqrt(self, round: Round) -> Self { + match self.category() { + // preserve zero sign + Category::Zero => return self, + // propagate NaN + // If the input is a signalling NaN, then IEEE 754 requires the result to be converted to a quiet NaN. + // On most CPUs that means the most significant bit of the significand field is 0 for signalling NaNs and 1 for quiet NaNs. + // On most CPUs they quiet a NaN by setting that bit to a 1, RISC-V instead returns the canonical NaN with positive sign, + // the most significant significand bit set and all other significand bits cleared. + // However, Rust and LLVM allow input NaNs to be returned unmodified as well as a few other options -- see Rust's rules for NaNs. + // https://doc.rust-lang.org/std/primitive.f32.html#nan-bit-patterns + // (Thanks @programmerjake for the comment) + Category::NaN => return IeeeDefaultExceptionHandling::result_from_nan(self).value, + // sqrt of negative number is NaN + _ if self.is_negative() => return Self::NAN, + // sqrt(inf) = inf + Category::Infinity => return Self::INFINITY, + Category::Normal => (), + } + + // Floating point precision, excluding the integer bit. + let prec = i32::try_from(Self::PRECISION).unwrap() - 1; + + // x = 2^(exp - prec) * mant + // where mant is an integer with prec+1 bits. + // mant is a u128, which is large enough for the largest prec (112 for f128). + let mut exp = self.ilogb(); + let mut mant = self.scalbn(prec - exp).to_u128(128).value; + + if exp % 2 != 0 { + // Make exponent even, so it can be divided by 2. + exp -= 1; + mant <<= 1; + } + + // Bit-by-bit (base-2 digit-by-digit) sqrt of mant. + // mant is treated here as a fixed point number with prec fractional bits. + // mant will be shifted left by one bit to have an extra fractional bit, which + // will be used to determine the rounding direction. + + // res is the truncated sqrt of mant, where one bit is added at each iteration. + let mut res = 0u128; + // rem is the remainder with the current res + // rem_i = 2^i * ((mant<<1) - res_i^2) + // starting with res = 0, rem = mant<<1 + let mut rem = mant << 1; + // s_i = 2*res_i + let mut s = 0u128; + // d is used to iterate over bits, from high to low (d_i = 2^(-i)) + let mut d = 1u128 << (prec + 1); + + // For iteration j=i+1, we need to find largest b_j = 0 or 1 such that + // (res_i + b_j * 2^(-j))^2 <= mant<<1 + // Expanding (a + b)^2 = a^2 + b^2 + 2*a*b: + // res_i^2 + (b_j * 2^(-j))^2 + 2 * res_i * b_j * 2^(-j) <= mant<<1 + // And rearranging the terms: + // b_j^2 * 2^(-j) + 2 * res_i * b_j <= 2^j * (mant<<1 - res_i^2) + // b_j^2 * 2^(-j) + 2 * res_i * b_j <= rem_i + + while d != 0 { + // Probe b_j^2 * 2^(-j) + 2 * res_i * b_j <= rem_i with b_j = 1: + // t = 2*res_i + 2^(-j) + let t = s + d; + if rem >= t { + // b_j should be 1, so make res_j = res_i + 2^(-j) and adjust rem + res += d; + s += d + d; + rem -= t; + } + // Adjust rem for next iteration + rem <<= 1; + // Shift iterator + d >>= 1; + } + + let mut status = Status::OK; + + // A nonzero remainder indicates that we could continue processing sqrt if we had + // more precision, potentially indefinitely. We don't because we have enough bits + // to fill our significand already, and only need the one extra bit to determine + // rounding. + if rem != 0 { + status = Status::INEXACT; + + match round { + // If the LSB is 0, we should round down and this 1 gets cut off. If the LSB + // is 1, it is either a tie (if all remaining bits would be 0) or something + // that should be rounded up. + // + // Square roots are either exact or irrational, so a `1` in the extra bit + // already implies an irrational result with more `1`s in the infinite + // precision tail that should be rounded up, which this does. We are in a + // `rem != 0` block but could technically add the `1` unconditionally, given + // that a 0 in the extra bit would imply an exact result to be rounded down + // (and the extra bit is just shifted out). + Round::NearestTiesToEven => res += 1, + // We know we have an inexact result that needs rounding up. If the round + // bit is 1, adding 1 is sufficient and adding 2 does nothing extra (the + // new LSB will get truncated). If the round bit is 0, we need to add + // two anyway to affect the significand. + Round::TowardPositive => res += 2, + // By default, shifting will round down. + Round::TowardNegative => (), + // Same as negative since the result of sqrt is positive. + Round::TowardZero => (), + Round::NearestTiesToAway => unimplemented!("unsupported rounding mode"), + }; + } + + // Remove the extra fractional bit. + res >>= 1; + + // Build resulting value with res as mantissa and exp/2 as exponent + status.and(Self::from_u128(res).value.scalbn(exp / 2 - prec)).value + } +} diff --git a/src/ieee.rs b/src/ieee.rs index 90b399e..f0203a8 100644 --- a/src/ieee.rs +++ b/src/ieee.rs @@ -824,9 +824,9 @@ impl fmt::Debug for IeeeFloat { // but it's a bit too long to keep repeating in the Rust port for all ops. // FIXME(eddyb) find a better name/organization for all of this functionality // (`IeeeDefaultExceptionHandling` doesn't have a counterpart in the C++ code). -struct IeeeDefaultExceptionHandling; +pub(crate) struct IeeeDefaultExceptionHandling; impl IeeeDefaultExceptionHandling { - fn result_from_nan(mut r: IeeeFloat) -> StatusAnd> { + pub fn result_from_nan(mut r: IeeeFloat) -> StatusAnd> { assert!(r.is_nan()); let status = if r.is_signaling() { @@ -865,7 +865,7 @@ impl IeeeDefaultExceptionHandling { status.and(r) } - fn binop_result_from_either_nan(a: IeeeFloat, b: IeeeFloat) -> StatusAnd> { + pub fn binop_result_from_either_nan(a: IeeeFloat, b: IeeeFloat) -> StatusAnd> { let r = match (a.category(), b.category()) { (Category::NaN, _) => a, (_, Category::NaN) => b, @@ -1563,122 +1563,6 @@ impl Float for IeeeFloat { } } - fn sqrt(self, round: Round) -> Self { - match self.category() { - // preserve zero sign - Category::Zero => return self, - // propagate NaN - // If the input is a signalling NaN, then IEEE 754 requires the result to be converted to a quiet NaN. - // On most CPUs that means the most significant bit of the significand field is 0 for signalling NaNs and 1 for quiet NaNs. - // On most CPUs they quiet a NaN by setting that bit to a 1, RISC-V instead returns the canonical NaN with positive sign, - // the most significant significand bit set and all other significand bits cleared. - // However, Rust and LLVM allow input NaNs to be returned unmodified as well as a few other options -- see Rust's rules for NaNs. - // https://doc.rust-lang.org/std/primitive.f32.html#nan-bit-patterns - // (Thanks @programmerjake for the comment) - Category::NaN => return IeeeDefaultExceptionHandling::result_from_nan(self).value, - // sqrt of negative number is NaN - _ if self.is_negative() => return Self::NAN, - // sqrt(inf) = inf - Category::Infinity => return Self::INFINITY, - Category::Normal => (), - } - - // Floating point precision, excluding the integer bit. - let prec = i32::try_from(Self::PRECISION).unwrap() - 1; - - // x = 2^(exp - prec) * mant - // where mant is an integer with prec+1 bits. - // mant is a u128, which is large enough for the largest prec (112 for f128). - let mut exp = self.ilogb(); - let mut mant = self.scalbn(prec - exp).to_u128(128).value; - - if exp % 2 != 0 { - // Make exponent even, so it can be divided by 2. - exp -= 1; - mant <<= 1; - } - - // Bit-by-bit (base-2 digit-by-digit) sqrt of mant. - // mant is treated here as a fixed point number with prec fractional bits. - // mant will be shifted left by one bit to have an extra fractional bit, which - // will be used to determine the rounding direction. - - // res is the truncated sqrt of mant, where one bit is added at each iteration. - let mut res = 0u128; - // rem is the remainder with the current res - // rem_i = 2^i * ((mant<<1) - res_i^2) - // starting with res = 0, rem = mant<<1 - let mut rem = mant << 1; - // s_i = 2*res_i - let mut s = 0u128; - // d is used to iterate over bits, from high to low (d_i = 2^(-i)) - let mut d = 1u128 << (prec + 1); - - // For iteration j=i+1, we need to find largest b_j = 0 or 1 such that - // (res_i + b_j * 2^(-j))^2 <= mant<<1 - // Expanding (a + b)^2 = a^2 + b^2 + 2*a*b: - // res_i^2 + (b_j * 2^(-j))^2 + 2 * res_i * b_j * 2^(-j) <= mant<<1 - // And rearranging the terms: - // b_j^2 * 2^(-j) + 2 * res_i * b_j <= 2^j * (mant<<1 - res_i^2) - // b_j^2 * 2^(-j) + 2 * res_i * b_j <= rem_i - - while d != 0 { - // Probe b_j^2 * 2^(-j) + 2 * res_i * b_j <= rem_i with b_j = 1: - // t = 2*res_i + 2^(-j) - let t = s + d; - if rem >= t { - // b_j should be 1, so make res_j = res_i + 2^(-j) and adjust rem - res += d; - s += d + d; - rem -= t; - } - // Adjust rem for next iteration - rem <<= 1; - // Shift iterator - d >>= 1; - } - - let mut status = Status::OK; - - // A nonzero remainder indicates that we could continue processing sqrt if we had - // more precision, potentially indefinitely. We don't because we have enough bits - // to fill our significand already, and only need the one extra bit to determine - // rounding. - if rem != 0 { - status = Status::INEXACT; - - match round { - // If the LSB is 0, we should round down and this 1 gets cut off. If the LSB - // is 1, it is either a tie (if all remaining bits would be 0) or something - // that should be rounded up. - // - // Square roots are either exact or irrational, so a `1` in the extra bit - // already implies an irrational result with more `1`s in the infinite - // precision tail that should be rounded up, which this does. We are in a - // `rem != 0` block but could technically add the `1` unconditionally, given - // that a 0 in the extra bit would imply an exact result to be rounded down - // (and the extra bit is just shifted out). - Round::NearestTiesToEven => res += 1, - // We know we have an inexact result that needs rounding up. If the round - // bit is 1, adding 1 is sufficient and adding 2 does nothing extra (the - // new LSB will get truncated). If the round bit is 0, we need to add - // two anyway to affect the significand. - Round::TowardPositive => res += 2, - // By default, shifting will round down. - Round::TowardNegative => (), - // Same as negative since the result of sqrt is positive. - Round::TowardZero => (), - Round::NearestTiesToAway => unimplemented!("unsupported rounding mode"), - }; - } - - // Remove the extra fractional bit. - res >>= 1; - - // Build resulting value with res as mantissa and exp/2 as exponent - status.and(Self::from_u128(res).value.scalbn(exp / 2 - prec)).value - } - fn from_bits(input: u128) -> Self { // Dispatch to semantics. S::from_bits(input) @@ -2008,6 +1892,10 @@ impl Float for IeeeFloat { } self.scalbn_r(-*exp, round) } + + fn sqrt(self, round: Round) -> Self { + self.ieee_sqrt(round) + } } impl FloatConvert> for IeeeFloat { diff --git a/src/lib.rs b/src/lib.rs index 0467bf0..74b3724 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -762,5 +762,6 @@ macro_rules! float_common_impls { }; } +pub mod downstream; pub mod ieee; pub mod ppc; diff --git a/tests/downstream.rs b/tests/downstream.rs index fc60a7c..fcb5a8c 100644 --- a/tests/downstream.rs +++ b/tests/downstream.rs @@ -1,7 +1,7 @@ //! Tests added to `rustc_apfloat`, that were not ported from the C++ code. -use rustc_apfloat::ieee::{Double, Single, X87DoubleExtended}; -use rustc_apfloat::Float; +use rustc_apfloat::ieee::{BFloat, Double, Float8E4M3FN, Float8E5M2, Half, Quad, Single, X87DoubleExtended}; +use rustc_apfloat::{Float, Round}; // `f32 -> i128 -> f32` previously-crashing bit-patterns (found by fuzzing). pub const FUZZ_IEEE32_ROUNDTRIP_THROUGH_I128_CASES: &[u32] = &[ @@ -408,3 +408,65 @@ fn fuzz_x87_f80_neg_with_expected_outputs() { assert_eq!((-X87DoubleExtended::from_bits(bits)).to_bits(), expected_bits); } } + +macro_rules! for_each_ieee_float_type { + (for<$ty_var:ident: Float> $e:expr) => {{ + { + type $ty_var = Half; + $e; + } + { + type $ty_var = Single; + $e; + } + { + type $ty_var = Double; + $e; + } + { + type $ty_var = Quad; + $e; + } + { + type $ty_var = BFloat; + $e; + } + { + type $ty_var = Float8E5M2; + $e; + } + { + type $ty_var = Float8E4M3FN; + $e; + } + { + type $ty_var = X87DoubleExtended; + $e; + } + }}; +} + +#[test] +fn sqrt() { + for_each_ieee_float_type!(for test::()); + fn test() { + for round in [ + Round::NearestTiesToEven, + Round::TowardPositive, + Round::TowardNegative, + Round::TowardZero, + ] { + assert!(F::ZERO.sqrt(round).bitwise_eq(F::ZERO)); + assert!((-F::ZERO).sqrt(round).bitwise_eq(-F::ZERO)); + assert!(F::INFINITY.sqrt(round).bitwise_eq(F::INFINITY)); + assert!(F::NAN.sqrt(round).is_nan()); + assert!((-F::INFINITY).sqrt(round).is_nan()); + assert!((-F::from_u128(5).value).sqrt(round).is_nan()); + let one = F::from_u128(1).value; + assert!(one.sqrt(round).bitwise_eq(one)); + let f1 = F::from_u128(64).value; + let f2 = F::from_u128(8).value; + assert!(f1.sqrt(round).bitwise_eq(f2)); + } + } +} diff --git a/tests/ieee.rs b/tests/ieee.rs index fc8ae3b..0c356d9 100644 --- a/tests/ieee.rs +++ b/tests/ieee.rs @@ -35,43 +35,6 @@ define_for_each_float_type! { rustc_apfloat::ppc::DoubleDouble, } -macro_rules! for_each_ieee_float_type { - (for<$ty_var:ident: Float> $e:expr) => {{ - { - type $ty_var = Half; - $e; - } - { - type $ty_var = Single; - $e; - } - { - type $ty_var = Double; - $e; - } - { - type $ty_var = Quad; - $e; - } - { - type $ty_var = BFloat; - $e; - } - { - type $ty_var = Float8E5M2; - $e; - } - { - type $ty_var = Float8E4M3FN; - $e; - } - { - type $ty_var = X87DoubleExtended; - $e; - } - }}; -} - trait SingleExt { fn from_f32(input: f32) -> Self; fn to_f32(self) -> f32; @@ -832,31 +795,6 @@ fn maximum() { assert!(nan.maximum(f1).to_f64().is_nan()); } -#[test] -fn sqrt() { - for_each_ieee_float_type!(for test::()); - fn test() { - for round in [ - Round::NearestTiesToEven, - Round::TowardPositive, - Round::TowardNegative, - Round::TowardZero, - ] { - assert!(F::ZERO.sqrt(round).bitwise_eq(F::ZERO)); - assert!((-F::ZERO).sqrt(round).bitwise_eq(-F::ZERO)); - assert!(F::INFINITY.sqrt(round).bitwise_eq(F::INFINITY)); - assert!(F::NAN.sqrt(round).is_nan()); - assert!((-F::INFINITY).sqrt(round).is_nan()); - assert!((-F::from_u128(5).value).sqrt(round).is_nan()); - let one = F::from_u128(1).value; - assert!(one.sqrt(round).bitwise_eq(one)); - let f1 = F::from_u128(64).value; - let f2 = F::from_u128(8).value; - assert!(f1.sqrt(round).bitwise_eq(f2)); - } - } -} - #[test] fn denormal() { // Test single precision From d0a2760be78e1921f0ed889b54339bd82b251278 Mon Sep 17 00:00:00 2001 From: Srayan Jana Date: Thu, 2 Apr 2026 20:05:08 -0700 Subject: [PATCH 15/31] nan improvement --- src/downstream.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/downstream.rs b/src/downstream.rs index 9728b3f..f428b5f 100644 --- a/src/downstream.rs +++ b/src/downstream.rs @@ -18,7 +18,7 @@ impl IeeeFloat { // However, Rust and LLVM allow input NaNs to be returned unmodified as well as a few other options -- see Rust's rules for NaNs. // https://doc.rust-lang.org/std/primitive.f32.html#nan-bit-patterns // (Thanks @programmerjake for the comment) - Category::NaN => return IeeeDefaultExceptionHandling::result_from_nan(self).value, + Category::NaN => return IeeeDefaultExceptionHandling::result_from_nan(self).value.copy_sign(Self::NAN), // sqrt of negative number is NaN _ if self.is_negative() => return Self::NAN, // sqrt(inf) = inf From d645d5f9f41898b501477de6565cdfb2a47a0f1c Mon Sep 17 00:00:00 2001 From: Srayan Jana Date: Thu, 2 Apr 2026 20:21:11 -0700 Subject: [PATCH 16/31] gonna push this and hope it fuzzes correctly --- src/downstream.rs | 14 +- src/ieee.rs | 2 +- src/lib.rs | 2 +- src/ppc.rs | 2 +- tests/downstream.rs | 357 ++++++++++++++++++++++++++++++++++++++++++-- 5 files changed, 358 insertions(+), 19 deletions(-) diff --git a/src/downstream.rs b/src/downstream.rs index f428b5f..adb3e88 100644 --- a/src/downstream.rs +++ b/src/downstream.rs @@ -1,15 +1,15 @@ use crate::{ ieee::{IeeeDefaultExceptionHandling, IeeeFloat, Semantics}, - Category, Float, Round, Status, + Category, Float, Round, Status, StatusAnd, }; impl IeeeFloat { /// This is a spec conformant implementation of the IEEE Float sqrt function /// This is put in downstream.rs because this function hasn't been implemented in the upstream C++ version yet. - pub(crate) fn ieee_sqrt(self, round: Round) -> Self { + pub(crate) fn ieee_sqrt(self, round: Round) -> StatusAnd { match self.category() { // preserve zero sign - Category::Zero => return self, + Category::Zero => return Status::OK.and(self), // propagate NaN // If the input is a signalling NaN, then IEEE 754 requires the result to be converted to a quiet NaN. // On most CPUs that means the most significant bit of the significand field is 0 for signalling NaNs and 1 for quiet NaNs. @@ -18,11 +18,11 @@ impl IeeeFloat { // However, Rust and LLVM allow input NaNs to be returned unmodified as well as a few other options -- see Rust's rules for NaNs. // https://doc.rust-lang.org/std/primitive.f32.html#nan-bit-patterns // (Thanks @programmerjake for the comment) - Category::NaN => return IeeeDefaultExceptionHandling::result_from_nan(self).value.copy_sign(Self::NAN), + Category::NaN => return IeeeDefaultExceptionHandling::result_from_nan(self.copy_sign(Self::NAN)), // sqrt of negative number is NaN - _ if self.is_negative() => return Self::NAN, + _ if self.is_negative() => return Status::INVALID_OP.and(Self::NAN), // sqrt(inf) = inf - Category::Infinity => return Self::INFINITY, + Category::Infinity => return Status::OK.and(Self::INFINITY), Category::Normal => (), } @@ -119,6 +119,6 @@ impl IeeeFloat { res >>= 1; // Build resulting value with res as mantissa and exp/2 as exponent - status.and(Self::from_u128(res).value.scalbn(exp / 2 - prec)).value + status.and(Self::from_u128(res).value.scalbn(exp / 2 - prec)) } } diff --git a/src/ieee.rs b/src/ieee.rs index f0203a8..57cf32e 100644 --- a/src/ieee.rs +++ b/src/ieee.rs @@ -1893,7 +1893,7 @@ impl Float for IeeeFloat { self.scalbn_r(-*exp, round) } - fn sqrt(self, round: Round) -> Self { + fn sqrt(self, round: Round) -> StatusAnd { self.ieee_sqrt(round) } } diff --git a/src/lib.rs b/src/lib.rs index 74b3724..fa48c52 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -487,7 +487,7 @@ pub trait Float: /// IEEE-754R sqrt: Returns the correctly rounded square root of the current value /// Note: we currently don't support raising any exceptions from sqrt, so the result is always exact and the status is always OK. #[allow(unused_variables)] - fn sqrt(self, round: Round) -> Self { + fn sqrt(self, round: Round) -> StatusAnd { unimplemented!() } diff --git a/src/ppc.rs b/src/ppc.rs index eb6872b..22d0d33 100644 --- a/src/ppc.rs +++ b/src/ppc.rs @@ -366,7 +366,7 @@ where } #[allow(unused_variables)] - fn sqrt(self, round: Round) -> Self { + fn sqrt(self, round: Round) -> StatusAnd { unimplemented!() } diff --git a/tests/downstream.rs b/tests/downstream.rs index fcb5a8c..7ab88d4 100644 --- a/tests/downstream.rs +++ b/tests/downstream.rs @@ -1,7 +1,7 @@ //! Tests added to `rustc_apfloat`, that were not ported from the C++ code. use rustc_apfloat::ieee::{BFloat, Double, Float8E4M3FN, Float8E5M2, Half, Quad, Single, X87DoubleExtended}; -use rustc_apfloat::{Float, Round}; +use rustc_apfloat::{Float, Round, Status, StatusAnd}; // `f32 -> i128 -> f32` previously-crashing bit-patterns (found by fuzzing). pub const FUZZ_IEEE32_ROUNDTRIP_THROUGH_I128_CASES: &[u32] = &[ @@ -456,17 +456,356 @@ fn sqrt() { Round::TowardNegative, Round::TowardZero, ] { - assert!(F::ZERO.sqrt(round).bitwise_eq(F::ZERO)); - assert!((-F::ZERO).sqrt(round).bitwise_eq(-F::ZERO)); - assert!(F::INFINITY.sqrt(round).bitwise_eq(F::INFINITY)); - assert!(F::NAN.sqrt(round).is_nan()); - assert!((-F::INFINITY).sqrt(round).is_nan()); - assert!((-F::from_u128(5).value).sqrt(round).is_nan()); + assert!(F::ZERO.sqrt(round).value.bitwise_eq(F::ZERO)); + assert!((-F::ZERO).sqrt(round).value.bitwise_eq(-F::ZERO)); + assert!(F::INFINITY.sqrt(round).value.bitwise_eq(F::INFINITY)); + assert!(F::NAN.sqrt(round).value.is_nan()); + assert!((-F::INFINITY).sqrt(round).value.is_nan()); + assert!((-F::from_u128(5).value).sqrt(round).value.is_nan()); let one = F::from_u128(1).value; - assert!(one.sqrt(round).bitwise_eq(one)); + assert!(one.sqrt(round).value.bitwise_eq(one)); let f1 = F::from_u128(64).value; let f2 = F::from_u128(8).value; - assert!(f1.sqrt(round).bitwise_eq(f2)); + assert!(f1.sqrt(round).value.bitwise_eq(f2)); } } } + +#[test] +fn fuzz_sqrt() { + // for round in [ + // Round::NearestTiesToEven, + // Round::TowardPositive, + // Round::TowardNegative, + // Round::TowardZero, + // ] { + // println!("checking rounding mode {round:?}"); + + // for xi in 0..=u16::MAX { + // if xi % 1_000 == 0 { + // println!("{xi}/{}", u16::MAX); + // } + // let x = f16::from_bits(xi); + // let a = apfloat_sqrtf16(x, round); + // let b = hardware_sqrt::sqrtf16(x, round); + + // // x86 preserves the sign on NaN inputs, others may not (e.g. aarch64 does not). + // let eq_value = if cfg!(target_arch = "x86_64") { + // a.value.bitwise_eq(b.value) + // } else { + // a.value.bitwise_eq(b.value) || (a.value.is_nan() && b.value.is_nan()) + // }; + + // if a.status != b.status || !eq_value { + // panic!( + // "\ + // incorrect result\n\ + // xi: {xi:#010x}\n\ + // x: {x:?}\n\ + // a: {a:?} {af}\n\ + // b: {b:?} {bf}\ + // ", + // af = f32::from_bits(a.value.to_bits().try_into().unwrap()), + // bf = f32::from_bits(a.value.to_bits().try_into().unwrap()), + // ); + // } + // } + // } + + for round in [ + Round::NearestTiesToEven, + Round::TowardPositive, + Round::TowardNegative, + Round::TowardZero, + ] { + println!("checking rounding mode {round:?}"); + + for xi in 0..=u32::MAX { + if xi % 100_000_000 == 0 { + println!("{xi}/{}", u32::MAX); + } + let x = f32::from_bits(xi); + let a = apfloat_sqrt(x, round); + let b = hardware_sqrt::sqrtf32(x, round); + + // x86 preserves the sign on NaN inputs, others may not (e.g. aarch64 does not). + let eq_value = if cfg!(target_arch = "x86_64") { + a.value.bitwise_eq(b.value) + } else { + a.value.bitwise_eq(b.value) || (a.value.is_nan() && b.value.is_nan()) + }; + + if a.status != b.status || !eq_value { + panic!( + "\ + incorrect result\n\ + xi: {xi:#010x}\n\ + x: {x:?}\n\ + a: {a:?} {af}\n\ + b: {b:?} {bf}\ + ", + af = f32::from_bits(a.value.to_bits().try_into().unwrap()), + bf = f32::from_bits(a.value.to_bits().try_into().unwrap()), + ); + } + } + } +} + +#[cfg(target_arch = "aarch64")] +fn apfloat_sqrtf16(x: f16, round: Round) -> StatusAnd { + Half::from_bits(x.to_bits().into()).sqrt(round) +} + +fn apfloat_sqrt(x: f32, round: Round) -> StatusAnd { + Single::from_bits(x.to_bits().into()).sqrt(round) +} + +// SQRTSS is in baseline SSE, SQRTSD needs SSE2 +#[cfg(any(target_arch = "x86_64", all(target_arch = "x86", target_feature = "sse2")))] +mod hardware_sqrt { + use super::*; + + pub fn sqrtf32(mut x: f32, round: Round) -> StatusAnd { + let mut csr_stash = 0u32; + let mut csr = make_mxcsr_cw(round); + + unsafe { + core::arch::asm!( + // stash the current control state + "stmxcsr [{csr_stash}]", + // set the control state we want, clears flags + "ldmxcsr [{csr}]", + // run sqrt + "sqrtss {x}, {x}", + // get the new control state + "stmxcsr [{csr}]", + // restore the original control state + "ldmxcsr [{csr_stash}]", + csr_stash = in(reg) &mut csr_stash, + csr = in(reg) &mut csr, + x = inout(xmm_reg) x, + options(nostack), + ); + } + + let status = check_exceptions(csr); + status.and(Single::from_bits(x.to_bits().into())) + } + + pub fn sqrtf64(mut x: f64, round: Round) -> StatusAnd { + let mut csr_stash = 0u32; + let mut csr = make_mxcsr_cw(round); + + unsafe { + core::arch::asm!( + // stash the current control state + "stmxcsr [{csr_stash}]", + // set the control state we want, clears flags + "ldmxcsr [{csr}]", + // run sqrt + "sqrtsd {x}, {x}", + // get the new control state + "stmxcsr [{csr}]", + // restore the original control state + "ldmxcsr [{csr_stash}]", + csr_stash = in(reg) &mut csr_stash, + csr = in(reg) &mut csr, + x = inout(xmm_reg) x, + options(nostack), + ); + } + + let status = check_exceptions(csr); + status.and(Double::from_bits(x.to_bits().into())) + } + + fn make_mxcsr_cw(round: Round) -> u32 { + // Default: Clear exception flags, no DAZ, no FTZ + let mut csr = 0u32; + // Set all masks so fp status doesn't turn into SIGFPE + csr |= 0b111111 << 7; + + let rc = match round { + Round::NearestTiesToEven => 0b00, + Round::TowardPositive => 0b10, + Round::TowardNegative => 0b01, + Round::TowardZero => 0b11, + Round::NearestTiesToAway => unimplemented!("unsupported rounding on x86"), + }; + + csr |= rc << 13; + csr + } + + fn check_exceptions(csr: u32) -> Status { + let mut status = Status::OK; + + if csr & (1 << 0) != 0 { + status |= Status::INVALID_OP; + } + if csr & (1 << 1) != 0 { + // denormal flag, not part of status + } + if csr & (1 << 2) != 0 { + status |= Status::DIV_BY_ZERO; + } + if csr & (1 << 3) != 0 { + status |= Status::OVERFLOW; + } + if csr & (1 << 4) != 0 { + status |= Status::UNDERFLOW; + } + if csr & (1 << 5) != 0 { + status |= Status::INEXACT; + } + + status + } +} + +#[cfg(target_arch = "aarch64")] +mod hardware_sqrt { + use super::*; + + #[cfg(target_feature = "fp16")] + pub fn sqrtf16(mut x: f16, round: Round) -> StatusAnd { + let fpcr = make_fpcr_cw(round); + let fpsr: u64; + + unsafe { + core::arch::asm!( + // stash the current control and status state + "mrs {fpcr_stash}, fpcr", + "mrs {fpsr_stash}, fpsr", + // zero the exception flags, set desired control state + "and {tmp}, {fpsr_stash}, #{CLEAR}", + "msr fpsr, {tmp}", + "msr fpcr, {fpcr}", + // run sqrt + "fsqrt {x:h}, {x:h}", + // get the status word + "mrs {fpcr}, fpsr", + // restore the original control and exception state + "msr fpcr, {fpcr_stash}", + "msr fpsr, {fpsr_stash}", + fpsr_stash = out(reg) _, + fpcr_stash = out(reg) _, + tmp = out(reg) _, + fpcr = inout(reg) fpcr => fpsr, + x = inout(vreg) x, + CLEAR = const !0x1f_u64, // + options(nomem, nostack), + ); + } + let status = check_exceptions(fpsr); + status.and(Half::from_bits(x.to_bits().into())) + } + + pub fn sqrtf32(mut x: f32, round: Round) -> StatusAnd { + let fpcr = make_fpcr_cw(round); + let fpsr: u64; + + unsafe { + core::arch::asm!( + // stash the current control and status state + "mrs {fpcr_stash}, fpcr", + "mrs {fpsr_stash}, fpsr", + // zero the exception flags, set desired control state + "and {tmp}, {fpsr_stash}, #{CLEAR}", + "msr fpsr, {tmp}", + "msr fpcr, {fpcr}", + // run sqrt + "fsqrt {x:s}, {x:s}", + // get the status word + "mrs {fpcr}, fpsr", + // restore the original control and exception state + "msr fpcr, {fpcr_stash}", + "msr fpsr, {fpsr_stash}", + fpsr_stash = out(reg) _, + fpcr_stash = out(reg) _, + tmp = out(reg) _, + fpcr = inout(reg) fpcr => fpsr, + x = inout(vreg) x, + CLEAR = const !0x1f_u64, // + options(nomem, nostack), + ); + } + let status = check_exceptions(fpsr); + status.and(Single::from_bits(x.to_bits().into())) + } + + pub fn sqrtf64(mut x: f64, round: Round) -> StatusAnd { + let fpcr = make_fpcr_cw(round); + let fpsr: u64; + + unsafe { + core::arch::asm!( + // stash the current control and status state + "mrs {fpcr_stash}, fpcr", + "mrs {fpsr_stash}, fpsr", + // zero the exception flags, set desired control state + "and {tmp}, {fpsr_stash}, #{CLEAR}", + "msr fpsr, {tmp}", + "msr fpcr, {fpcr}", + // run sqrt + "fsqrt {x:d}, {x:d}", + // get the status word + "mrs {fpcr}, fpsr", + // restore the original control and exception state + "msr fpcr, {fpcr_stash}", + "msr fpsr, {fpsr_stash}", + fpsr_stash = out(reg) _, + fpcr_stash = out(reg) _, + tmp = out(reg) _, + fpcr = inout(reg) fpcr => fpsr, + x = inout(vreg) x, + CLEAR = const !0x1f_u64, // + options(nomem, nostack), + ); + } + let status = check_exceptions(fpsr); + status.and(Double::from_bits(x.to_bits().into())) + } + + fn make_fpcr_cw(round: Round) -> u64 { + // Default: Clear exception flags, no DAZ, no FTZ + let mut csr = 0u64; + + // Disable traps on all 5 floating point exceptions + csr |= 0b11111 << 8; + + let rc = match round { + Round::NearestTiesToEven => 0b00, + Round::TowardPositive => 0b01, + Round::TowardNegative => 0b10, + Round::TowardZero => 0b11, + Round::NearestTiesToAway => unimplemented!("unsupported rounding on aarch64"), + }; + + csr |= rc << 22; + csr + } + + fn check_exceptions(fpcr: u64) -> Status { + let mut status = Status::OK; + + if fpcr & (1 << 0) != 0 { + status |= Status::INVALID_OP; + } + if fpcr & (1 << 1) != 0 { + status |= Status::DIV_BY_ZERO; + } + if fpcr & (1 << 2) != 0 { + status |= Status::OVERFLOW; + } + if fpcr & (1 << 3) != 0 { + status |= Status::UNDERFLOW; + } + if fpcr & (1 << 4) != 0 { + status |= Status::INEXACT; + } + + status + } +} From fa1f2607f029299a13c13e74cbe244adb87b69e9 Mon Sep 17 00:00:00 2001 From: Srayan Jana Date: Thu, 2 Apr 2026 20:59:14 -0700 Subject: [PATCH 17/31] add dead code check --- tests/downstream.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/downstream.rs b/tests/downstream.rs index 7ab88d4..fb67077 100644 --- a/tests/downstream.rs +++ b/tests/downstream.rs @@ -593,6 +593,7 @@ mod hardware_sqrt { status.and(Single::from_bits(x.to_bits().into())) } + #[allow(dead_code)] pub fn sqrtf64(mut x: f64, round: Round) -> StatusAnd { let mut csr_stash = 0u32; let mut csr = make_mxcsr_cw(round); From 23db47a1929f27eb41d16820f28433aaa35fe8f7 Mon Sep 17 00:00:00 2001 From: Srayan Jana Date: Thu, 2 Apr 2026 22:49:45 -0700 Subject: [PATCH 18/31] remove copy sign, hope this fuzzes correctly --- src/downstream.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/downstream.rs b/src/downstream.rs index adb3e88..f63595c 100644 --- a/src/downstream.rs +++ b/src/downstream.rs @@ -18,7 +18,7 @@ impl IeeeFloat { // However, Rust and LLVM allow input NaNs to be returned unmodified as well as a few other options -- see Rust's rules for NaNs. // https://doc.rust-lang.org/std/primitive.f32.html#nan-bit-patterns // (Thanks @programmerjake for the comment) - Category::NaN => return IeeeDefaultExceptionHandling::result_from_nan(self.copy_sign(Self::NAN)), + Category::NaN => return IeeeDefaultExceptionHandling::result_from_nan(self), // sqrt of negative number is NaN _ if self.is_negative() => return Status::INVALID_OP.and(Self::NAN), // sqrt(inf) = inf From a14e94ec6c25cbfb14833b2ebace464cd704bf40 Mon Sep 17 00:00:00 2001 From: miro Date: Wed, 7 Oct 2026 15:22:14 +0200 Subject: [PATCH 19/31] fixed compiler errors --- fuzz/src/host.rs | 15 +++++++++++++++ fuzz/src/main.rs | 2 +- 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/fuzz/src/host.rs b/fuzz/src/host.rs index 01ebe2f..7c5a65f 100644 --- a/fuzz/src/host.rs +++ b/fuzz/src/host.rs @@ -342,6 +342,11 @@ mod x86 { no_fp_env(rm, || x as Self) } } + + // TODO: write the ASM for other rounding modes. + fn sqrt(self, rm: Round) -> Option> { + no_fp_env(rm, || self.sqrt()) + } } impl HostFloat for f32 { @@ -479,6 +484,11 @@ mod x86 { fn from_single_r(x: f32, _rm: Round) -> Option> { Some(Status::OK.and(x)) } + + // TODO: write the ASM for other rounding modes. + fn sqrt(self, rm: Round) -> Option> { + no_fp_env(rm, || self.sqrt()) + } } impl HostFloat for f64 { @@ -616,6 +626,11 @@ mod x86 { Some(status.and(dst)) } } + + // TODO: write the ASM for other rounding modes. + fn sqrt(self, rm: Round) -> Option> { + no_fp_env(rm, || self.sqrt()) + } } /// Make a control word or return `None` if the rounding mode is not supported. diff --git a/fuzz/src/main.rs b/fuzz/src/main.rs index b97da44..44096c9 100644 --- a/fuzz/src/main.rs +++ b/fuzz/src/main.rs @@ -1114,13 +1114,13 @@ fn eval_host( // FIXME: rem disregards rounding mode Op::Rem => a.rem(b, rm)?, Op::MulAdd => a.mul_add_r(b, c, rm)?, + Op::Sqrt => a.sqrt(rm)?, // FIXME: the below operations discard a status. We should turn them into // unidirectional operations. Op::FToI128ToF => F::from_i128_r(a.to_i128_r(rm)?.value, rm)?, Op::FToU128ToF => F::from_u128_r(a.to_u128_r(rm)?.value, rm)?, Op::FToSingleToF => F::from_single_r(a.to_single_r(rm)?.value, rm)?, Op::FToDoubleToF => F::from_double_r(a.to_double_r(rm)?.value, rm)?, - Op::Sqrt => F::s }; Some(res.map(F::to_bits)) From f0c10783745ad0a0ec256de5fc5ee2398cd60f8d Mon Sep 17 00:00:00 2001 From: miro Date: Thu, 8 Oct 2026 00:53:26 +0200 Subject: [PATCH 20/31] made fuzz_sqrt parallel --- Cargo.toml | 1 + tests/downstream.rs | 72 +++++++++++++++++++++++++++------------------ 2 files changed, 45 insertions(+), 28 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index f6285f4..be952fc 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -31,6 +31,7 @@ smallvec = { version = "1.11.0", features = ["const_generics", "union"] } [dev-dependencies] criterion = { version = "0.5.1", features = ["html_reports"] } +rayon = "1.10.0" # Speed up some things like decoding and corpus generation [profile.dev.package.rustc_apfloat-fuzz] diff --git a/tests/downstream.rs b/tests/downstream.rs index fb67077..272257b 100644 --- a/tests/downstream.rs +++ b/tests/downstream.rs @@ -1,5 +1,8 @@ //! Tests added to `rustc_apfloat`, that were not ported from the C++ code. +use std::panic::panic_any; + +use rayon::iter::{IntoParallelIterator, ParallelIterator}; use rustc_apfloat::ieee::{BFloat, Double, Float8E4M3FN, Float8E5M2, Half, Quad, Single, X87DoubleExtended}; use rustc_apfloat::{Float, Round, Status, StatusAnd}; @@ -520,34 +523,47 @@ fn fuzz_sqrt() { ] { println!("checking rounding mode {round:?}"); - for xi in 0..=u32::MAX { - if xi % 100_000_000 == 0 { - println!("{xi}/{}", u32::MAX); - } - let x = f32::from_bits(xi); - let a = apfloat_sqrt(x, round); - let b = hardware_sqrt::sqrtf32(x, round); - - // x86 preserves the sign on NaN inputs, others may not (e.g. aarch64 does not). - let eq_value = if cfg!(target_arch = "x86_64") { - a.value.bitwise_eq(b.value) - } else { - a.value.bitwise_eq(b.value) || (a.value.is_nan() && b.value.is_nan()) - }; - - if a.status != b.status || !eq_value { - panic!( - "\ - incorrect result\n\ - xi: {xi:#010x}\n\ - x: {x:?}\n\ - a: {a:?} {af}\n\ - b: {b:?} {bf}\ - ", - af = f32::from_bits(a.value.to_bits().try_into().unwrap()), - bf = f32::from_bits(a.value.to_bits().try_into().unwrap()), - ); - } + let result = (0..u32::MAX) + .into_par_iter() + .map(|xi| { + // TODO: figure out if this still makes sense. + if xi % 100_000_000 == 0 { + println!("{xi}/{}", u32::MAX); + } + + let x = f32::from_bits(xi); + let software = apfloat_sqrt(x, round); + let hardware = hardware_sqrt::sqrtf32(x, round); + + // x86 preserves the sign on NaN inputs, others may not (e.g. aarch64 does not). + let eq_value = if cfg!(target_arch = "x86_64") { + software.value.bitwise_eq(hardware.value) + } else { + software.value.bitwise_eq(hardware.value) || (software.value.is_nan() && hardware.value.is_nan()) + }; + + if software.status != hardware.status || !eq_value { + let str = format!( + "\ + incorrect result\n\ + xi: {xi:#010x}\n\ + x: {x:?}\n\ + software: {software:?} {af}\n\ + hardware: {hardware:?} {bf}\ + ", + af = f32::from_bits(software.value.to_bits().try_into().unwrap()), + bf = f32::from_bits(software.value.to_bits().try_into().unwrap()), + ); + + Err(str) + } else { + Ok(()) + } + }) + .collect::>(); + + if let Err(err) = result { + panic!("{}", err); } } } From cfb4c01eb194e72e0641182f6e54dd9879b197f4 Mon Sep 17 00:00:00 2001 From: miro Date: Thu, 8 Oct 2026 01:23:34 +0200 Subject: [PATCH 21/31] fixed sqrt of a negative number NaNs having the wrong NaN sign --- src/downstream.rs | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/downstream.rs b/src/downstream.rs index f63595c..89fd6ff 100644 --- a/src/downstream.rs +++ b/src/downstream.rs @@ -19,8 +19,10 @@ impl IeeeFloat { // https://doc.rust-lang.org/std/primitive.f32.html#nan-bit-patterns // (Thanks @programmerjake for the comment) Category::NaN => return IeeeDefaultExceptionHandling::result_from_nan(self), - // sqrt of negative number is NaN - _ if self.is_negative() => return Status::INVALID_OP.and(Self::NAN), + // sqrt of negative number is NaN. + // + // according to fuzzing against x86_64 it's negative NaN. + _ if self.is_negative() => return Status::INVALID_OP.and(-Self::NAN), // sqrt(inf) = inf Category::Infinity => return Status::OK.and(Self::INFINITY), Category::Normal => (), From 097213276d5a4e5abeb59a0a38a269d14a811bc6 Mon Sep 17 00:00:00 2001 From: miro Date: Thu, 8 Oct 2026 15:01:02 +0200 Subject: [PATCH 22/31] undid making NaNs negative as that's just x86 behaviour, fixed fuzz_sqrt comparing NaN signs on non-NaN inputs, and fixed Op::ALL not containing all Op variants --- fuzz/src/main.rs | 1 + justfile | 4 ++-- src/downstream.rs | 2 +- tests/downstream.rs | 5 ++++- 4 files changed, 8 insertions(+), 4 deletions(-) diff --git a/fuzz/src/main.rs b/fuzz/src/main.rs index 44096c9..ab92293 100644 --- a/fuzz/src/main.rs +++ b/fuzz/src/main.rs @@ -422,6 +422,7 @@ impl Op { Self::FToU128ToF, Self::FToSingleToF, Self::FToDoubleToF, + Self::Sqrt, ]; pub fn airity(self) -> Arity { diff --git a/justfile b/justfile index 353af38..46bdcb9 100644 --- a/justfile +++ b/justfile @@ -16,8 +16,8 @@ _default: {{ just_executable() }} --list # Run non-fuzzing tests -test: - cargo test --workspace +test *names: + cargo test --workspace {{names}} # Create directories and build the executable, but don't start fuzzing. fuzz-build: diff --git a/src/downstream.rs b/src/downstream.rs index 89fd6ff..76d7509 100644 --- a/src/downstream.rs +++ b/src/downstream.rs @@ -22,7 +22,7 @@ impl IeeeFloat { // sqrt of negative number is NaN. // // according to fuzzing against x86_64 it's negative NaN. - _ if self.is_negative() => return Status::INVALID_OP.and(-Self::NAN), + _ if self.is_negative() => return Status::INVALID_OP.and(Self::NAN), // sqrt(inf) = inf Category::Infinity => return Status::OK.and(Self::INFINITY), Category::Normal => (), diff --git a/tests/downstream.rs b/tests/downstream.rs index 272257b..a1ae4fd 100644 --- a/tests/downstream.rs +++ b/tests/downstream.rs @@ -536,7 +536,10 @@ fn fuzz_sqrt() { let hardware = hardware_sqrt::sqrtf32(x, round); // x86 preserves the sign on NaN inputs, others may not (e.g. aarch64 does not). - let eq_value = if cfg!(target_arch = "x86_64") { + // + // non-NaN inputs result in a newly-created NaN if negative, + // and newly-created NaN bit patterns from x86 have a negative sign while this library has a positive sign. + let eq_value = if cfg!(target_arch = "x86_64") && x.is_nan() { software.value.bitwise_eq(hardware.value) } else { software.value.bitwise_eq(hardware.value) || (software.value.is_nan() && hardware.value.is_nan()) From bed837a89d9a3d671ac4bb2a8df21b703489c77c Mon Sep 17 00:00:00 2001 From: miro Date: Thu, 8 Oct 2026 17:49:53 +0200 Subject: [PATCH 23/31] fixed the intentional panic in ieee_sqrt counting as a crash when fuzzing and the fuzzer attempting to run the C++ implementation of sqrt (which doesn't exist) --- fuzz/src/main.rs | 115 +++++++++++++++++++++++++++++------------------ 1 file changed, 72 insertions(+), 43 deletions(-) diff --git a/fuzz/src/main.rs b/fuzz/src/main.rs index ab92293..681487d 100644 --- a/fuzz/src/main.rs +++ b/fuzz/src/main.rs @@ -10,6 +10,9 @@ mod host; use io::IsTerminal; use io::Read; use std::io; +use std::panic::AssertUnwindSafe; +use std::panic::catch_unwind; +use std::panic::resume_unwind; use std::path::PathBuf; use std::{fmt, fs}; @@ -441,6 +444,10 @@ impl Op { Self::Sqrt => Arity::Binary, } } + + fn no_cxx(self) -> bool { + self == Self::Sqrt + } } /// Number of inputs to an operation. @@ -685,8 +692,7 @@ where // Evaluate the APFloat version as well as all possible references. FuzzOpEvalOutputs { rs_apf: eval_rust_ap(cfg.op, cfg.rm, a, b, c), - cxx_apf: cfg - .run_cxx + cxx_apf: (cfg.run_cxx && !cfg.op.no_cxx()) .then(|| F::cxx_apf_eval_fuzz_op(cfg.op, cfg.rm, a, b, c)), host: cfg .run_host @@ -722,7 +728,7 @@ fn decode_operands(op: Op, data: &[u8]) -> Result<(F, F, F), Decod /// Collected outputs from all input sources. #[derive(Clone, Debug)] struct FuzzOpEvalOutputs { - rs_apf: StatusAnd, + rs_apf: Option>, cxx_apf: Option>, host: Option>, } @@ -799,26 +805,33 @@ impl FuzzOpEvalOutputs { host_stat_ignore: (!F::host_supports_fp_env()).then_some("host does not support fpenv"), }; + // If there is no rust result then there's nothing to fuzz against. + // currently this case happens if rust intentionally panics. + let Some(rs_apf) = self.rs_apf else { + println!("no rust result to fuzz against"); + return Ok(()); + }; + if let Some(cxx_res) = self.cxx_apf { - res.cxx_error = self.rs_apf.value != cxx_res.value; - res.cxx_stat_error = self.rs_apf.status != cxx_res.status; + res.cxx_error = rs_apf.value != cxx_res.value; + res.cxx_stat_error = rs_apf.status != cxx_res.status; res.cxx_ignore = res .cxx_ignore - .or_else(|| ignore_cxx(cfg, a, b, c, self.rs_apf.value, cxx_res.value)); + .or_else(|| ignore_cxx(cfg, a, b, c, rs_apf.value, cxx_res.value)); res.cxx_stat_ignore = res .cxx_stat_ignore - .or_else(|| ignore_cxx_status(cfg, a, b, c, self.rs_apf, cxx_res)); + .or_else(|| ignore_cxx_status(cfg, a, b, c, rs_apf, cxx_res)); } if let Some(host_res) = self.host { - res.host_error = self.rs_apf.value != host_res.value; - res.host_stat_error = self.rs_apf.status != host_res.status; + res.host_error = rs_apf.value != host_res.value; + res.host_stat_error = rs_apf.status != host_res.status; res.host_ignore = res .host_ignore - .or_else(|| ignore_host(cfg, a, b, c, self.rs_apf.value, host_res.value)); + .or_else(|| ignore_host(cfg, a, b, c, rs_apf.value, host_res.value)); res.host_stat_ignore = res .host_stat_ignore - .or_else(|| ignore_host_status(cfg, a, b, c, self.rs_apf, host_res)); + .or_else(|| ignore_host_status(cfg, a, b, c, rs_apf, host_res)); } let failure = (res.cxx_error && res.cxx_ignore.is_none()) @@ -840,7 +853,7 @@ impl FuzzOpEvalOutputs { println!(", {:?})", cfg.rm); print_float( - self.rs_apf, + rs_apf, "Rust / rustc_apfloat", false, false, @@ -1055,42 +1068,58 @@ fn ignore_host_status( } /// Execute the requested operation as an AP float with the given rounding mode. -fn eval_rust_ap(op: Op, rm: Round, a: F, b: F, c: F) -> StatusAnd +fn eval_rust_ap(op: Op, rm: Round, a: F, b: F, c: F) -> Option> where Single: FloatConvert<::RustcApFloat>, Double: FloatConvert<::RustcApFloat>, { - let res = match op { - Op::Neg => Status::OK.and(-a.to_ap()), - Op::Add => a.to_ap().add_r(b.to_ap(), rm), - Op::Sub => a.to_ap().sub_r(b.to_ap(), rm), - Op::Mul => a.to_ap().mul_r(b.to_ap(), rm), - Op::Div => a.to_ap().div_r(b.to_ap(), rm), - // FIXME: rem disregards rounding mode - Op::Rem => a.to_ap() % b.to_ap(), - Op::MulAdd => a.to_ap().mul_add_r(b.to_ap(), c.to_ap(), rm), - // FIXME: the below operations discard a status. We should turn them into - // unidirectional operations. - Op::FToI128ToF => { - F::RustcApFloat::from_i128_r(a.to_ap().to_i128_r(128, rm, &mut false).value, rm) - } - Op::FToU128ToF => { - F::RustcApFloat::from_u128_r(a.to_ap().to_u128_r(128, rm, &mut false).value, rm) - } - Op::FToSingleToF => FloatConvert::::convert_r( - FloatConvert::::convert_r(a.to_ap(), rm, &mut false).value, - rm, - &mut false, - ), - Op::FToDoubleToF => FloatConvert::::convert_r( - FloatConvert::::convert_r(a.to_ap(), rm, &mut false).value, - rm, - &mut false, - ), - Op::Sqrt => a.to_ap().sqrt(rm), - }; + let result = catch_unwind(AssertUnwindSafe(|| { + let res = match op { + Op::Neg => Status::OK.and(-a.to_ap()), + Op::Add => a.to_ap().add_r(b.to_ap(), rm), + Op::Sub => a.to_ap().sub_r(b.to_ap(), rm), + Op::Mul => a.to_ap().mul_r(b.to_ap(), rm), + Op::Div => a.to_ap().div_r(b.to_ap(), rm), + // FIXME: rem disregards rounding mode + Op::Rem => a.to_ap() % b.to_ap(), + Op::MulAdd => a.to_ap().mul_add_r(b.to_ap(), c.to_ap(), rm), + // FIXME: the below operations discard a status. We should turn them into + // unidirectional operations. + Op::FToI128ToF => { + F::RustcApFloat::from_i128_r(a.to_ap().to_i128_r(128, rm, &mut false).value, rm) + } + Op::FToU128ToF => { + F::RustcApFloat::from_u128_r(a.to_ap().to_u128_r(128, rm, &mut false).value, rm) + } + Op::FToSingleToF => FloatConvert::::convert_r( + FloatConvert::::convert_r(a.to_ap(), rm, &mut false).value, + rm, + &mut false, + ), + Op::FToDoubleToF => FloatConvert::::convert_r( + FloatConvert::::convert_r(a.to_ap(), rm, &mut false).value, + rm, + &mut false, + ), + Op::Sqrt => a.to_ap().sqrt(rm), + }; - res.map(F::from_ap) + res.map(F::from_ap) + })); + + match result { + Ok(ok) => Some(ok), + // HACK: this attempts to prevent ieee_sqrt's intentional panic on unsupported rounding modes + // from causing false-positive crashes. + Err(payload) if let Some(&message) = payload.downcast_ref::<&'static str>() => { + if message == "not implemented: unsupported rounding mode" { + None + } else { + resume_unwind(Box::new(message)) + } + } + Err(payload) => resume_unwind(payload), + } } /// Execute the requested operation on the host with the given rounding mode, if possible. If From 22c7757047cc516038133dbf068fb442e57497f4 Mon Sep 17 00:00:00 2001 From: miro Date: Sat, 10 Oct 2026 01:50:56 +0200 Subject: [PATCH 24/31] added correct host sqrt fuzzing using asm --- fuzz/src/host.rs | 27 ++++++++++++++++++++++++--- 1 file changed, 24 insertions(+), 3 deletions(-) diff --git a/fuzz/src/host.rs b/fuzz/src/host.rs index 7c5a65f..ffe6cdd 100644 --- a/fuzz/src/host.rs +++ b/fuzz/src/host.rs @@ -485,9 +485,19 @@ mod x86 { Some(Status::OK.and(x)) } - // TODO: write the ASM for other rounding modes. fn sqrt(self, rm: Round) -> Option> { - no_fp_env(rm, || self.sqrt()) + // WARNING(miro): this is my first time writing inline ASM. + unsafe { + let output: f32; + + let status = with_fp_env!( + rm, + "sqrtss {x}, {y}", + x = out(xmm_reg) output, + y = in(xmm_reg) self, + ); + Some(status.and(output)) + } } } @@ -629,7 +639,18 @@ mod x86 { // TODO: write the ASM for other rounding modes. fn sqrt(self, rm: Round) -> Option> { - no_fp_env(rm, || self.sqrt()) + // WARNING(miro): this is my second time writing inline ASM. + unsafe { + let output: f64; + + let status = with_fp_env!( + rm, + "sqrtsd {x}, {y}", + x = out(xmm_reg) output, + y = in(xmm_reg) self, + ); + Some(status.and(output)) + } } } From 9a8771d585d1b78224e88a0f98b00e0a8b61c992 Mon Sep 17 00:00:00 2001 From: miro Date: Sat, 10 Oct 2026 01:56:05 +0200 Subject: [PATCH 25/31] removed files I used to make the fuzzer build --- .envrc | 1 - rust-toolchain.toml | 2 -- shell.nix | 8 -------- 3 files changed, 11 deletions(-) delete mode 100644 .envrc delete mode 100644 rust-toolchain.toml delete mode 100644 shell.nix diff --git a/.envrc b/.envrc deleted file mode 100644 index 1d953f4..0000000 --- a/.envrc +++ /dev/null @@ -1 +0,0 @@ -use nix diff --git a/rust-toolchain.toml b/rust-toolchain.toml deleted file mode 100644 index 5d56faf..0000000 --- a/rust-toolchain.toml +++ /dev/null @@ -1,2 +0,0 @@ -[toolchain] -channel = "nightly" diff --git a/shell.nix b/shell.nix deleted file mode 100644 index cd1c3ac..0000000 --- a/shell.nix +++ /dev/null @@ -1,8 +0,0 @@ -{pkgs ? import {}}: -pkgs.mkShell { - buildInputs = with pkgs; [ - python3 - ]; - - LD_LIBRARY_PATH = "${pkgs.stdenv.cc.cc.lib}/lib"; -} From 4d5a35109dbb8e2bcebaa0d47874dccacb2caca2 Mon Sep 17 00:00:00 2001 From: miro Date: Sat, 10 Oct 2026 01:58:01 +0200 Subject: [PATCH 26/31] this one also needs to be deleted --- .cargo/config.toml | 2 -- 1 file changed, 2 deletions(-) delete mode 100644 .cargo/config.toml diff --git a/.cargo/config.toml b/.cargo/config.toml deleted file mode 100644 index a6b014e..0000000 --- a/.cargo/config.toml +++ /dev/null @@ -1,2 +0,0 @@ -[build] -target-dir = "target" From 7c00892a4bc15ec57ddc48191ed5c43b3cb1c4c3 Mon Sep 17 00:00:00 2001 From: miro Date: Sat, 10 Oct 2026 02:09:53 +0200 Subject: [PATCH 27/31] rustfmt --- fuzz/src/host.rs | 4 ++-- fuzz/src/main.rs | 13 +++---------- src/downstream.rs | 2 +- src/lib.rs | 2 +- tests/downstream.rs | 2 +- 5 files changed, 8 insertions(+), 15 deletions(-) diff --git a/fuzz/src/host.rs b/fuzz/src/host.rs index ffe6cdd..cf5a654 100644 --- a/fuzz/src/host.rs +++ b/fuzz/src/host.rs @@ -489,7 +489,7 @@ mod x86 { // WARNING(miro): this is my first time writing inline ASM. unsafe { let output: f32; - + let status = with_fp_env!( rm, "sqrtss {x}, {y}", @@ -642,7 +642,7 @@ mod x86 { // WARNING(miro): this is my second time writing inline ASM. unsafe { let output: f64; - + let status = with_fp_env!( rm, "sqrtsd {x}, {y}", diff --git a/fuzz/src/main.rs b/fuzz/src/main.rs index 681487d..9a0c497 100644 --- a/fuzz/src/main.rs +++ b/fuzz/src/main.rs @@ -811,7 +811,7 @@ impl FuzzOpEvalOutputs { println!("no rust result to fuzz against"); return Ok(()); }; - + if let Some(cxx_res) = self.cxx_apf { res.cxx_error = rs_apf.value != cxx_res.value; res.cxx_stat_error = rs_apf.status != cxx_res.status; @@ -852,14 +852,7 @@ impl FuzzOpEvalOutputs { } println!(", {:?})", cfg.rm); - print_float( - rs_apf, - "Rust / rustc_apfloat", - false, - false, - None, - None, - ); + print_float(rs_apf, "Rust / rustc_apfloat", false, false, None, None); if let Some(cxx_res) = self.cxx_apf { print_float( cxx_res, @@ -1106,7 +1099,7 @@ where res.map(F::from_ap) })); - + match result { Ok(ok) => Some(ok), // HACK: this attempts to prevent ieee_sqrt's intentional panic on unsupported rounding modes diff --git a/src/downstream.rs b/src/downstream.rs index 76d7509..a56f2e2 100644 --- a/src/downstream.rs +++ b/src/downstream.rs @@ -20,7 +20,7 @@ impl IeeeFloat { // (Thanks @programmerjake for the comment) Category::NaN => return IeeeDefaultExceptionHandling::result_from_nan(self), // sqrt of negative number is NaN. - // + // // according to fuzzing against x86_64 it's negative NaN. _ if self.is_negative() => return Status::INVALID_OP.and(Self::NAN), // sqrt(inf) = inf diff --git a/src/lib.rs b/src/lib.rs index 37dd95a..5024294 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -516,7 +516,7 @@ pub trait Float: fn sqrt(self, round: Round) -> StatusAnd { unimplemented!() } - /// + /// /// IEEE-754 _isSignMinus_: Returns true if and only if the current value has /// a negative sign. /// diff --git a/tests/downstream.rs b/tests/downstream.rs index a1ae4fd..eb855d4 100644 --- a/tests/downstream.rs +++ b/tests/downstream.rs @@ -536,7 +536,7 @@ fn fuzz_sqrt() { let hardware = hardware_sqrt::sqrtf32(x, round); // x86 preserves the sign on NaN inputs, others may not (e.g. aarch64 does not). - // + // // non-NaN inputs result in a newly-created NaN if negative, // and newly-created NaN bit patterns from x86 have a negative sign while this library has a positive sign. let eq_value = if cfg!(target_arch = "x86_64") && x.is_nan() { From a70a540ba2c9eea1cd17808e4280d665ee21db2d Mon Sep 17 00:00:00 2001 From: miro Date: Sat, 10 Oct 2026 13:36:40 +0200 Subject: [PATCH 28/31] added f16 sqrt --- fuzz/src/host.rs | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/fuzz/src/host.rs b/fuzz/src/host.rs index cf5a654..8ce2765 100644 --- a/fuzz/src/host.rs +++ b/fuzz/src/host.rs @@ -120,7 +120,6 @@ impl_host_float!(f64, u64); #[cfg(target_has_reliable_f128)] impl_host_float!(f128, u128); -// FIXME: this is missing sqrt impls. /// Assembly implementations on x86 which respect rounding mode. #[cfg(x86_sse2)] #[cfg(target_has_reliable_f16)] @@ -345,7 +344,22 @@ mod x86 { // TODO: write the ASM for other rounding modes. fn sqrt(self, rm: Round) -> Option> { - no_fp_env(rm, || self.sqrt()) + if is_x86_feature_detected!("avx512fp16") { + // WARNING(miro): i'm unable to test this aa my CPU lacks AVX512 + unsafe { + let output: f16; + + let status = with_fp_env!( + rm, + "vsqrtsh {dest}, {src}, {src}", + dest = out(xmm_reg) output, + src = in(xmm_reg) self, + ); + Some(status.and(output)) + } + } else { + no_fp_env(rm, || self.sqrt()) + } } } From c08dbde5b517361588d602d36a802d6150fa7f01 Mon Sep 17 00:00:00 2001 From: miro Date: Sat, 10 Oct 2026 13:47:42 +0200 Subject: [PATCH 29/31] renamed sqrt asm registers --- fuzz/src/host.rs | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/fuzz/src/host.rs b/fuzz/src/host.rs index 8ce2765..3eab1ca 100644 --- a/fuzz/src/host.rs +++ b/fuzz/src/host.rs @@ -506,9 +506,9 @@ mod x86 { let status = with_fp_env!( rm, - "sqrtss {x}, {y}", - x = out(xmm_reg) output, - y = in(xmm_reg) self, + "sqrtss {dest}, {src}", + dest = out(xmm_reg) output, + src = in(xmm_reg) self, ); Some(status.and(output)) } @@ -659,9 +659,9 @@ mod x86 { let status = with_fp_env!( rm, - "sqrtsd {x}, {y}", - x = out(xmm_reg) output, - y = in(xmm_reg) self, + "sqrtsd {dest}, {src}", + dest = out(xmm_reg) output, + src = in(xmm_reg) self, ); Some(status.and(output)) } From 40a41e9058805e46f876585eab180b9c35449509 Mon Sep 17 00:00:00 2001 From: miro Date: Sat, 10 Oct 2026 14:52:22 +0200 Subject: [PATCH 30/31] removed unused import --- tests/downstream.rs | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/downstream.rs b/tests/downstream.rs index eb855d4..945fd6b 100644 --- a/tests/downstream.rs +++ b/tests/downstream.rs @@ -1,7 +1,5 @@ //! Tests added to `rustc_apfloat`, that were not ported from the C++ code. -use std::panic::panic_any; - use rayon::iter::{IntoParallelIterator, ParallelIterator}; use rustc_apfloat::ieee::{BFloat, Double, Float8E4M3FN, Float8E5M2, Half, Quad, Single, X87DoubleExtended}; use rustc_apfloat::{Float, Round, Status, StatusAnd}; From 28d1339149844c696d4678c917041d3daffb3fc4 Mon Sep 17 00:00:00 2001 From: miro Date: Sat, 10 Oct 2026 16:21:10 +0200 Subject: [PATCH 31/31] ignored fuzz_sqrt by default as its very slow --- src/downstream.rs | 3 +++ tests/downstream.rs | 40 +--------------------------------------- 2 files changed, 4 insertions(+), 39 deletions(-) diff --git a/src/downstream.rs b/src/downstream.rs index a56f2e2..a467cb6 100644 --- a/src/downstream.rs +++ b/src/downstream.rs @@ -4,6 +4,9 @@ use crate::{ }; impl IeeeFloat { + // WARNING: make sure to manually run the fuzz_sqrt test when editing this function, + // it's very slow so it's ignored by default. + // `cargo test fuzz_sqrt -- --include-ignored` /// This is a spec conformant implementation of the IEEE Float sqrt function /// This is put in downstream.rs because this function hasn't been implemented in the upstream C++ version yet. pub(crate) fn ieee_sqrt(self, round: Round) -> StatusAnd { diff --git a/tests/downstream.rs b/tests/downstream.rs index 945fd6b..eb096b4 100644 --- a/tests/downstream.rs +++ b/tests/downstream.rs @@ -473,46 +473,8 @@ fn sqrt() { } #[test] +#[ignore = "this test is very slow and causes a timeout in the CI"] fn fuzz_sqrt() { - // for round in [ - // Round::NearestTiesToEven, - // Round::TowardPositive, - // Round::TowardNegative, - // Round::TowardZero, - // ] { - // println!("checking rounding mode {round:?}"); - - // for xi in 0..=u16::MAX { - // if xi % 1_000 == 0 { - // println!("{xi}/{}", u16::MAX); - // } - // let x = f16::from_bits(xi); - // let a = apfloat_sqrtf16(x, round); - // let b = hardware_sqrt::sqrtf16(x, round); - - // // x86 preserves the sign on NaN inputs, others may not (e.g. aarch64 does not). - // let eq_value = if cfg!(target_arch = "x86_64") { - // a.value.bitwise_eq(b.value) - // } else { - // a.value.bitwise_eq(b.value) || (a.value.is_nan() && b.value.is_nan()) - // }; - - // if a.status != b.status || !eq_value { - // panic!( - // "\ - // incorrect result\n\ - // xi: {xi:#010x}\n\ - // x: {x:?}\n\ - // a: {a:?} {af}\n\ - // b: {b:?} {bf}\ - // ", - // af = f32::from_bits(a.value.to_bits().try_into().unwrap()), - // bf = f32::from_bits(a.value.to_bits().try_into().unwrap()), - // ); - // } - // } - // } - for round in [ Round::NearestTiesToEven, Round::TowardPositive,