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 2020/01/17 18:45:53 UTC
[commons-numbers] branch master updated (a79dcaa -> c4dc88c)
This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a change to branch master
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git.
from a79dcaa Change x2y2m1 summation algorithm.
new 62ee60b Modify hypotenuse scaling formula for x^2+y^2 in log().
new c4dc88c Ensure divide correctly scales numbers to avoid over/underflow.
The 2 revisions listed above as "new" are entirely new to this
repository and will be described in separate emails. The revisions
listed as "add" were already present in the repository and have only
been added to this reference.
Summary of changes:
.../apache/commons/numbers/complex/Complex.java | 96 ++++++++++++++++++++--
.../numbers/complex/ComplexEdgeCaseTest.java | 45 +++++++++-
2 files changed, 131 insertions(+), 10 deletions(-)
[commons-numbers] 02/02: Ensure divide correctly scales numbers to
avoid over/underflow.
Posted by 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 c4dc88cab80dd8c4448f43b2930373df9f92c118
Author: aherbert <ah...@apache.org>
AuthorDate: Fri Jan 17 18:45:48 2020 +0000
Ensure divide correctly scales numbers to avoid over/underflow.
Added tested to ensure the smallest and largest possible value divided
by itself is (1, 0). Previously it would compute zero for small and
overflow to infinite for large.
---
.../apache/commons/numbers/complex/Complex.java | 88 +++++++++++++++++++++-
.../numbers/complex/ComplexEdgeCaseTest.java | 45 ++++++++++-
2 files changed, 126 insertions(+), 7 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 8987031..d5a2ca7 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
@@ -112,6 +112,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);
+ /** Mask to remove the sign bit from a long. */
+ private static final long UNSIGN_MASK = 0x7fffffffffffffffL;
+ /** Mask to extract the 52-bit mantissa from a long representation of a double. */
+ private static final long MANTISSA_MASK = 0x000fffffffffffffL;
/** 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}.*/
@@ -646,14 +650,22 @@ public final class Complex implements Serializable {
double c = re2;
double d = im2;
int ilogbw = 0;
- // Get the exponent to scale the divisor.
- final int exponent = getMaxExponent(c, d);
+ // Get the exponent to scale the divisor parts to the range [1, 2).
+ final int exponent = getScale(c, d);
if (exponent <= Double.MAX_EXPONENT) {
ilogbw = exponent;
c = Math.scalb(c, -ilogbw);
d = Math.scalb(d, -ilogbw);
}
final double denom = c * c + d * d;
+ // Divisor is in the range [1, 2).
+ // Avoid overflow if a or b are very big. Here we do not need to
+ // handle sub-normal exponents so just get the maximum exponent.
+ if (getMaxExponent(a, b) > Double.MAX_EXPONENT - 2) {
+ ilogbw -= 2;
+ a /= 4;
+ b /= 4;
+ }
double x = Math.scalb((a * c + b * d) / denom, -ilogbw);
double y = Math.scalb((b * c - a * d) / denom, -ilogbw);
// Recover infinities and zeros that computed as NaN+iNaN
@@ -2834,7 +2846,7 @@ public final class Complex implements Serializable {
} else {
// Over/underflow
// scale so that abs(x) is near 1, with even exponent.
- final int scale = getMaxExponent(x, y) & MASK_INT_TO_EVEN;
+ final int scale = getScale(x, y) & MASK_INT_TO_EVEN;
final double sx = Math.scalb(x, -scale);
final double sy = Math.scalb(y, -scale);
final double st = Math.sqrt(2 * (Math.sqrt(sx * sx + sy * sy) + sx));
@@ -3233,6 +3245,70 @@ public final class Complex implements Serializable {
}
/**
+ * Returns a scale suitable for use with {@link Math#scalb(double, int)} to normalise
+ * the number to the interval {@code [1, 2)}.
+ *
+ * <p>The scale is typically the largest unbiased exponent used in the representation of the
+ * two numbers. In contrast to {@link Math#getExponent(double)} this handles
+ * sub-normal numbers by computing the number of leading zeros in the mantissa
+ * and shifting the unbiased exponent. The result is that for all finite, non-zero,
+ * numbers {@code a, b}, the magnitude of {@code scalb(x, -getScale(x))} is
+ * always in the range {@code [1, 2)}, where {@code x = max(|a|, |b|)}.
+ *
+ * <p>This method is a functional equivalent of the c function ilogb(double) adapted for
+ * two input arguments.
+ *
+ * <p>The result is to be used to scale a complex number using {@link Math#scalb(double, int)}.
+ * Hence the special case of both zero arguments is handled using the return value for NaN
+ * as zero cannot be scaled. This is different from {@link Math#getExponent(double)}
+ * or {@link #getMaxExponent(double, double)}.
+ *
+ * <p>Special cases:
+ *
+ * <ul>
+ * <li>If either argument is NaN or infinite, then the result is
+ * {@link Double#MAX_EXPONENT} + 1.
+ * <li>If both arguments are zero, then the result is
+ * {@link Double#MAX_EXPONENT} + 1.
+ * </ul>
+ *
+ * @param a the first value
+ * @param b the second value
+ * @return The maximum unbiased exponent of the values to be used for scaling
+ * @see Math#getExponent(double)
+ * @see Math#scalb(double, int)
+ * @see <a href="http://www.cplusplus.com/reference/cmath/ilogb/">ilogb</a>
+ */
+ private static int getScale(double a, double b) {
+ // Only interested in the exponent and mantissa so remove the sign bit
+ final long x = Double.doubleToRawLongBits(a) & UNSIGN_MASK;
+ final long y = Double.doubleToRawLongBits(b) & UNSIGN_MASK;
+ // Only interested in the maximum
+ final long max = Math.max(x, y);
+ // Get the unbiased exponent
+ int exp = (int) ((max >>> 52) - EXPONENT_OFFSET);
+
+ // Do not distinguish nan/inf
+ if (exp == Double.MAX_EXPONENT + 1) {
+ return exp;
+ }
+ // Handle sub-normal numbers
+ if (exp == Double.MIN_EXPONENT - 1) {
+ // Special case for zero, return as nan/inf to indicate scaling is not possible
+ if (max == 0) {
+ return Double.MAX_EXPONENT + 1;
+ }
+ // A sub-normal number has an exponent below -1022. The amount below
+ // is defined by the number of shifts of the most significant bit in
+ // the mantissa that is required to get a 1 at position 53 (i.e. as
+ // if it were a normal number with assumed leading bit)
+ final long mantissa = max & MANTISSA_MASK;
+ exp -= Long.numberOfLeadingZeros(mantissa << 12);
+ }
+ return exp;
+ }
+
+ /**
* Returns the largest unbiased exponent used in the representation of the
* two numbers. Special cases:
*
@@ -3243,16 +3319,20 @@ public final class Complex implements Serializable {
* {@link Double#MIN_EXPONENT} -1.
* </ul>
*
+ * <p>This is used by {@link #divide(double, double, double, double)} as
+ * a simple detection that a number may overflow if multiplied
+ * by a value in the interval [1, 2).
+ *
* @param a the first value
* @param b the second value
* @return The maximum unbiased exponent of the values.
* @see Math#getExponent(double)
+ * @see #divide(double, double, double, double)
*/
private static int getMaxExponent(double a, double b) {
// This could return:
// Math.getExponent(Math.max(Math.abs(a), Math.abs(b)))
// A speed test is required to determine performance.
-
return Math.max(Math.getExponent(a), Math.getExponent(b));
}
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 9b2317e..5e4e1f4 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
@@ -614,9 +614,48 @@ public class ComplexEdgeCaseTest {
assertComplex(a, -a, name, operation, 1.6388720948399111e-154, -6.7884304867749655e-155);
}
- // Note:
- // multiply is tested in CStandardTest
- // divide is tested in CStandardTest
+ // Note: inf/nan edge cases for
+ // multiply/divide are tested in CStandardTest
+
+ @Test
+ public void testDivide() {
+ final String name = "divide";
+ final BiFunction<Complex, Complex, Complex> operation = Complex::divide;
+
+ // Should be able to divide by a complex whose absolute (c*c+d*d)
+ // overflows or underflows including all sub-normal numbers.
+
+ // Worst case is using Double.MIN_VALUE
+ // Should normalise c and d to range [1, 2) resulting in:
+ // c = d = 1
+ // c * c + d * d = 2
+ // scaled x = (a * c + b * d) / denom = Double.MIN_VALUE
+ // scaled y = (b * c - a * d) / denom = 0
+ // The values are rescaled by 1023 + 51 (shift the last bit of the 52 bit mantissa)
+ double x = Math.scalb(Double.MIN_VALUE, 1023 + 51);
+ Assertions.assertEquals(1.0, x);
+ // In other words the result is (x+iy) / (x+iy) = (1+i0)
+ // The result is the same if imaginary is zero (i.e. a real only divide)
+
+ assertComplex(Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE, name, operation, 1.0, 0.0);
+ assertComplex(Double.MAX_VALUE, 0.0, Double.MAX_VALUE, 0.0, name, operation, 1.0, 0.0);
+
+ assertComplex(1.0, 1.0, 1.0, 1.0, name, operation, 1.0, 0.0);
+ assertComplex(1.0, 0.0, 1.0, 0.0, name, operation, 1.0, 0.0);
+ // Should work for all small values
+ x = Double.MIN_NORMAL;
+ while (x != 0) {
+ assertComplex(x, x, x, x, name, operation, 1.0, 0.0);
+ assertComplex(x, 0, x, 0, name, operation, 1.0, 0.0);
+ x /= 2;
+ }
+
+ // Some cases of not self-divide
+ assertComplex(1, 1, Double.MIN_VALUE, Double.MIN_VALUE, name, operation, inf, 0);
+ // As computed by GNU g++
+ assertComplex(Double.MIN_NORMAL, Double.MIN_NORMAL, Double.MIN_VALUE, Double.MIN_VALUE, name, operation, 4503599627370496L, 0);
+ assertComplex(Double.MIN_VALUE, Double.MIN_VALUE, Double.MIN_NORMAL, Double.MIN_NORMAL, name, operation, 2.2204460492503131e-16, 0);
+ }
@Test
public void testPow() {
[commons-numbers] 01/02: Modify hypotenuse scaling formula for
x^2+y^2 in log().
Posted by 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 62ee60b9911c7ebaa94db494dc93deb47db82819
Author: aherbert <ah...@apache.org>
AuthorDate: Thu Jan 16 16:42:49 2020 +0000
Modify hypotenuse scaling formula for x^2+y^2 in log().
The modification change uses the method from FastMath.hypot().
---
.../src/main/java/org/apache/commons/numbers/complex/Complex.java | 8 +++++---
1 file changed, 5 insertions(+), 3 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 7a4dd6c..8987031 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
@@ -94,8 +94,8 @@ public final class Complex implements Serializable {
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;
+ /** Half the number of bits of precision of the mantissa of a {@code double} rounded up: {@code 27}. */
+ private static final double HALF_PRECISION = 27;
/** The bit representation of {@code -0.0}. */
private static final long NEGATIVE_ZERO_LONG_BITS = Double.doubleToLongBits(-0.0);
/** Exponent offset in IEEE754 representation. */
@@ -2340,7 +2340,9 @@ public final class Complex implements Serializable {
// Do scaling
final int expx = Math.getExponent(x);
final int expy = Math.getExponent(y);
- if (2 * (expx - expy) > PRECISION_1) {
+ // Hull et al: 2 * (expx - expy) > precision + 1
+ // Modified to use the pre-computed half precision
+ if ((expx - expy) > HALF_PRECISION) {
// y can be ignored
re = log.apply(x);
} else {