From e179d6f7d4045ad9381249afb849fcb48831616c Mon Sep 17 00:00:00 2001 From: Lee Dogeon Date: Wed, 18 Mar 2026 16:27:52 +0900 Subject: [PATCH] Follow C99 Annex G --- src/lib.rs | 791 ++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 788 insertions(+), 3 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 661b67b..a1331cb 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -203,6 +203,151 @@ impl Complex { } } +/// `copysign` for `Float` types (num_traits `Float` does not provide it). +#[cfg(any(feature = "std", feature = "libm"))] +#[inline] +fn copysign(magnitude: T, sign: T) -> T { + if sign.is_sign_negative() { + -magnitude.abs() + } else { + magnitude.abs() + } +} + +/// C99 Annex G.5.1 complex multiplication with NaN recovery. +/// +/// Computes `(a + bi) * (c + di)` using the standard formula, then +/// recovers meaningful infinite results when the naive formula produces +/// `NaN + NaN·i` due to `∞ × 0` intermediate products. +#[cfg(any(feature = "std", feature = "libm"))] +fn fmul_impl(a: T, b: T, c: T, d: T) -> Complex { + let mut x = a * c - b * d; + let mut y = a * d + b * c; + + if x.is_nan() && y.is_nan() { + let mut ac = a; + let mut bc = b; + let mut cc = c; + let mut dc = d; + let mut recalc = false; + + // Save original infinity flags before any modification. + let a_inf = a.is_infinite(); + let b_inf = b.is_infinite(); + let c_inf = c.is_infinite(); + let d_inf = d.is_infinite(); + + if a_inf || b_inf { + // (a + bi) is infinite: box infinite components to ±1, finite to ±0 + ac = copysign(if a_inf { T::one() } else { T::zero() }, a); + bc = copysign(if b_inf { T::one() } else { T::zero() }, b); + // Replace NaN in the other operand with signed zero + if cc.is_nan() { + cc = copysign(T::zero(), c); + } + if dc.is_nan() { + dc = copysign(T::zero(), d); + } + recalc = true; + } + + if c_inf || d_inf { + // (c + di) is infinite: box infinite components to ±1, finite to ±0 + cc = copysign(if c_inf { T::one() } else { T::zero() }, c); + dc = copysign(if d_inf { T::one() } else { T::zero() }, d); + // Replace NaN in the other operand with signed zero + if ac.is_nan() { + ac = copysign(T::zero(), a); + } + if bc.is_nan() { + bc = copysign(T::zero(), b); + } + recalc = true; + } + + // If neither operand is infinite but intermediate products overflowed + if !recalc + && ((a * c).is_infinite() + || (b * d).is_infinite() + || (a * d).is_infinite() + || (b * c).is_infinite()) + { + if ac.is_nan() { + ac = copysign(T::zero(), a); + } + if bc.is_nan() { + bc = copysign(T::zero(), b); + } + if cc.is_nan() { + cc = copysign(T::zero(), c); + } + if dc.is_nan() { + dc = copysign(T::zero(), d); + } + recalc = true; + } + + if recalc { + x = T::infinity() * (ac * cc - bc * dc); + y = T::infinity() * (ac * dc + bc * cc); + } + } + + Complex::new(x, y) +} + +/// C99 Annex G.5.2 complex division with NaN recovery. +/// +/// Computes `(a + bi) / (c + di)` using Smith's method for numerical +/// stability, then recovers meaningful results (infinities or zeros) when +/// the naive formula produces `NaN + NaN·i`. +#[cfg(any(feature = "std", feature = "libm"))] +fn fdiv_impl(a: T, b: T, c: T, d: T) -> Complex { + // C99 Annex G.5.2 Case C: finite / infinite → signed zero. + // Handled up-front because Smith's method may produce non-NaN zeros + // with incorrect signs, preventing the post-hoc NaN recovery from + // triggering. + if (c.is_infinite() || d.is_infinite()) && a.is_finite() && b.is_finite() { + let cc = copysign(if c.is_infinite() { T::one() } else { T::zero() }, c); + let dc = copysign(if d.is_infinite() { T::one() } else { T::zero() }, d); + let x = T::zero() * (a * cc + b * dc); + let y = T::zero() * (b * cc - a * dc); + return Complex::new(x, y); + } + + // Smith's method: avoids overflow/underflow by dividing through + // by the larger component of the denominator. + let (mut x, mut y); + if d.abs() <= c.abs() { + let ratio = d / c; + let denom = c + d * ratio; + x = (a + b * ratio) / denom; + y = (b - a * ratio) / denom; + } else { + let ratio = c / d; + let denom = d + c * ratio; + x = (a * ratio + b) / denom; + y = (b * ratio - a) / denom; + } + + // C99 Annex G.5.2 recovery for remaining NaN cases + if x.is_nan() && y.is_nan() { + if c.is_zero() && d.is_zero() && (!a.is_nan() || !b.is_nan()) { + // Case A: nonzero / zero → directed infinity + x = copysign(T::infinity(), c) * a; + y = copysign(T::infinity(), c) * b; + } else if (a.is_infinite() || b.is_infinite()) && c.is_finite() && d.is_finite() { + // Case B: infinite / finite → infinity + let ac = copysign(if a.is_infinite() { T::one() } else { T::zero() }, a); + let bc = copysign(if b.is_infinite() { T::one() } else { T::zero() }, b); + x = T::infinity() * (ac * c + bc * d); + y = T::infinity() * (bc * c - ac * d); + } + } + + Complex::new(x, y) +} + #[cfg(any(feature = "std", feature = "libm"))] impl Complex { /// Create a new Complex with a given phase: `exp(i * phase)`. @@ -569,6 +714,33 @@ impl Complex { ((one + self).ln() - (one - self).ln()) / two } + /// Multiplies `self` by `other` using floating-point operations. + /// + /// This handles infinities more gracefully than the generic `Mul` + /// implementation, following the C99 Annex G.5.1 recovery for cases + /// where the naive formula `(ac-bd) + i(ad+bc)` produces NaN results + /// due to `∞ × 0` intermediate products. + /// + /// # Examples + /// + /// ``` + /// use num_complex::Complex64; + /// + /// // Generic multiplication propagates NaN from ∞ × 0. + /// let a = Complex64::new(f64::INFINITY, f64::INFINITY); + /// let b = Complex64::new(1.0, 0.0); + /// assert!((a * b).re.is_nan()); + /// + /// // But `fmul` recovers the correct infinite result. + /// let result = a.fmul(b); + /// assert_eq!(result.re, f64::INFINITY); + /// assert_eq!(result.im, f64::INFINITY); + /// ``` + #[inline] + pub fn fmul(self, other: Complex) -> Complex { + fmul_impl(self.re, self.im, other.re, other.im) + } + /// Returns `1/self` using floating-point operations. /// /// This may be more accurate than the generic `self.inv()` in cases @@ -599,8 +771,11 @@ impl Complex { /// Returns `self/other` using floating-point operations. /// - /// This may be more accurate than the generic `Div` implementation in cases - /// where `other.norm_sqr()` would overflow to ∞ or underflow to 0. + /// This uses Smith's method for numerical stability and C99 Annex G.5.2 + /// recovery for infinity/NaN edge cases. It is more accurate than the + /// generic `Div` implementation in cases where `other.norm_sqr()` would + /// overflow to ∞ or underflow to 0, and correctly handles division + /// involving infinities and zeros. /// /// # Examples /// @@ -622,7 +797,7 @@ impl Complex { /// ``` #[inline] pub fn fdiv(self, other: Complex) -> Complex { - self * other.finv() + fdiv_impl(self.re, self.im, other.re, other.im) } } @@ -1641,6 +1816,8 @@ pub(crate) mod test { use super::{ComplexErrorKind, ParseComplexError}; use core::f64; use core::str::FromStr; + extern crate std; + use std::vec::Vec; use std::string::{String, ToString}; @@ -2432,6 +2609,614 @@ pub(crate) mod test { )); } } + + /// Assert that two Complex64 values are equal per C99 Annex G semantics: + /// - NaN matches NaN + /// - Signed zero: +0 and -0 are distinguished + /// - Infinities must match sign + fn assert_c99_eq(actual: Complex64, expected: Complex64, msg: &str) { + fn feq(a: f64, b: f64) -> bool { + if b.is_nan() { + a.is_nan() + } else if b == 0.0 && !b.is_nan() { + a == 0.0 && a.is_sign_positive() == b.is_sign_positive() + } else { + a == b + } + } + assert!( + feq(actual.re, expected.re) && feq(actual.im, expected.im), + "{msg}\n actual: ({:+}, {:+})\n expected: ({:+}, {:+})", + actual.re, + actual.im, + expected.re, + expected.im, + ); + } + + // ====================================================================== + // C99 Annex G.5.1: Complex multiplication (fmul) + // ====================================================================== + + #[test] + fn test_fmul_finite() { + // Normal finite multiplication should still work correctly. + assert!(close(_05_05i.fmul(_05_05i), _0_1i.unscale(2.0))); + assert!(close(_1_1i.fmul(_0_1i), _neg1_1i)); + assert!(close(_0_1i.fmul(_0_1i), -_1_0i)); + for &c in all_consts.iter() { + assert!(close(c.fmul(_1_0i), c)); + assert!(close(_1_0i.fmul(c), c)); + } + } + + #[test] + fn test_fmul_step1_first_operand_infinite() { + // Step 1: (a+bi) is infinite (a or b is ∞), other operand finite. + let c = |r, i| Complex64::new(r, i); + let inf = f64::INFINITY; + let ninf = -f64::INFINITY; + let nan = f64::NAN; + + // +∞ components + assert_c99_eq(c(inf, inf).fmul(c(1.0, 0.0)), c(inf, inf), + "(inf+infi)*(1+0i)"); + assert_c99_eq(c(inf, -inf).fmul(c(1.0, 0.0)), c(inf, -inf), + "(inf-infi)*(1+0i)"); + assert_c99_eq(c(inf, 0.0).fmul(c(1.0, 0.0)), c(inf, nan), + "(inf+0i)*(1+0i)"); + assert_c99_eq(c(0.0, inf).fmul(c(1.0, 0.0)), c(nan, inf), + "(0+infi)*(1+0i)"); + assert_c99_eq(c(inf, inf).fmul(c(0.0, 1.0)), c(-inf, inf), + "(inf+infi)*(0+1i)"); + + // -∞ components + assert_c99_eq(c(ninf, 0.0).fmul(c(1.0, 1.0)), c(ninf, ninf), + "(-inf+0i)*(1+1i)"); + assert_c99_eq(c(ninf, inf).fmul(c(1.0, 1.0)), c(ninf, nan), + "(-inf+infi)*(1+1i)"); + assert_c99_eq(c(ninf, ninf).fmul(c(1.0, 0.0)), c(ninf, ninf), + "(-inf-infi)*(1+0i)"); + + // Finite nonzero imaginary part with ∞ real + assert_c99_eq(c(inf, 2.0).fmul(c(1.0, 1.0)), c(inf, inf), + "(inf+2i)*(1+1i)"); + assert_c99_eq(c(inf, 2.0).fmul(c(0.0, 1.0)), c(nan, inf), + "(inf+2i)*(0+1i)"); + assert_c99_eq(c(2.0, inf).fmul(c(1.0, 1.0)), c(-inf, inf), + "(2+infi)*(1+1i)"); + + // NaN component with ∞ + assert_c99_eq(c(inf, nan).fmul(c(1.0, 0.0)), c(inf, nan), + "(inf+nani)*(1+0i)"); + assert_c99_eq(c(inf, nan).fmul(c(1.0, 1.0)), c(inf, inf), + "(inf+nani)*(1+1i)"); + assert_c99_eq(c(nan, inf).fmul(c(1.0, 0.0)), c(nan, inf), + "(nan+infi)*(1+0i)"); + assert_c99_eq(c(ninf, nan).fmul(c(1.0, 1.0)), c(ninf, ninf), + "(-inf+nani)*(1+1i)"); + assert_c99_eq(c(nan, ninf).fmul(c(1.0, 1.0)), c(inf, ninf), + "(nan-infi)*(1+1i)"); + } + + #[test] + fn test_fmul_step2_second_operand_infinite() { + // Step 2: (c+di) is infinite, first operand finite. + let c = |r, i| Complex64::new(r, i); + let inf = f64::INFINITY; + let ninf = -f64::INFINITY; + let nan = f64::NAN; + + assert_c99_eq(c(1.0, 0.0).fmul(c(inf, inf)), c(inf, inf), + "(1+0i)*(inf+infi)"); + assert_c99_eq(c(1.0, 0.0).fmul(c(nan, inf)), c(nan, inf), + "(1+0i)*(nan+infi)"); + assert_c99_eq(c(1.0, 1.0).fmul(c(ninf, 0.0)), c(ninf, ninf), + "(1+1i)*(-inf+0i)"); + assert_c99_eq(c(1.0, 1.0).fmul(c(ninf, inf)), c(ninf, nan), + "(1+1i)*(-inf+infi)"); + assert_c99_eq(c(1.0, 0.0).fmul(c(0.0, ninf)), c(nan, ninf), + "(1+0i)*(0-infi)"); + assert_c99_eq(c(1.0, 1.0).fmul(c(0.0, ninf)), c(inf, ninf), + "(1+1i)*(0-infi)"); + assert_c99_eq(c(2.0, 3.0).fmul(c(ninf, inf)), c(ninf, nan), + "(2+3i)*(-inf+infi)"); + } + + #[test] + fn test_fmul_both_infinite() { + // Steps 1+2: both operands infinite. + let c = |r, i| Complex64::new(r, i); + let inf = f64::INFINITY; + let ninf = -f64::INFINITY; + let nan = f64::NAN; + + assert_c99_eq(c(inf, 0.0).fmul(c(inf, 0.0)), c(inf, nan), + "(inf+0i)*(inf+0i)"); + assert_c99_eq(c(inf, 0.0).fmul(c(-inf, 0.0)), c(-inf, nan), + "(inf+0i)*(-inf+0i)"); + assert_c99_eq(c(inf, inf).fmul(c(inf, inf)), c(nan, inf), + "(inf+infi)*(inf+infi)"); + assert_c99_eq(c(inf, nan).fmul(c(inf, nan)), c(inf, nan), + "(inf+nani)*(inf+nani)"); + assert_c99_eq(c(inf, inf).fmul(c(ninf, inf)), c(ninf, nan), + "(inf+infi)*(-inf+infi)"); + assert_c99_eq(c(ninf, inf).fmul(c(ninf, ninf)), c(inf, nan), + "(-inf+infi)*(-inf-infi)"); + assert_c99_eq(c(inf, ninf).fmul(c(ninf, inf)), c(nan, inf), + "(inf-infi)*(-inf+infi)"); + assert_c99_eq(c(ninf, nan).fmul(c(inf, nan)), c(ninf, nan), + "(-inf+nani)*(inf+nani)"); + assert_c99_eq(c(nan, inf).fmul(c(nan, ninf)), c(inf, nan), + "(nan+infi)*(nan-infi)"); + assert_c99_eq(c(inf, 0.0).fmul(c(0.0, inf)), c(nan, inf), + "(inf+0i)*(0+infi)"); + assert_c99_eq(c(inf, 0.0).fmul(c(0.0, ninf)), c(nan, ninf), + "(inf+0i)*(0-infi)"); + assert_c99_eq(c(ninf, 0.0).fmul(c(0.0, inf)), c(nan, ninf), + "(-inf+0i)*(0+infi)"); + } + + #[test] + fn test_fmul_step3_intermediate_overflow() { + // Step 3: finite operands but intermediate products overflow to ∞. + let c = |r, i| Complex64::new(r, i); + let inf = f64::INFINITY; + let big = f64::MAX; + + // big*big overflows to ∞, so at least one intermediate is ∞. + // (big+big*i)*(big+0i): x=big²-0=∞, y=0+big²=∞ → no recovery needed (not NaN) + assert_c99_eq(c(big, big).fmul(c(big, 0.0)), c(inf, inf), + "(MAX+MAXi)*(MAX+0i)"); + // (big+big*i)*(big+big*i): x=big²-big²=∞-∞=NaN, y=big²+big²=∞ → only re NaN + assert_c99_eq(c(big, big).fmul(c(big, big)), c(f64::NAN, inf), + "(MAX+MAXi)*(MAX+MAXi)"); + } + + #[test] + fn test_fmul_inf_times_zero() { + // ∞ × 0: indeterminate, recovery produces NaN. + let c = |r, i| Complex64::new(r, i); + let inf = f64::INFINITY; + let ninf = -f64::INFINITY; + let nan = f64::NAN; + + assert_c99_eq(c(inf, 0.0).fmul(c(0.0, 0.0)), c(nan, nan), + "(inf+0i)*(0+0i)"); + assert_c99_eq(c(ninf, 0.0).fmul(c(0.0, 0.0)), c(nan, nan), + "(-inf+0i)*(0+0i)"); + assert_c99_eq(c(0.0, inf).fmul(c(0.0, 0.0)), c(nan, nan), + "(0+infi)*(0+0i)"); + assert_c99_eq(c(inf, inf).fmul(c(0.0, 0.0)), c(nan, nan), + "(inf+infi)*(0+0i)"); + } + + #[test] + fn test_fmul_nan_no_recovery() { + // Pure NaN operands: no recovery possible. + let c = |r, i| Complex64::new(r, i); + let nan = f64::NAN; + + assert_c99_eq(c(nan, nan).fmul(c(1.0, 0.0)), c(nan, nan), + "(nan+nani)*(1+0i)"); + assert_c99_eq(c(nan, nan).fmul(c(nan, nan)), c(nan, nan), + "(nan+nani)*(nan+nani)"); + assert_c99_eq(c(nan, 0.0).fmul(c(0.0, 0.0)), c(nan, nan), + "(nan+0i)*(0+0i)"); + assert_c99_eq(c(nan, 0.0).fmul(c(1.0, 0.0)), c(nan, nan), + "(nan+0i)*(1+0i)"); + assert_c99_eq(c(0.0, nan).fmul(c(0.0, 1.0)), c(nan, nan), + "(0+nani)*(0+1i)"); + } + + // ====================================================================== + // C99 Annex G.5.2: Complex division (fdiv) + // ====================================================================== + + #[test] + fn test_fdiv_finite() { + // Normal finite division. + assert!(close(_neg1_1i.fdiv(_0_1i), _1_1i)); + for &c in all_consts.iter() { + if !c.is_zero() { + assert!(close(c.fdiv(c), _1_0i)); + } + } + } + + #[test] + fn test_fdiv_case_a_nonzero_over_zero() { + // Case A: c=0, d=0, numerator not all-NaN → directed ∞. + let c = |r, i| Complex64::new(r, i); + let inf = f64::INFINITY; + let ninf = -f64::INFINITY; + let nan = f64::NAN; + let nz: f64 = -0.0; + + assert_c99_eq(c(1.0, 0.0).fdiv(c(0.0, 0.0)), c(inf, nan), + "(1+0i)/(0+0i)"); + assert_c99_eq(c(1.0, 1.0).fdiv(c(0.0, 0.0)), c(inf, inf), + "(1+1i)/(0+0i)"); + assert_c99_eq(c(-1.0, 0.0).fdiv(c(0.0, 0.0)), c(ninf, nan), + "(-1+0i)/(0+0i)"); + + // Negative zero in denominator changes sign of ∞ + assert_c99_eq(c(1.0, 0.0).fdiv(c(nz, 0.0)), c(ninf, nan), + "(1+0i)/(-0+0i)"); + // -0 in imaginary part of denominator doesn't affect (copysign uses c) + assert_c99_eq(c(1.0, 0.0).fdiv(c(0.0, nz)), c(inf, nan), + "(1+0i)/(0-0i)"); + + // Imaginary-only numerator + assert_c99_eq(c(0.0, 1.0).fdiv(c(0.0, 0.0)), c(nan, inf), + "(0+1i)/(0+0i)"); + // Negative numerator + assert_c99_eq(c(-1.0, -1.0).fdiv(c(0.0, 0.0)), c(ninf, ninf), + "(-1-1i)/(0+0i)"); + + // 0/0 → NaN (numerator is zero, ∞*0=NaN) + assert_c99_eq(c(0.0, 0.0).fdiv(c(0.0, 0.0)), c(nan, nan), + "(0+0i)/(0+0i)"); + // Partial NaN in numerator: !nan.is_nan() || !0.is_nan() → true + assert_c99_eq(c(nan, 0.0).fdiv(c(0.0, 0.0)), c(nan, nan), + "(nan+0i)/(0+0i)"); + // All-NaN numerator: !nan.is_nan() || !nan.is_nan() → false → no recovery + assert_c99_eq(c(nan, nan).fdiv(c(0.0, 0.0)), c(nan, nan), + "(nan+nani)/(0+0i)"); + + // ∞/0: ∞ is not NaN, so recovery gives ∞*∞=∞ + assert_c99_eq(c(inf, 0.0).fdiv(c(0.0, 0.0)), c(inf, nan), + "(inf+0i)/(0+0i)"); + } + + #[test] + fn test_fdiv_case_b_infinite_over_finite() { + // Case B: numerator infinite, denominator finite → ∞. + let c = |r, i| Complex64::new(r, i); + let inf = f64::INFINITY; + let ninf = -f64::INFINITY; + let nan = f64::NAN; + + assert_c99_eq(c(inf, 0.0).fdiv(c(1.0, 0.0)), c(inf, nan), + "(inf+0i)/(1+0i)"); + assert_c99_eq(c(ninf, 0.0).fdiv(c(1.0, 0.0)), c(ninf, nan), + "(-inf+0i)/(1+0i)"); + assert_c99_eq(c(0.0, inf).fdiv(c(1.0, 0.0)), c(nan, inf), + "(0+infi)/(1+0i)"); + assert_c99_eq(c(0.0, ninf).fdiv(c(1.0, 0.0)), c(nan, ninf), + "(0-infi)/(1+0i)"); + assert_c99_eq(c(inf, inf).fdiv(c(1.0, 0.0)), c(inf, inf), + "(inf+infi)/(1+0i)"); + assert_c99_eq(c(ninf, inf).fdiv(c(1.0, 0.0)), c(ninf, inf), + "(-inf+infi)/(1+0i)"); + assert_c99_eq(c(inf, inf).fdiv(c(1.0, 1.0)), c(inf, nan), + "(inf+infi)/(1+1i)"); + assert_c99_eq(c(ninf, ninf).fdiv(c(1.0, 1.0)), c(ninf, nan), + "(-inf-infi)/(1+1i)"); + assert_c99_eq(c(inf, inf).fdiv(c(-1.0, 1.0)), c(nan, ninf), + "(inf+infi)/(-1+1i)"); + assert_c99_eq(c(inf, nan).fdiv(c(1.0, 0.0)), c(inf, nan), + "(inf+nani)/(1+0i)"); + assert_c99_eq(c(inf, nan).fdiv(c(1.0, 1.0)), c(inf, ninf), + "(inf+nani)/(1+1i)"); + assert_c99_eq(c(nan, inf).fdiv(c(1.0, 0.0)), c(nan, inf), + "(nan+infi)/(1+0i)"); + assert_c99_eq(c(nan, ninf).fdiv(c(0.0, 1.0)), c(ninf, nan), + "(nan-infi)/(0+1i)"); + } + + #[test] + fn test_fdiv_case_c_finite_over_infinite() { + // Case C: numerator finite, denominator infinite → signed zero. + let c = |r, i| Complex64::new(r, i); + let inf = f64::INFINITY; + let ninf = -f64::INFINITY; + + assert_c99_eq(c(1.0, 0.0).fdiv(c(inf, 0.0)), c(0.0, 0.0), + "(1+0i)/(inf+0i)"); + assert_c99_eq(c(1.0, 0.0).fdiv(c(ninf, 0.0)), c(-0.0, -0.0), + "(1+0i)/(-inf+0i)"); + assert_c99_eq(c(1.0, 1.0).fdiv(c(inf, inf)), c(0.0, 0.0), + "(1+1i)/(inf+infi)"); + assert_c99_eq(c(1.0, 0.0).fdiv(c(inf, inf)), c(0.0, -0.0), + "(1+0i)/(inf+infi)"); + assert_c99_eq(c(0.0, 1.0).fdiv(c(inf, inf)), c(0.0, 0.0), + "(0+1i)/(inf+infi)"); + assert_c99_eq(c(1.0, 0.0).fdiv(c(0.0, inf)), c(0.0, -0.0), + "(1+0i)/(0+infi)"); + // Negative infinity in denominator + assert_c99_eq(c(1.0, 0.0).fdiv(c(0.0, ninf)), c(0.0, 0.0), + "(1+0i)/(0-infi)"); + assert_c99_eq(c(-1.0, -1.0).fdiv(c(inf, inf)), c(-0.0, 0.0), + "(-1-1i)/(inf+infi)"); + assert_c99_eq(c(1.0, 0.0).fdiv(c(inf, ninf)), c(0.0, 0.0), + "(1+0i)/(inf-infi)"); + assert_c99_eq(c(1.0, 1.0).fdiv(c(ninf, inf)), c(0.0, -0.0), + "(1+1i)/(-inf+infi)"); + // Zero numerator / infinite denominator → +0 + assert_c99_eq(c(0.0, 0.0).fdiv(c(inf, 0.0)), c(0.0, 0.0), + "(0+0i)/(inf+0i)"); + } + + #[test] + fn test_fdiv_infinite_over_infinite() { + // ∞/∞ is indeterminate → NaN+NaN·i. + let c = |r, i| Complex64::new(r, i); + let inf = f64::INFINITY; + let ninf = -f64::INFINITY; + let nan = f64::NAN; + + assert_c99_eq(c(inf, 0.0).fdiv(c(inf, 0.0)), c(nan, nan), + "(inf+0i)/(inf+0i)"); + assert_c99_eq(c(inf, inf).fdiv(c(inf, inf)), c(nan, nan), + "(inf+infi)/(inf+infi)"); + assert_c99_eq(c(inf, inf).fdiv(c(ninf, inf)), c(nan, nan), + "(inf+infi)/(-inf+infi)"); + } + + #[test] + fn test_fdiv_nan_no_recovery() { + // Pure NaN: no recovery possible. + let c = |r, i| Complex64::new(r, i); + let nan = f64::NAN; + + assert_c99_eq(c(nan, nan).fdiv(c(1.0, 0.0)), c(nan, nan), + "(nan+nani)/(1+0i)"); + assert_c99_eq(c(1.0, 0.0).fdiv(c(nan, nan)), c(nan, nan), + "(1+0i)/(nan+nani)"); + assert_c99_eq(c(nan, 0.0).fdiv(c(nan, 0.0)), c(nan, nan), + "(nan+0i)/(nan+0i)"); + } + + #[test] + fn test_fdiv_overflow_stability() { + // Verify fdiv handles large values better than generic div. + let a = Complex64::new(2.0, 3.0); + let b = Complex64::new(1e300, 1e300); + + // Generic division overflows. + assert!(!(a / b).is_normal()); + + // fdiv with Smith's method handles it. + let quotient = a.fdiv(b); + assert!(quotient.is_normal()); + let expected = Complex64::new(2.5e-300, 5e-301); + assert!((quotient - expected).norm() < 1e-315); + } + + // ====================================================================== + // LLVM compiler-rt style classification test + // + // Ported from: + // llvm-project/compiler-rt/test/builtins/Unit/muldc3_test.c + // llvm-project/compiler-rt/test/builtins/Unit/divdc3_test.c + // + // Exhaustively tests 121×121 = 14,641 pairs for both fmul and fdiv, + // validating that the result category (zero/nonzero/inf/NaN/nonzero_nan) + // matches the expected 5×5 matrix derived from C99 Annex G. + // ====================================================================== + + #[derive(Debug, Clone, Copy, PartialEq, Eq)] + enum Category { + Zero, + NonZero, + Inf, + NaN, + NonZeroNan, + } + + fn classify(z: Complex64) -> Category { + let (re, im) = (z.re, z.im); + if re == 0.0 && im == 0.0 { + return Category::Zero; + } + if re.is_infinite() || im.is_infinite() { + return Category::Inf; + } + if re.is_nan() && im.is_nan() { + return Category::NaN; + } + if re.is_nan() { + if im == 0.0 { + return Category::NaN; + } + return Category::NonZeroNan; + } + if im.is_nan() { + if re == 0.0 { + return Category::NaN; + } + return Category::NonZeroNan; + } + Category::NonZero + } + + /// 11 representative component values: NaN, -∞, finite, ±0, +∞ + fn test_components() -> [f64; 11] { + use core::f64::{INFINITY, NAN}; + [ + NAN, + -INFINITY, + -2.0, + -1.0, + -0.5, + -0.0, + 0.0, + 0.5, + 1.0, + 2.0, + INFINITY, + ] + } + + /// All 121 complex numbers formed from the 11 components. + fn test_values() -> Vec { + let comp = test_components(); + let mut v: Vec = Vec::with_capacity(121); + for &re in &comp { + for &im in &comp { + v.push(Complex64::new(re, im)); + } + } + v + } + + /// C99 Annex G expected result category for multiplication. + /// + /// Rows = classify(lhs), Columns = classify(rhs). + fn expected_fmul_category(a: Category, b: Category) -> Category { + use Category::*; + match (a, b) { + (Zero, Zero) => Zero, + (Zero, NonZero) | (NonZero, Zero) => Zero, + (Zero, Inf) | (Inf, Zero) => NaN, + (Zero, NaN) | (NaN, Zero) => NaN, + (Zero, NonZeroNan) | (NonZeroNan, Zero) => NaN, + + (NonZero, NonZero) => NonZero, + (NonZero, Inf) | (Inf, NonZero) => Inf, + (NonZero, NaN) | (NaN, NonZero) => NaN, + (NonZero, NonZeroNan) | (NonZeroNan, NonZero) => NaN, + + (Inf, Inf) => Inf, + (Inf, NaN) | (NaN, Inf) => NaN, + (Inf, NonZeroNan) | (NonZeroNan, Inf) => Inf, + + (NaN, NaN) => NaN, + (NaN, NonZeroNan) | (NonZeroNan, NaN) => NaN, + + (NonZeroNan, NonZeroNan) => NaN, + } + } + + /// C99 Annex G expected result category for division. + /// + /// Rows = classify(numerator), Columns = classify(denominator). + fn expected_fdiv_category(a: Category, b: Category) -> Category { + use Category::*; + match (a, b) { + (Zero, Zero) => NaN, + (Zero, NonZero) => Zero, + (Zero, Inf) => Zero, + (Zero, NaN) => NaN, + (Zero, NonZeroNan) => NaN, + + (NonZero, Zero) => Inf, + (NonZero, NonZero) => NonZero, + (NonZero, Inf) => Zero, + (NonZero, NaN) => NaN, + (NonZero, NonZeroNan) => NaN, + + (Inf, Zero) => Inf, + (Inf, NonZero) => Inf, + (Inf, Inf) => NaN, + (Inf, NaN) => NaN, + (Inf, NonZeroNan) => NaN, + + (NaN, Zero) => NaN, + (NaN, NonZero) => NaN, + (NaN, Inf) => NaN, + (NaN, NaN) => NaN, + (NaN, NonZeroNan) => NaN, + + (NonZeroNan, Zero) => Inf, + (NonZeroNan, NonZero) => NaN, + (NonZeroNan, Inf) => NaN, + (NonZeroNan, NaN) => NaN, + (NonZeroNan, NonZeroNan) => NaN, + } + } + + #[test] + fn test_fmul_classification_exhaustive() { + // LLVM compiler-rt muldc3_test: 121×121 = 14,641 pairs + let vals = test_values(); + let mut failures: Vec = Vec::new(); + for &a in &vals { + for &b in &vals { + let result: Complex64 = a.fmul(b); + let cat_a = classify(a); + let cat_b = classify(b); + let expected = expected_fmul_category(cat_a, cat_b); + let actual = classify(result); + if actual != expected { + failures.push(format!( + "fmul({:?}, {:?}): {:?}*{:?}={:?}, expected {:?}, got {:?}", + a, b, cat_a, cat_b, result, expected, actual + )); + } + // For NonZero * NonZero, also check approximate correctness + if cat_a == Category::NonZero && cat_b == Category::NonZero { + let naive = Complex64::new( + a.re * b.re - a.im * b.im, + a.re * b.im + a.im * b.re, + ); + let rel = (result - naive).norm() / naive.norm(); + if rel > 1e-9 { + failures.push(format!( + "fmul({:?}, {:?}): accuracy {:.2e} > 1e-9", + a, b, rel + )); + } + } + } + } + if !failures.is_empty() { + let count = failures.len(); + let mut msg = std::format!("{count} classification failures in fmul (showing first 20):\n"); + for f in failures.iter().take(20) { + msg.push_str(f); + msg.push('\n'); + } + panic!("{}", msg); + } + } + + #[test] + fn test_fdiv_classification_exhaustive() { + // LLVM compiler-rt divdc3_test: 121×121 = 14,641 pairs + let vals = test_values(); + let mut failures: Vec = Vec::new(); + for &a in &vals { + for &b in &vals { + let result: Complex64 = a.fdiv(b); + let cat_a = classify(a); + let cat_b = classify(b); + let expected = expected_fdiv_category(cat_a, cat_b); + let actual = classify(result); + if actual != expected { + failures.push(format!( + "fdiv({:?}, {:?}): {:?}/{:?}={:?}, expected {:?}, got {:?}", + a, b, cat_a, cat_b, result, expected, actual + )); + } + // For NonZero / NonZero, also check approximate correctness + if cat_a == Category::NonZero && cat_b == Category::NonZero { + let denom = b.re * b.re + b.im * b.im; + let naive = Complex64::new( + (a.re * b.re + a.im * b.im) / denom, + (a.im * b.re - a.re * b.im) / denom, + ); + let rel = (result - naive).norm() / result.norm(); + if rel > 1e-6 { + failures.push(format!( + "fdiv({:?}, {:?}): accuracy {:.2e} > 1e-6", + a, b, rel + )); + } + } + } + } + if !failures.is_empty() { + let count = failures.len(); + let mut msg = std::format!("{count} classification failures in fdiv (showing first 20):\n"); + for f in failures.iter().take(20) { + msg.push_str(f); + msg.push('\n'); + } + panic!("{}", msg); + } + } } // Test both a + b and a += b