commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From aherb...@apache.org
Subject [commons-numbers] branch master updated: Improve documentation of the high precision sum in Complex.
Date Sat, 22 Feb 2020 00:04:18 GMT
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


The following commit(s) were added to refs/heads/master by this push:
     new 2b13b8f  Improve documentation of the high precision sum in Complex.
2b13b8f is described below

commit 2b13b8f7b98571d07eba452c157d0d98263990cc
Author: Alex Herbert <aherbert@apache.org>
AuthorDate: Sat Feb 22 00:04:13 2020 +0000

    Improve documentation of the high precision sum in Complex.
    
    Added note to state that the current method is consistently more
    accurate than a faster alternative suggested by Shewchuk. To achieve the
    same consistency requires a distillation sum of the output expansion of
    fast-two-sum.
---
 .../apache/commons/numbers/complex/Complex.java    | 53 +++++++++++++++-------
 .../numbers/complex/ComplexEdgeCaseTest.java       | 10 ++--
 2 files changed, 41 insertions(+), 22 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 cae5385..559d7e8 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
@@ -2522,7 +2522,7 @@ public final class Complex implements Serializable  {
     /**
      * Sum x^2 + y^2 - 1. It is assumed that {@code y <= x < 1}.
      *
-     * <p>Implement Shewchuk's expansion-sum algorithm: [x2Low, x2High, -1] + [y2Low,
y2High].
+     * <p>Implement Shewchuk's expansion-sum algorithm: [x2Low, x2High] + [-1] + [y2Low,
y2High].
      *
      * @param x2High High part of x^2.
      * @param x2Low Low part of x^2.
@@ -2537,10 +2537,10 @@ public final class Complex implements Serializable  {
         // The following algorithm will produce a non-overlapping expansion h where the
         // sum h_i = e + f and components of h are in increasing order of magnitude.
 
-        // Expansion sum proceeds by a grow-expansion of the first part from one expansion
+        // Expansion-sum proceeds by a grow-expansion of the first part from one expansion
         // into the other, extending its length by 1. The process repeats for the next part
-        // but the grow expansion starts at the previous merge position + 1.
-        // Thus expansion sum requires mn two-sum operations to merge length m into length
n
+        // but the grow-expansion starts at the previous merge position + 1.
+        // Thus expansion-sum requires mn two-sum operations to merge length m into length
n
         // resulting in length m+n-1.
 
         // Variables numbered from 1 as per Figure 7 (p.12). The output expansion h is placed
@@ -2549,7 +2549,21 @@ public final class Complex implements Serializable  {
         // We have two expansions for x^2 and y^2 and the whole number -1.
         // Expecting (x^2 + y^2) close to 1 we generate first the intermediate expansion
         // (x^2 - 1) moving the result away from 1 where there are sparse floating point
-        // representations. This is then added to a similar magnitude y^2.
+        // representations. This is then added to a similar magnitude y^2. Leaving the -1
+        // until last suffers from 1 ulp rounding errors more often and the requirement
+        // for a distillation sum to reduce rounding error frequency.
+
+        // Note: Do not use the alternative fast-expansion-sum of the parts sorted by magnitude.
+        // The parts can be ordered with a single comparison into:
+        // [y2Low, (y2High|x2Low), x2High, -1]
+        // The fast-two-sum saves 1 fast-two-sum and 3 two-sum operations (21 additions)
and
+        // adds a penalty of a single branch condition.
+        // However the order in not "strongly non-overlapping" and the fast-expansion-sum
+        // output will not be strongly non-overlapping. The sum of the output has 1 ulp error
+        // on random cis numbers approximately 1 in 160 events. This can be removed by a
+        // distillation two-sum pass over the final expansion as a cost of 1 fast-two-sum
and
+        // 3 two-sum operations! So we use the expansion sum with the same operations and
+        // no branches.
 
         // q=running sum
         double q = x2Low - 1;
@@ -2563,22 +2577,27 @@ public final class Complex implements Serializable  {
         // Grow expansion of f1 into e
         q = f1 + e1;
         e1 = sumLow(f1, e1, q);
-        double p = q;
-        q += e2;
-        e2 = sumLow(p, e2, q);
-        double e4 = q + e3;
-        e3 = sumLow(q, e3, e4);
+        double p = q + e2;
+        e2 = sumLow(q, e2, p);
+        double e4 = p + e3;
+        e3 = sumLow(p, e3, e4);
 
         // Grow expansion of f2 into e (only required to start at e2)
         q = f2 + e2;
         e2 = sumLow(f2, e2, q);
-        p = q;
-        q += e3;
-        e3 = sumLow(p, e3, q);
-        final double e5 = q + e4;
-        e4 = sumLow(q, e4, e5);
-
-        // Final summation
+        p = q + e3;
+        e3 = sumLow(q, e3, p);
+        final double e5 = p + e4;
+        e4 = sumLow(p, e4, e5);
+
+        // Final summation:
+        // The sum of the parts is within 1 ulp of the true expansion value e:
+        // |e - sum| < ulp(sum).
+        // To achieve the exact result requires iteration of a distillation two-sum through
+        // the expansion until convergence, i.e. no smaller term changes higher terms.
+        // This requires (n-1) iterations for length n. Here we neglect this as
+        // although the method is not ensured to be exact is it robust on random
+        // cis numbers.
         return e1 + e2 + e3 + e4 + e5;
     }
 
diff --git a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexEdgeCaseTest.java
b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexEdgeCaseTest.java
index ac20b48..4f76f4d 100644
--- a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexEdgeCaseTest.java
+++ b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexEdgeCaseTest.java
@@ -489,14 +489,14 @@ public class ComplexEdgeCaseTest {
         // Radius 0.9999
         assertLog(0.97, Math.sqrt(0.9999 - 0.97 * 0.97), 0);
 
-        // cis numbers on a 1/8 circle with a set radius.
+        // polar numbers on a 1/8 circle with a magnitude close to 1.
         final int steps = 20;
-        final double[] radius = {0.99, 1.0, 1.01};
+        final double[] magnitude = {0.999, 1.0, 1.001};
         final int[] ulps = {0, 0, 1};
-        for (int j = 0; j < radius.length; j++) {
+        for (int j = 0; j < magnitude.length; j++) {
             for (int i = 1; i <= steps; i++) {
                 final double theta = i * Math.PI / (4 * steps);
-                assertLog(radius[j] * Math.sin(theta), radius[j] * Math.cos(theta), ulps[j]);
+                assertLog(magnitude[j] * Math.sin(theta), magnitude[j] * Math.cos(theta),
ulps[j]);
             }
         }
 
@@ -504,7 +504,7 @@ public class ComplexEdgeCaseTest {
         double theta = Math.PI / (4 * steps);
         while (theta > 0) {
             theta /= 2;
-            assertLog(Math.sin(theta), Math.cos(theta), 1);
+            assertLog(Math.sin(theta), Math.cos(theta), 0);
         }
 
         // Extreme cases.


Mime
View raw message