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 ac7e98d29be8180a862fe6e2878775474b57b17f Author: Alex Herbert <aherb...@apache.org> AuthorDate: Sun Jan 5 22:48:46 2020 +0000 Use high precision x^2 + y^2 - 1 in atanh. --- .../apache/commons/numbers/complex/Complex.java | 29 ++++++++++++++++++---- 1 file changed, 24 insertions(+), 5 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 fbf0777..4b4d98b 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 @@ -1642,8 +1642,8 @@ public final class Complex implements Serializable { private static Complex atanh(final double real, final double imaginary, final ComplexConstructor constructor) { // Compute with positive values and determine sign at the end - final double x = Math.abs(real); - final double y = Math.abs(imaginary); + double x = Math.abs(real); + double y = Math.abs(imaginary); // The result (without sign correction) double re; double im; @@ -1680,7 +1680,7 @@ public final class Complex implements Serializable { // minus x plus 1: (-x+1) final double mxp1 = 1 - x; - final double yy = y * y; + double yy = y * y; // The definition of real component is: // real = log( ((x+1)^2+y^2) / ((1-x)^2+y^2) ) / 4 // This simplifies by adding 1 and subtracting 1 as a fraction: @@ -1688,9 +1688,28 @@ public final class Complex implements Serializable { // // real(atanh(z)) == log(1 + 4*x / ((1-x)^2+y^2)) / 4 // imag(atanh(z)) == tan^-1 (2y, (1-x)(1+x) - y^2) / 2 + // imag(atanh(z)) == tan^-1 (2y, (1 - x^2 - y^2) / 2 // The division is done at the end of the function. re = Math.log1p(4 * x / (mxp1 * mxp1 + yy)); - im = Math.atan2(2 * y, mxp1 * (1 + x) - yy); + // Modified from boost which does not switch the magnitude of x and y. + // The denominator for atan2 is 1 - x^2 - y^2. + // This can be made more precise if |x| > |y|. + final double numerator = 2 * y; + double denominator; + if (x < y) { + final double tmp = x; + x = y; + y = tmp; + } + // 1 - x is precise if |x| >= 1 + if (x >= 1) { + denominator = (1 - x) * (1 + x) - y * y; + } else { + // |x| < 1: Use high precision if possible: + // 1 - x^2 - y^2 = -(x^2 + y^2 - 1) + denominator = -x2y2m1(x, y); + } + im = Math.atan2(numerator, denominator); } else { // This section handles exception cases that would normally cause // underflow or overflow in the main formulas. @@ -1772,7 +1791,7 @@ public final class Complex implements Serializable { } else { // Medium x, small y. // Modified from boost which checks (y == 0) && (x == 1) and sets re = 0. - // This is same as the result from calling atan2(0, 0) so just do that. + // This is same as the result from calling atan2(0, 0) so exclude this case. // 1 - y^2 = 1 so ignore subtracting y^2 im = Math.atan2(2 * y, (1 - x) * (1 + x)); }