This is an automated email from the ASF dual-hosted git repository. aherbert pushed a commit to branch master in repository https://gitbox.apache.org/repos/asf/commons-numbers.git
commit 1d63ae2bd28ccce68832f6348ad37961504908a6 Author: aherbert <[email protected]> AuthorDate: Mon Dec 9 15:46:34 2019 +0000 Ensure sqrt() always functions on finite input. Any finite complex can be square rooted. If computation of the new magnitude overflows then the computation can be scaled down, computed then scaled back. This ensures infinity is never returned for a finite complex. --- .../apache/commons/numbers/complex/Complex.java | 49 +++++++++++++++++----- .../commons/numbers/complex/CReferenceTest.java | 17 +++++++- 2 files changed, 54 insertions(+), 12 deletions(-) diff --git a/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java b/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java index 8ce7574..8e57ac0 100644 --- a/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java +++ b/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java @@ -1704,19 +1704,33 @@ public final class Complex implements Serializable { } return new Complex(sqrtAbs, imaginary); } - final double abs = getAbsolute(real, imaginary); - final double av = (Math.abs(real) + abs) / 2; - if (av == Double.POSITIVE_INFINITY) { - // Compute in polar coords. - // This handles extreme values that fail in the cartesian representation. - final double sqrtAbs = Math.sqrt(abs); - final double halfArg = getArgument(real, imaginary) / 2; - final double re = sqrtAbs * Math.cos(halfArg); - final double im = sqrtAbs * Math.sin(halfArg); - return new Complex(re, im); + // Get the absolute of the real + double absA = Math.abs(real); + // Compute |a + b i| + double absC = getAbsolute(real, imaginary); + + // t = sqrt((|a| + |a + b i|) / 2) + // This is always representable as this complex is finite. + double t; + + // Overflow safe + if (absC == Double.POSITIVE_INFINITY) { + // Complex is too large. + // Divide by the largest absolute component, + // compute the required sqrt and then scale back. + // Use the equality: sqrt(n) = sqrt(s) * sqrt(n/s) + // t = sqrt(max) * sqrt((|a|/max + |a + b i|/max) / 2) + // Note: The function may be non-monotonic at the junction. + // The alternative of returning infinity for a finite input is worse. + final double max = Math.max(absA, Math.abs(imaginary)); + absA /= max; + absC = getAbsolute(absA, imaginary / max); + t = Math.sqrt(max) * Math.sqrt((absA + absC) / 2); + } else { + // Over-flow safe average + t = Math.sqrt(average(absA, absC)); } - final double t = Math.sqrt(av); if (real >= 0) { return new Complex(t, imaginary / (2 * t)); } @@ -1740,6 +1754,19 @@ public final class Complex implements Serializable { } /** + * Compute the average of two positive finite values in an overflow safe manner. + * + * @param a the first value + * @param b the second value + * @return the average + */ + private static double average(double a, double b) { + return (a < b) ? + a + (b - a) / 2 : + b + (a - b) / 2; + } + + /** * Compute the * <a href="http://mathworld.wolfram.com/Tangent.html";> * tangent</a> of this complex number. diff --git a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java index 2e817b7..662a972 100644 --- a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java +++ b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java @@ -35,7 +35,6 @@ import java.util.function.UnaryOperator; */ public class CReferenceTest { private static final double inf = Double.POSITIVE_INFINITY; - private static final double nan = Double.NaN; /** * Assert the two numbers are equal within the provided units of least precision. @@ -401,6 +400,7 @@ public class CReferenceTest { public void testSqrt() { // Note: When computed in polar coordinates: // real = (x^2 + y^2)^0.25 * cos(0.5 * atan2(y, x)) + // imag = (x^2 + y^2)^0.25 * sin(0.5 * atan2(y, x)) // If x is positive and y is +/-0.0 atan2 returns +/-0. // If x is negative and y is +/-0.0 atan2 returns +/-PI. // This causes problems as @@ -410,6 +410,21 @@ public class CReferenceTest { // there is no imaginary component and real is negative. // The computation should be done for real only numbers separately. + // Check overflow safe. + final double a = Double.MAX_VALUE; + final double b = a / 4; + Assertions.assertEquals(inf, Complex.ofCartesian(a, b).abs(), "Expected overflow"); + // Compute the expected new magnitude by expressing b as a scale factor of a: + // (x^2 + y^2)^0.25 + // = sqrt(sqrt(a^2 + (b/a)^2 * a^2)) + // = sqrt(sqrt((1+(b/a)^2) * a^2)) + // = sqrt(sqrt((1+(b/a)^2))) * sqrt(a) + final double newAbs = Math.sqrt(Math.sqrt(1 + (b / a) * (b / a))) * Math.sqrt(a); + assertComplex(a, b, Complex::sqrt, newAbs * Math.cos(0.5 * Math.atan2(b, a)), + newAbs * Math.sin(0.5 * Math.atan2(b, a)), 3); + assertComplex(b, a, Complex::sqrt, newAbs * Math.cos(0.5 * Math.atan2(a, b)), + newAbs * Math.sin(0.5 * Math.atan2(a, b)), 2); + assertComplex(-2, 0.0, Complex::sqrt, 0, 1.4142135623730951); assertComplex(-2, 0.5, Complex::sqrt, 0.17543205637629397, 1.425053124063947, 5); assertComplex(-2, 1, Complex::sqrt, 0.3435607497225126, 1.4553466902253549, 3);
