Skip to content

Commit

Permalink
Complex.Sqrt comply with IEEE 754 and use double.Hypot (dotnet#104262)
Browse files Browse the repository at this point in the history 
* 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
  • Loading branch information
@lilinus
lilinus committed Jul 3, 2024
1 parent 89d63f9 commit 898ffc5
Show file tree
Hide file tree
Showing 2 changed files with 48 additions and 53 deletions.
76 changes: 29 additions & 47 deletions src/libraries/System.Runtime.Numerics/src/System/Numerics/Complex.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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));
}

Expand All @@ -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);
}
Expand Down
25 changes: 19 additions & 6 deletions src/libraries/System.Runtime.Numerics/tests/ComplexTests.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -225,7 +225,7 @@ public static IEnumerable<object[]> 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
}
}

Expand All @@ -245,11 +245,16 @@ public static IEnumerable<object[]> 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 };
}

Expand Down Expand Up @@ -1560,12 +1565,20 @@ public static IEnumerable<object[]> 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 };

}
Expand Down

0 comments on commit 898ffc5

Please sign in to comment.