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 f61e1136f57db2429182c64e3bf2ef9fef3712f6
Author: Alex Herbert <aherb...@apache.org>
AuthorDate: Thu Jan 2 15:03:47 2020 +0000
Use high precision computation in Complex::log.
This avoids the enormous relative error documented by Hull et al when 4
* y^2 >= |x^2 – 1|.
---
.../apache/commons/numbers/complex/Complex.java | 138 +++++++++++++++++++--
.../numbers/complex/ComplexEdgeCaseTest.java | 77 +++++++-----
pom.xml | 5 +
3 files changed, 179 insertions(+), 41 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 202cc1d..fbf0777 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
@@ -99,6 +99,10 @@ public final class Complex implements Serializable {
* @see <a href="http://en.wikipedia.org/wiki/Machine_epsilon";>Machine
epsilon</a>
*/
private static final double EPSILON =
Double.longBitsToDouble((EXPONENT_OFFSET - 53L) << 52);
+ /** The multiplier used to split the double value into hi and low parts.
This must be odd
+ * and a value of 2^s + 1 in the range {@code p/2 <= s <= p-1} where p is
the number of
+ * bits of precision of the floating point number. Here {@code s = 27}.*/
+ private static final double MULTIPLIER = (1 << 27) + 1.0;
/**
* Crossover point to switch computation for asin/acos factor A.
@@ -2221,16 +2225,132 @@ public final class Complex implements Serializable {
// Otherwise there can be serious cancellation and the relative error
in the real part
// could be enormous.
- // TODO - investigate the computation in high precision using
- // LinearCombination#value(double, double, double, double, double,
double)
- // from o.a.c.numbers.arrays.
- final double xm1xp1 = (x - 1) * (x + 1);
final double yy = y * y;
- if (x < 1 && 4 * yy > Math.abs(xm1xp1)) {
- // Large relative error...
- return xm1xp1 + yy;
+ // Modify to use high precision before the threshold set by Hull et al.
+ // This is to preserve the monotonic output of the computation at the
switch.
+ // Set the threshold when x^2 + y^2 is above 0.5 thus subtracting 1
results in a number
+ // that can be expressed with a higher precision than any number in
the range 0.5-1.0
+ // due to the variable exponent used below 0.5.
+ if (x < 1 && x * x + yy > 0.5) {
+ // Large relative error.
+ // This does not use o.a.c.numbers.LinearCombination.value(x, x,
y, y, 1, -1).
+ // It is optimised knowing that:
+ // - the products are squares
+ // - the final term is -1 (which does not require split
multiplication and addition)
+ // - The answer will not be NaN as the terms are not NaN components
+ // - The order is known to be 1 > |x| >= |y|
+ // The squares are computed using a split multiply algorithm and
+ // the summation using a fast-expansion-sum algorithm.
+
+ // Split x and y as one 26 bits number and one 27 bits number
+ final double xHigh = splitHigh(x);
+ final double xLow = x - xHigh;
+ final double yHigh = splitHigh(y);
+ final double yLow = y - yHigh;
+
+ // Accurate split multiplication x * x and y * y
+ final double x2High = x * x;
+ final double x2Low = squareRoundOff(xLow, xHigh, x2High);
+ final double y2High = y * y;
+ final double y2Low = squareRoundOff(yLow, yHigh, y2High);
+
+ return fastExpansionSumx2y2m1(x2High, x2Low, y2High, y2Low);
}
- return xm1xp1 + yy;
+ return (x - 1) * (x + 1) + yy;
+ }
+
+ /**
+ * Implement Dekker's method to split a value into two parts. Multiplying
by (2^s + 1) create
+ * a big value from which to derive the two split parts.
+ * <pre>
+ * c = (2^s + 1) * a
+ * a_big = c - a
+ * a_hi = c - a_big
+ * a_lo = a - a_hi
+ * a = a_hi + a_lo
+ * </pre>
+ *
+ * <p>The multiplicand must be odd allowing a p-bit value to be split into
+ * (p-s)-bit value {@code a_hi} and a nonoverlapping (s-1)-bit value
{@code a_lo}.
+ * Combined they have (p-1) bits of significand but the sign bit of
{@code a_lo}
+ * contains a bit of information.
+ *
+ * @param a Value.
+ * @return the high part of the value.
+ * @see <a
href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps";>
+ * Shewchuk (1997) Theorum 17</a>
+ */
+ private static double splitHigh(double a) {
+ final double c = MULTIPLIER * a;
+ return c - (c - a);
+ }
+
+ /**
+ * Compute the round-off from the square of a split number with {@code
low} and {@code high}
+ * components. Uses Dekker's algorithm for split multiplication modified
for a square product.
+ *
+ * @param low Low part of number.
+ * @param high High part of number.
+ * @param square Square of the number.
+ * @return <code>low * low - ((product - high * high) - 2 * low *
high)</code>
+ */
+ private static double squareRoundOff(double low, double high, double
square) {
+ return low * low - ((square - high * high) - 2 * low * high);
+ }
+
+ /**
+ * Compute the round-off from the sum of two numbers {@code a} and {@code
b} using
+ * Knuth's two-sum algorithm. The values are not required to be ordered by
magnitude.
+ *
+ * @param a First part of sum.
+ * @param b Second part of sum.
+ * @param sum Sum.
+ * @return <code>(b - (sum - (sum - b))) + (a - (sum - b))</code>
+ */
+ private static double sumRoundOff(double a, double b, double sum) {
+ final double aPrime = sum - b;
+ // sum - aPrime == bPrime.
+ // a - aPrime == a round-off
+ // b - bPrime == b round-off
+ return (a - aPrime) + (b - (sum - aPrime));
+ }
+
+ /**
+ * Sum x^2 + y^2 - 1.
+ *
+ * <p>Uses Shewchuk's fast expansion sum under the assumption that {@code
1 > |x| >= |y|}
+ * thus the high and low products from {@code x^2} and {@code y^2} can be
ordered into
+ * a single sequence g, in order of nondecreasing magnitude with no
sorting.
+ *
+ * @param x2High High part of x^2.
+ * @param x2Low Low part of x^2.
+ * @param y2High High part of y^2.
+ * @param y2Low Low part of y^2.
+ * @return x^2 + y^2 - 1
+ * @see <a
href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps";>
+ * Shewchuk (1997) Adaptive Precision Floating-Point Arithmetic and Fast
Robust Geometric Predicates: Theorum 13</a>
+ */
+ private static double fastExpansionSumx2y2m1(double x2High, double x2Low,
double y2High, double y2Low) {
+ // Create a single sequence g that is in order of nondecreasing
magnitude.
+ // Since |x| >= |y| then x2High >= y2High and the two smallest
components are the
+ // round-off from x^2 and y^2. The final term is -1 which will never
have a bit overlap
+ // with the products x^2 or y^2 as 1 has no bits.
+ // Thus this implements fast-two-sum: e + f as {x2Low, x2High} +
{y2Low, y2High, -1}.
+ // In Shewchuk the two smallest components are added using a
+ // fast-two-sum which requires ordering. Here we use a two-sum to
avoid order comparators.
+ // The remaining components are added with a two-sum that does not
require ordering.
+ // Q_i is the running total of summing sequence g, {y2Low, x2Low,
y2High, x2High, -1}.
+ // h_i is the output sequence; this is guaranteed to produce a
strongly non-overlapping
+ // output if its inputs are strongly non-overlapping. The sum of h is
the sum of e + f.
+ final double q1 = x2Low + y2Low;
+ final double h1 = sumRoundOff(x2Low, y2Low, q1);
+ final double q2 = q1 + y2High;
+ final double h2 = sumRoundOff(q1, y2High, q2);
+ final double q3 = q2 + x2High;
+ final double h3 = sumRoundOff(q2, x2High, q3);
+ final double h5 = q3 - 1;
+ final double h4 = sumRoundOff(q3, -1, h5);
+ return h1 + h2 + h3 + h4 + h5;
}
/**
@@ -2534,8 +2654,6 @@ public final class Complex implements Serializable {
return tanh(-imaginary, real, Complex::multiplyNegativeI);
}
- // TODO
-
/**
* Returns the
* <a href="http://mathworld.wolfram.com/HyperbolicTangent.html";>
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 90ca920..b736f86 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
@@ -367,46 +367,60 @@ public class ComplexEdgeCaseTest {
// No cancellation error when max > 1
assertLog(1.0001, Math.sqrt(1.2 - 1.0001 * 1.0001), 1);
- assertLog(1.0001, Math.sqrt(1.1 - 1.0001 * 1.0001), 2);
- assertLog(1.0001, Math.sqrt(1.02 - 1.0001 * 1.0001), 6);
+ assertLog(1.0001, Math.sqrt(1.1 - 1.0001 * 1.0001), 1);
+ assertLog(1.0001, Math.sqrt(1.02 - 1.0001 * 1.0001), 0);
+ assertLog(1.0001, Math.sqrt(1.01 - 1.0001 * 1.0001), 0);
// Cancellation error when max < 1.
- // The ULPs may need to be revised if the log() method is improved.
// Hard: 4 * min^2 < |max^2 - 1|
- assertLog(0.9, 0.00001, 3);
- assertLog(0.95, 0.00001, 8);
+ // Gets harder as max is further from 1
+ assertLog(0.99, 0.00001, 0);
+ assertLog(0.95, 0.00001, 0);
+ assertLog(0.9, 0.00001, 0);
+ assertLog(0.85, 0.00001, 0);
+ assertLog(0.8, 0.00001, 0);
+ assertLog(0.75, 0.00001, 0);
+ // At this point the log function does not use high precision
computation
+ assertLog(0.7, 0.00001, 2);
// Very hard: 4 * min^2 > |max^2 - 1|
// Radius 0.99
- assertLog(0.97, Math.sqrt(0.99 - 0.97 * 0.97), 30);
+ assertLog(0.97, Math.sqrt(0.99 - 0.97 * 0.97), 0);
// Radius 1.01
- assertLog(0.97, Math.sqrt(1.01 - 0.97 * 0.97), 34);
+ assertLog(0.97, Math.sqrt(1.01 - 0.97 * 0.97), 0);
// Massive relative error
// Radius 0.9999
- assertLog(0.97, Math.sqrt(0.9999 - 0.97 * 0.97), 3639);
-
- // This code is for demonstration purposes.
-
-// // Demonstrate relative error using
-// // cis numbers on a 1/8 circle with a set radius.
-// int steps = 10;
-// for (double radius : new double[] {0.99, 1.0, 1.01}) {
-// for (int i = 1; i < steps; i++) {
-// double theta = i * Math.PI / (8 * steps);
-// assertLog(radius * Math.sin(theta), radius *
Math.cos(theta), -1);
-// }
-// }
-//
-// // Extreme. Here for documentation of the relative error.
-// double up1 = Math.nextUp(1.0);
-// double down1 = Math.nextDown(1.0);
-// assertLog(up1, Double.MIN_NORMAL, -1);
-// assertLog(up1, Double.MIN_VALUE, -1);
-// assertLog(down1, Double.MIN_NORMAL, -1);
-// assertLog(down1, Double.MIN_VALUE, -1);
+ assertLog(0.97, Math.sqrt(0.9999 - 0.97 * 0.97), 0);
+
+ // cis numbers on a 1/8 circle with a set radius.
+ final int steps = 20;
+ final double[] radius = {0.99, 1.0, 1.01};
+ final int[] ulps = {0, 1, 1};
+ for (int j = 0; j < radius.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]);
+ }
+ }
+
+ // cis numbers using an increasingly smaller angle
+ double theta = Math.PI / (4 * steps);
+ while (theta > 0) {
+ theta /= 2;
+ assertLog(Math.sin(theta), Math.cos(theta), 1);
+ }
+
+ // Extreme cases.
+ final double up1 = Math.nextUp(1.0);
+ final double down1 = Math.nextDown(1.0);
+ assertLog(down1, Double.MIN_NORMAL, 0);
+ assertLog(down1, Double.MIN_VALUE, 0);
+ // No use of high-precision computation
+ assertLog(up1, Double.MIN_NORMAL, 2);
+ assertLog(up1, Double.MIN_VALUE, 2);
}
/**
@@ -422,9 +436,10 @@ public class ComplexEdgeCaseTest {
*/
private static void assertLog(double x, double y, long maxUlps) {
// Compute the best value we can
- final BigDecimal bx = BigDecimal.valueOf(x);
- final BigDecimal by = BigDecimal.valueOf(y);
- final double real = 0.5 *
Math.log1p(bx.multiply(bx).add(by.multiply(by)).subtract(BigDecimal.ONE).doubleValue());
+ final BigDecimal bx = new BigDecimal(x);
+ final BigDecimal by = new BigDecimal(y);
+ final BigDecimal exact =
bx.multiply(bx).add(by.multiply(by)).subtract(BigDecimal.ONE);
+ final double real = 0.5 * Math.log1p(exact.doubleValue());
final double imag = Math.atan2(y, x);
assertComplex(x, y, "log", Complex::log, real, imag, maxUlps);
}
diff --git a/pom.xml b/pom.xml
index a1fb00f..eb32172 100644
--- a/pom.xml
+++ b/pom.xml
@@ -152,6 +152,11 @@
</dependency>
<dependency>
<groupId>org.apache.commons</groupId>
+ <artifactId>commons-numbers-arrays</artifactId>
+ <version>${project.version}</version>
+ </dependency>
+ <dependency>
+ <groupId>org.apache.commons</groupId>
<artifactId>commons-numbers-core</artifactId>
<version>${project.version}</version>
<type>test-jar</type>
