From 898ffc5f74fad1593ed4036a69da039937224f99 Mon Sep 17 00:00:00 2001 From: Linus Hamlin <78953007+lilinus@users.noreply.github.com> Date: Wed, 3 Jul 2024 22:41:35 +0200 Subject: [PATCH] Complex.Sqrt comply with IEEE 754 and use double.Hypot (#104262) * Use double.Hypot in Complex * Replace Complex.Hypot with double.Hypot * Make Complex.Sqrt IEEE 754 compliant * Fix Complex.Abs and Complex.Sqrt for IEEE 754 * Update Complex.cs --- .../src/System/Numerics/Complex.cs | 76 +++++++------------ .../tests/ComplexTests.cs | 25 ++++-- 2 files changed, 48 insertions(+), 53 deletions(-) diff --git a/src/libraries/System.Runtime.Numerics/src/System/Numerics/Complex.cs b/src/libraries/System.Runtime.Numerics/src/System/Numerics/Complex.cs index 9c4f2d640aa11..2eb4458ce38fc 100644 --- a/src/libraries/System.Runtime.Numerics/src/System/Numerics/Complex.cs +++ b/src/libraries/System.Runtime.Numerics/src/System/Numerics/Complex.cs @@ -289,46 +289,7 @@ public static Complex Divide(double dividend, Complex divisor) public static double Abs(Complex value) { - return Hypot(value.m_real, value.m_imaginary); - } - - private static double Hypot(double a, double b) - { - // Using - // sqrt(a^2 + b^2) = |a| * sqrt(1 + (b/a)^2) - // we can factor out the larger component to dodge overflow even when a * a would overflow. - - a = Math.Abs(a); - b = Math.Abs(b); - - double small, large; - if (a < b) - { - small = a; - large = b; - } - else - { - small = b; - large = a; - } - - if (small == 0.0) - { - return (large); - } - else if (double.IsPositiveInfinity(large) && !double.IsNaN(small)) - { - // The NaN test is necessary so we don't return +inf when small=NaN and large=+inf. - // NaN in any other place returns NaN without any special handling. - return (double.PositiveInfinity); - } - else - { - double ratio = small / large; - return (large * Math.Sqrt(1.0 + ratio * ratio)); - } - + return double.Hypot(value.m_real, value.m_imaginary); } private static double Log1P(double x) @@ -590,8 +551,8 @@ private static void Asin_Internal(double x, double y, out double b, out double b } else { - double r = Hypot((x + 1.0), y); - double s = Hypot((x - 1.0), y); + double r = double.Hypot((x + 1.0), y); + double s = double.Hypot((x - 1.0), y); double a = (r + s) * 0.5; b = x / a; @@ -673,6 +634,28 @@ public static Complex Exp(Complex value) public static Complex Sqrt(Complex value) { + // Handle NaN input cases according to IEEE 754 + if (double.IsNaN(value.m_real)) + { + if (double.IsInfinity(value.m_imaginary)) + { + return new Complex(double.PositiveInfinity, value.m_imaginary); + } + return new Complex(double.NaN, double.NaN); + } + if (double.IsNaN(value.m_imaginary)) + { + if (double.IsPositiveInfinity(value.m_real)) + { + return new Complex(double.NaN, double.PositiveInfinity); + } + if (double.IsNegativeInfinity(value.m_real)) + { + return new Complex(double.PositiveInfinity, double.NaN); + } + return new Complex(double.NaN, double.NaN); + } + if (value.m_imaginary == 0.0) { // Handle the trivial case quickly. @@ -715,11 +698,10 @@ public static Complex Sqrt(Complex value) double imaginaryCopy = value.m_imaginary; if ((Math.Abs(realCopy) >= s_sqrtRescaleThreshold) || (Math.Abs(imaginaryCopy) >= s_sqrtRescaleThreshold)) { - if (double.IsInfinity(value.m_imaginary) && !double.IsNaN(value.m_real)) + if (double.IsInfinity(value.m_imaginary)) { // We need to handle infinite imaginary parts specially because otherwise - // our formulas below produce inf/inf = NaN. The NaN test is necessary - // so that we return NaN rather than (+inf,inf) for (NaN,inf). + // our formulas below produce inf/inf = NaN. return (new Complex(double.PositiveInfinity, imaginaryCopy)); } @@ -732,12 +714,12 @@ public static Complex Sqrt(Complex value) double x, y; if (realCopy >= 0.0) { - x = Math.Sqrt((Hypot(realCopy, imaginaryCopy) + realCopy) * 0.5); + x = Math.Sqrt((double.Hypot(realCopy, imaginaryCopy) + realCopy) * 0.5); y = imaginaryCopy / (2.0 * x); } else { - y = Math.Sqrt((Hypot(realCopy, imaginaryCopy) - realCopy) * 0.5); + y = Math.Sqrt((double.Hypot(realCopy, imaginaryCopy) - realCopy) * 0.5); if (imaginaryCopy < 0.0) y = -y; x = imaginaryCopy / (2.0 * y); } diff --git a/src/libraries/System.Runtime.Numerics/tests/ComplexTests.cs b/src/libraries/System.Runtime.Numerics/tests/ComplexTests.cs index bb5497df33962..fc843295b2e51 100644 --- a/src/libraries/System.Runtime.Numerics/tests/ComplexTests.cs +++ b/src/libraries/System.Runtime.Numerics/tests/ComplexTests.cs @@ -225,7 +225,7 @@ public static IEnumerable Abs_Advanced_TestData() foreach (double invalidImaginary in s_invalidDoubleValues) { yield return new object[] { RandomPositiveDouble(), invalidImaginary, Math.Abs(invalidImaginary) }; // Invalid imaginary - yield return new object[] { invalidReal, invalidImaginary, (double.IsNaN(invalidReal) || double.IsNaN(invalidImaginary)) ? double.NaN : double.PositiveInfinity }; // Invalid real, invalid imaginary + yield return new object[] { invalidReal, invalidImaginary, (double.IsNaN(invalidReal) && double.IsNaN(invalidImaginary)) ? double.NaN : double.PositiveInfinity }; // Invalid real, invalid imaginary } } @@ -245,11 +245,16 @@ public static IEnumerable Abs_Advanced_TestData() yield return new object[] { double.MaxValue, double.NegativeInfinity, double.PositiveInfinity }; yield return new object[] { double.PositiveInfinity, double.NegativeInfinity, double.PositiveInfinity }; - // NaN in any slot returns NaN. + // NaN and +-inf returns inf + yield return new object[] { double.NaN, double.NegativeInfinity, double.PositiveInfinity }; + yield return new object[] { double.NaN, double.PositiveInfinity, double.PositiveInfinity }; + yield return new object[] { double.PositiveInfinity, double.NaN, double.PositiveInfinity }; + yield return new object[] { double.NegativeInfinity, double.NaN, double.PositiveInfinity }; + + // Otherwise, NaN in any slot returns NaN. yield return new object[] { double.NaN, 0, double.NaN }; // Regression test: Complex.Abs() is inconsistent on NaN / Complex yield return new object[] { 0.0, double.NaN, double.NaN }; yield return new object[] { double.MaxValue, double.NaN, double.NaN }; - yield return new object[] { double.NaN, double.NegativeInfinity, double.NaN }; yield return new object[] { double.NaN, double.NaN, double.NaN }; } @@ -1560,12 +1565,20 @@ public static IEnumerable Sqrt_AdvancedTestData () yield return new object[] { RandomPositiveDouble(), double.NegativeInfinity, double.PositiveInfinity, double.NegativeInfinity }; yield return new object[] { RandomNegativeDouble(), double.NegativeInfinity, double.PositiveInfinity, double.NegativeInfinity }; - // NaN in any component produces NaNs in both components (except arguably on real line). + // (Nan, +-inf) returns (inf, +-inf) + yield return new object[] { double.NaN, double.PositiveInfinity, double.PositiveInfinity, double.PositiveInfinity }; + yield return new object[] { double.NaN, double.NegativeInfinity, double.PositiveInfinity, double.NegativeInfinity }; + + // (inf, NaN) returns (inf, NaN) + yield return new object[] { double.PositiveInfinity, double.NaN, double.NaN, double.PositiveInfinity }; + + // (-inf, NaN) returns (NaN, inf) + yield return new object[] { double.NegativeInfinity, double.NaN, double.PositiveInfinity, double.NaN }; + + // Otherwise, NaN in any component produces NaNs in both components. yield return new object[] { 0.0, double.NaN, double.NaN, double.NaN }; yield return new object[] { double.MaxValue, double.NaN, double.NaN, double.NaN }; - yield return new object[] { double.NegativeInfinity, double.NaN, double.NaN, double.NaN }; yield return new object[] { double.NaN, -double.MaxValue, double.NaN, double.NaN }; - yield return new object[] { double.NaN, double.PositiveInfinity, double.NaN, double.NaN }; yield return new object[] { double.NaN, double.NaN, double.NaN, double.NaN }; }