You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2019/12/24 22:31:14 UTC
[commons-numbers] 06/09: Update log() to use method from Hull et al.
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 c283420af7165e09729002978ce195ec48fbab25
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Tue Dec 24 14:44:44 2019 +0000
Update log() to use method from Hull et al.
This requires changes to the edge case test as small values are computed
more accurately.
---
.../apache/commons/numbers/complex/Complex.java | 163 ++++++++++++++++-----
.../numbers/complex/ComplexEdgeCaseTest.java | 86 ++++++++++-
2 files changed, 212 insertions(+), 37 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 68dc607..a223f2a 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
@@ -73,8 +73,16 @@ public final class Complex implements Serializable {
private static final int MASK_INT_TO_EVEN = ~0x1;
/** Natural logarithm of 2 (ln(2)). */
private static final double LN_2 = Math.log(2);
+ /** Base 10 logarithm of 10 divided by 2 (log10(e)/2). */
+ private static final double LOG_10E_O_2 = Math.log10(Math.E) / 2;
/** Base 10 logarithm of 2 (log10(2)). */
private static final double LOG10_2 = Math.log10(2);
+ /** {@code 1/2}. */
+ private static final double HALF = 0.5;
+ /** {@code sqrt(2)}. */
+ private static final double ROOT2 = Math.sqrt(2);
+ /** The number of bits of precision of the mantissa of a {@code double} + 1: {@code 54}. */
+ private static final double PRECISION_1 = 54;
/**
* Crossover point to switch computation for asin/acos factor A.
@@ -2076,6 +2084,13 @@ public final class Complex implements Serializable {
* ln(a + b i) = ln(|a + b i|) + i arg(a + b i)
* </pre>
*
+ * <p>The implementation is based on the method described in:</p>
+ * <blockquote>
+ * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1994)
+ * Implementing complex elementary functions using exception handling.
+ * ACM Transactions on Mathematical Software, Vol 20, No 2, pp 215-244.
+ * </blockquote>
+ *
* @return the natural logarithm of {@code this}.
* @see Math#log(double)
* @see #abs()
@@ -2083,7 +2098,7 @@ public final class Complex implements Serializable {
* @see <a href="http://functions.wolfram.com/ElementaryFunctions/Log/">Log</a>
*/
public Complex log() {
- return log(Math::log, LN_2, Complex::ofCartesian);
+ return log(Math::log, HALF, LN_2, Complex::ofCartesian);
}
/**
@@ -2101,7 +2116,7 @@ public final class Complex implements Serializable {
* @see #arg()
*/
public Complex log10() {
- return log(Math::log10, LOG10_2, Complex::ofCartesian);
+ return log(Math::log10, LOG_10E_O_2, LOG10_2, Complex::ofCartesian);
}
/**
@@ -2116,51 +2131,131 @@ public final class Complex implements Serializable {
* will be incorrect. This is provided as an internal optimisation.
*
* @param log Log function.
+ * @param logOfeOver2 The log function applied to e, then divided by 2.
* @param logOf2 The log function applied to 2.
* @param constructor Constructor for the returned complex.
* @return the logarithm of {@code this}.
* @see #abs()
* @see #arg()
*/
- private Complex log(UnaryOperation log, double logOf2, ComplexConstructor constructor) {
- // TODO - Add edge cases with values for abs() close to 1
- // These should be handled correctly.
+ private Complex log(UnaryOperation log, double logOfeOver2, double logOf2, ComplexConstructor constructor) {
+ // Handle NaN
+ if (Double.isNaN(real) || Double.isNaN(imaginary)) {
+ // Return NaN unless infinite
+ if (isInfinite()) {
+ return constructor.create(Double.POSITIVE_INFINITY, Double.NaN);
+ }
+ // No-use of the input constructor
+ return NAN;
+ }
- // Change this to use a method based on glibc with detection of values of
- // abs() close to 1 and switch to using log1p(abs - 1).
+ // Compute the real part:
+ // log(sqrt(x^2 + y^2))
+ // log(x^2 + y^2) / 2
- // All ISO C99 edge cases satisfied by the Math library.
- // Make computation overflow safe.
+ // Compute with positive values
+ double x = Math.abs(real);
+ double y = Math.abs(imaginary);
- // Note:
- // log(|a + b i|) = log(sqrt(a^2 + b^2)) = 0.5 * log(a^2 + b^2)
- // If real and imaginary are with a safe region then omit the sqrt().
- final double x = Math.abs(real);
- final double y = Math.abs(imaginary);
+ // Find the larger magnitude.
+ if (x < y) {
+ double tmp = x;
+ x = y;
+ y = tmp;
+ }
- // Use the safe region defined for atanh to avoid over/underflow for x^2
- if (inRegion(x, y, SAFE_LOWER, SAFE_UPPER)) {
- //return constructor.create(log.apply(abs()), arg());
- return constructor.create(0.5 * log.apply(x * x + y * y), arg());
+ if (x == 0) {
+ // Handle zero: raises the ‘‘divide-by-zero’’ floating-point exception.
+ return constructor.create(Double.NEGATIVE_INFINITY,
+ negative(real) ? Math.copySign(Math.PI, imaginary) : imaginary);
}
- final double abs = abs();
- if (abs == Double.POSITIVE_INFINITY && isFinite()) {
- // Edge-case where the |a + b i| overflows.
- // |a + b i| = sqrt(a^2 + b^2)
- // This can be scaled linearly.
- // Scale the absolute and exploit:
- // ln(abs / scale) = ln(abs) - ln(scale)
- // ln(abs) = ln(abs / scale) + ln(scale)
- // Use precise scaling with:
- // scale ~ 2^exponent
- final int exponent = getMaxExponent(real, imaginary);
- // Implement scaling using 2^-exponent
- final double absOs = Math.hypot(Math.scalb(real, -exponent), Math.scalb(imaginary, -exponent));
- // log(2^exponent) = ln2(2^exponent) * log(2)
- return constructor.create(log.apply(absOs) + exponent * logOf2, arg());
+ double re;
+
+ if (x > HALF && x < ROOT2) {
+ // x^2+y^2 close to 1. Use log1p(x^2+y^2 - 1) / 2.
+ re = Math.log1p(x2y2m1(x, y)) * logOfeOver2;
+ } else if (y == 0) {
+ // Handle real only number
+ re = log.apply(x);
+ } else if (x > SAFE_MAX || x < SAFE_MIN || y < SAFE_MIN) {
+ // Over/underflow of sqrt(x^2+y^2)
+ if (isPosInfinite(x)) {
+ // Handle infinity
+ re = x;
+ } else {
+ // Do scaling
+ int expx = Math.getExponent(x);
+ int expy = Math.getExponent(y);
+ if (2 * (expx - expy) > PRECISION_1) {
+ // y can be ignored
+ re = log.apply(x);
+ } else {
+ // Hull et al:
+ // "It is important that the scaling be chosen so
+ // that there is no possibility of cancellation in this addition"
+ // i.e. sx^2 + sy^2 should not be close to 1.
+ // Their paper uses expx + 2 for underflow but expx for overflow.
+ // It has been modified here to use expx - 2.
+ int scale;
+ if (x > SAFE_MAX) {
+ // overflow
+ scale = expx - 2;
+ } else {
+ // underflow
+ scale = expx + 2;
+ }
+ double sx = Math.scalb(x, -scale);
+ double sy = Math.scalb(y, -scale);
+ re = scale * logOf2 + 0.5 * log.apply(sx * sx + sy * sy);
+ }
+ }
+ } else {
+ // Safe region that avoids under/overflow
+ re = 0.5 * log.apply(x * x + y * y);
+ }
+
+ // All ISO C99 edge cases for the imaginary are satisfied by the Math library.
+ return constructor.create(re, arg());
+ }
+
+ /**
+ * Compute {@code x^2 + y^2 - 1} in high precision.
+ * Assumes that the values x and y can be multiplied without overflow; that
+ * {@code x >= y}; and both values are positive.
+ *
+ * @param x the x value
+ * @param y the y value
+ * @return {@code x^2 + y^2 - 1}
+ */
+ private static double x2y2m1(double x, double y) {
+ // Hull et al used (x-1)*(x+1)+y*y.
+ // From the paper on page 236:
+
+ // If x == 1 there is no cancellation.
+
+ // If x > 1, there is also no cancellation, but the argument is now accurate
+ // only to within a factor of 1 + 3 EPSILSON (note that x – 1 is exact),
+ // so that error = 3 EPSILON.
+
+ // If x < 1, there can be serious cancellation:
+
+ // If 4 y^2 < |x^2 – 1| the cancellation is not serious ... the argument is accurate
+ // only to within a factor of 1 + 4 EPSILSON so that error = 4 EPSILON.
+
+ // 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;
}
- return constructor.create(log.apply(abs), arg());
+ return xm1xp1 + yy;
}
/**
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 86d8e00..df20dc9 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
@@ -20,6 +20,7 @@ package org.apache.commons.numbers.complex;
import org.junit.jupiter.api.Assertions;
import org.junit.jupiter.api.Test;
+import java.math.BigDecimal;
import java.util.function.BiFunction;
import java.util.function.UnaryOperator;
@@ -344,9 +345,88 @@ public class ComplexEdgeCaseTest {
assertComplex(-Double.MAX_VALUE, Double.MAX_VALUE / 4, name, operation, 7.098130252042921e2, 2.896613990462929);
assertComplex(Double.MAX_VALUE, Double.MAX_VALUE / 4, name, operation, 7.098130252042921e2, 2.449786631268641e-1, 2);
- // Underflow if sqrt(a^2 + b^2) == 0
- assertComplex(-Double.MIN_VALUE, Double.MIN_VALUE, name, operation, -744.44007192138122, 2.3561944901923448);
- assertComplex(Double.MIN_VALUE, Double.MIN_VALUE, name, operation, -744.44007192138122, 0.78539816339744828);
+ // Underflow if sqrt(a^2 + b^2) -> 0
+ assertComplex(-Double.MIN_NORMAL, Double.MIN_NORMAL, name, operation, -708.04984494198413, 2.3561944901923448);
+ assertComplex(Double.MIN_NORMAL, Double.MIN_NORMAL, name, operation, -708.04984494198413, 0.78539816339744828);
+ // Math.hypot(min, min) = min.
+ // To compute the expected result do scaling of the actual hypot = sqrt(2).
+ // log(a/n) = log(a) - log(n)
+ // n = 2^1074 => log(a) - log(2) * 1074
+ double expected = Math.log(Math.sqrt(2)) - Math.log(2) * 1074;
+ assertComplex(-Double.MIN_VALUE, Double.MIN_VALUE, name, operation, expected, Math.atan2(1, -1));
+ expected = Math.log(Math.sqrt(5)) - Math.log(2) * 1074;
+ assertComplex(-Double.MIN_VALUE, 2 * Double.MIN_VALUE, name, operation, expected, Math.atan2(2, -1));
+
+ // Imprecision if sqrt(a^2 + b^2) == 1 as log(1) is 0.
+ // Method should switch to using log1p(x^2 + x^2 - 1) * 0.5.
+
+ // In the following:
+ // max = max(real, imaginary)
+ // min = min(real, imaginary)
+
+ // 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);
+
+ // 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);
+
+ // Very hard: 4 * min^2 > |max^2 - 1|
+
+ // Radius 0.99
+ assertLog(0.97, Math.sqrt(0.99 - 0.97 * 0.97), 30);
+ // Radius 1.01
+ assertLog(0.97, Math.sqrt(1.01 - 0.97 * 0.97), 34);
+
+ // 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);
+ }
+
+ /**
+ * Assert the Complex log function using BigDecimal to compute the field norm
+ * {@code x*x + y*y} and then {@link Math#log1p(double)} to compute the log of
+ * the modulus \ using {@code 0.5 * log1p(x*x + y*y - 1)}. This test is for the
+ * extreme case for performance around {@code sqrt(x*x + y*y) = 1} where using
+ * {@link Math#log(double)} will fail dramatically.
+ *
+ * @param x the real value of the complex
+ * @param y the imaginary value of the complex
+ * @param maxUlps the maximum units of least precision between the two values
+ */
+ private static void assertLog(double x, double y, long maxUlps) {
+ // Compute the best value we can
+ BigDecimal bx = BigDecimal.valueOf(x);
+ BigDecimal by = BigDecimal.valueOf(y);
+ double real = 0.5 * Math.log1p(bx.multiply(bx).add(by.multiply(by)).subtract(BigDecimal.ONE).doubleValue());
+ double imag = Math.atan2(y, x);
+ assertComplex(x, y, "log", Complex::log, real, imag, maxUlps);
}
@Test