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 0cfc94de303b6cc56c3dbc6157db8d70bfa8c975 Author: Alex Herbert <[email protected]> AuthorDate: Sun Dec 8 23:18:40 2019 +0000 Compute extreme values of sqrt() in polar coordinates. --- .../apache/commons/numbers/complex/Complex.java | 22 ++++++++++++++++------ .../commons/numbers/complex/CReferenceTest.java | 12 +++++++++++- 2 files changed, 27 insertions(+), 7 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 2d88c2f..9f9ac9e 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 @@ -1711,13 +1711,23 @@ public final class Complex implements Serializable { if (real == 0 && imaginary == 0) { return new Complex(0, imaginary); } - final double t = Math.sqrt(average(Math.abs(real), Math.hypot(real, imaginary))); - // t is positive: - // To prevent overflow dividing by (2 * t) divide by t then by 2. + final double abs = getAbsolute(real, imaginary); + final double av = average(Math.abs(real), abs); + 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); + } + + final double t = Math.sqrt(av); if (real >= 0) { - return new Complex(t, (imaginary / t) / 2); + return new Complex(t, imaginary / (2 * t)); } - return new Complex((Math.abs(imaginary) / t) / 2, + return new Complex(Math.abs(imaginary) / (2 * t), Math.copySign(1.0, imaginary) * t); } // Imaginary is nan @@ -1744,7 +1754,7 @@ public final class Complex implements Serializable { * @return the average */ private static double average(double a, double b) { - return (b < a) ? + return (a < b) ? a + (b - a) / 2 : b + (a - b) / 2; } 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 2dbf394..22b8cf2 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 @@ -399,18 +399,28 @@ public class CReferenceTest { @Test public void testSqrt() { + // TODO: GNU g++ computes in polar coordinates. + // Test other reference implementations. + // This gives strange results for real negative and imaginary 0. + + // Polar result + //assertComplex(-2, 0.0, Complex::sqrt, 8.6595605623549341e-17, 1.4142135623730951); 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); assertComplex(-2, 2, Complex::sqrt, 0.64359425290558281, 1.5537739740300374, 2); + // Polar result + //assertComplex(-1, 0.0, Complex::sqrt, 6.123233995736766e-17, 1); assertComplex(-1, 0.0, Complex::sqrt, 0, 1); assertComplex(-1, 0.5, Complex::sqrt, 0.24293413587832291, 1.0290855136357462, 3); assertComplex(-1, 1, Complex::sqrt, 0.45508986056222739, 1.0986841134678098); assertComplex(-1, 2, Complex::sqrt, 0.78615137775742339, 1.2720196495140688); + // Polar result + //assertComplex(-0.5, 0.0, Complex::sqrt, 4.329780281177467e-17, 0.70710678118654757); assertComplex(-0.5, 0.0, Complex::sqrt, 0, 0.70710678118654757); assertComplex(-0.5, 0.5, Complex::sqrt, 0.3217971264527914, 0.77688698701501868, 2); assertComplex(-0.5, 1, Complex::sqrt, 0.55589297025142126, 0.89945371997393353); - assertComplex(-0.5, 2, Complex::sqrt, 0.88361553087551337, 1.1317139242778693); + assertComplex(-0.5, 2, Complex::sqrt, 0.88361553087551337, 1.1317139242778693, 2); assertComplex(-0.0, 0.0, Complex::sqrt, 0.0, 0.0); assertComplex(-0.0, 0.5, Complex::sqrt, 0.50000000000000011, 0.5); assertComplex(-0.0, 1, Complex::sqrt, 0.70710678118654757, 0.70710678118654746);
