[commons-numbers] 02/04: Use high precision x^2 + y^2 - 1 in atanh.

[ah...@apache.org]

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 <ah...@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));
                 }