You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2019/12/09 00:03:43 UTC
[commons-numbers] 05/09: Compute extreme values of sqrt() in polar
coordinates.
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 <ah...@apache.org>
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);