[commons-numbers] 02/07: Reciprocal to use divide(Complex).

[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 610e70ca7174067b047bbdb8a3a00675a4c70d89
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Fri Dec 27 21:30:29 2019 +0000

    Reciprocal to use divide(Complex).
    
    This has edge case recovery for infinity and NaN missing from the
    previous direct implementation.
    
    It also computes values when the previous method did not due to use of
    scaling. This allows computation for (Double.MAX_VALUE,
    Double.MAX_VALUE).reciprocal().
---
 .../apache/commons/numbers/complex/Complex.java    | 30 ++++++++------------
 .../commons/numbers/complex/ComplexTest.java       | 33 ++++++++++++++++------
 2 files changed, 35 insertions(+), 28 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 01f87cd..ec2d2ae 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
@@ -748,30 +748,22 @@ public final class Complex implements Serializable  {
     }
 
     /**
-     * Returns the multiplicative inverse of this instance.
+     * Returns the multiplicative inverse of this instance \( z \).
+     * \[ \frac{1}{z} = \frac{1}{a + i b} = \frac{a}{a^2+b^2} - i \frac{b}{a^2+b^2} \]
+     * This method produces the same result as:
+     * <pre>
+     *  Complex.ONE.divide(this)
+     * </pre>
      *
      * @return {@code 1 / this}.
      * @see #divide(Complex)
      */
     public Complex reciprocal() {
-        if (Math.abs(real) < Math.abs(imaginary)) {
-            final double q = real / imaginary;
-            final double scale = 1.0 / (real * q + imaginary);
-            double scaleQ = 0;
-            if (q != 0 &&
-                scale != 0) {
-                scaleQ = scale * q;
-            }
-            return new Complex(scaleQ, -scale);
-        }
-        final double q = imaginary / real;
-        final double scale = 1.0 / (imaginary * q + real);
-        double scaleQ = 0;
-        if (q != 0 &&
-            scale != 0) {
-            scaleQ = scale * q;
-        }
-        return new Complex(scale, -scaleQ);
+        // Note that this cannot be optimised assuming a=1 and b=0.
+        // These preserve compatibility with the divide method:
+        // 1. create NaNs for infinite parts c or d: 0.0 * inf = nan
+        // 2. maintain signs when either c or d are negative signed and the other part is zero
+        return divide(1.0, 0.0, real, imaginary);
     }
 
     /**
diff --git a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexTest.java b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexTest.java
index 3d1fd66..ad07c10 100644
--- a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexTest.java
+++ b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexTest.java
@@ -571,8 +571,8 @@ public class ComplexTest {
         final Complex act = z.reciprocal();
         final double expRe = 5.0 / 61.0;
         final double expIm = -6.0 / 61.0;
-        Assertions.assertEquals(expRe, act.getReal(), Math.ulp(expRe));
-        Assertions.assertEquals(expIm, act.getImaginary(), Math.ulp(expIm));
+        Assertions.assertEquals(expRe, act.getReal());
+        Assertions.assertEquals(expIm, act.getImaginary());
     }
 
     @Test
@@ -602,13 +602,28 @@ public class ComplexTest {
     }
 
     @Test
-    public void testReciprocalMax() {
-        // This hits the edge case in reciprocal() for when q != 0 but scale == 0
-        final double smaller = Math.nextDown(Double.MAX_VALUE);
-        Complex z = Complex.ofCartesian(smaller, Double.MAX_VALUE);
-        Assertions.assertEquals(Complex.ofCartesian(0.0, -0.0), z.reciprocal());
-        z = Complex.ofCartesian(Double.MAX_VALUE, smaller);
-        Assertions.assertEquals(Complex.ofCartesian(0.0, -0.0), z.reciprocal());
+    public void testReciprocalZero() {
+        final Complex z = Complex.ZERO.reciprocal();
+        Assertions.assertEquals(inf, z.getReal());
+        Assertions.assertEquals(nan, z.getImaginary());
+    }
+
+    @Test
+    public void testReciprocalMatchesDivide() {
+        final double[] parts = {Double.NEGATIVE_INFINITY, -1, -0.0, 0.0, 1, Math.PI, Double.POSITIVE_INFINITY, Double.NaN};
+        for (final double x : parts) {
+            for (final double y : parts) {
+                final Complex z = Complex.ofCartesian(x, y);
+                Assertions.assertEquals(Complex.ONE.divide(z), z.reciprocal(), () -> z.toString());
+            }
+        }
+        final UniformRandomProvider rng = RandomSource.create(RandomSource.SPLIT_MIX_64);
+        for (int i = 0; i < 10; i++) {
+            final double x = -1 + rng.nextDouble() * 2;
+            final double y = -1 + rng.nextDouble() * 2;
+            final Complex z = Complex.ofCartesian(x, y);
+            Assertions.assertEquals(Complex.ONE.divide(z), z.reciprocal());
+        }
     }
 
     @Test