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/02/03 16:53:44 UTC
[commons-numbers] branch master updated (79fbf9d -> 8af1a5d)
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 79fbf9d Added more random distribution types to the Complex performance test.
new b2622ef Complex: log() and sqrt() should use the same magnitude as abs().
new c79bdff sqrt() to use multiply by sqrt(2) / 2 for imaginary only numbers.
new 8af1a5d Change tanh to use the identities for all computations.
The 3 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 | 278 ++++++++++++---------
.../commons/numbers/complex/CReferenceTest.java | 2 +-
.../numbers/complex/ComplexEdgeCaseTest.java | 76 ++++--
.../commons/numbers/complex/ComplexTest.java | 143 +++++++++--
4 files changed, 333 insertions(+), 166 deletions(-)
[commons-numbers] 02/03: sqrt() to use multiply by sqrt(2) / 2 for
imaginary only numbers.
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 c79bdffb5c79fe19a662a42b8c69ad3869ed75f8
Author: aherbert <ah...@apache.org>
AuthorDate: Mon Feb 3 15:51:07 2020 +0000
sqrt() to use multiply by sqrt(2) / 2 for imaginary only numbers.
The change alters the divide by sqrt(2) to a multiply by the closest
double representation of 1 / sqrt(2). This is sqrt(2) / 2.
This adds tests to ensure the polar representation is not used when
either a real or imaginary only number:
real = sqrt(abs()) * Math.cos(arg() / 2)
imag = sqrt(abs()) * Math.sin(arg() / 2)
Due to floating-point errors in the trignomic functions and the
inexactness of Math.PI the result is not the correct result. It
explicitly checks conditions that would fail if using the polar
computation with the trigonomic functions cos and sin.
---
.../apache/commons/numbers/complex/Complex.java | 26 +++++---
.../numbers/complex/ComplexEdgeCaseTest.java | 75 ++++++++++++++--------
2 files changed, 67 insertions(+), 34 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 fb6174e..4498530 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
@@ -91,7 +91,12 @@ public final class Complex implements Serializable {
/** {@code 1/2}. */
private static final double HALF = 0.5;
/** {@code sqrt(2)}. */
- private static final double ROOT2 = Math.sqrt(2);
+ private static final double ROOT2 = 1.4142135623730951;
+ /** {@code 1.0 / sqrt(2)}.
+ * This is pre-computed to the closest double from the exact result.
+ * It is 1 ULP different from 1.0 / Math.sqrt(2) but equal to Math.sqrt(2) / 2.
+ */
+ private static final double ONE_OVER_ROOT2 = 0.7071067811865476;
/** The bit representation of {@code -0.0}. */
private static final long NEGATIVE_ZERO_LONG_BITS = Double.doubleToLongBits(-0.0);
/** Exponent offset in IEEE754 representation. */
@@ -2859,8 +2864,13 @@ public final class Complex implements Serializable {
}
return new Complex(sqrtAbs, imaginary);
} else if (x == 0) {
- // Imaginary only
- final double sqrtAbs = Math.sqrt(y) / ROOT2;
+ // Imaginary only. This sets the two components to the same magnitude.
+ // Note: In polar coordinates this does not happen:
+ // real = sqrt(abs()) * Math.cos(arg() / 2)
+ // imag = sqrt(abs()) * Math.sin(arg() / 2)
+ // arg() / 2 = pi/4 and cos and sin should both return sqrt(2)/2 but
+ // are different by 1 ULP.
+ final double sqrtAbs = Math.sqrt(y) * ONE_OVER_ROOT2;
return new Complex(sqrtAbs, Math.copySign(sqrtAbs, imaginary));
} else {
// Over/underflow.
@@ -2874,13 +2884,13 @@ public final class Complex implements Serializable {
double sx;
double sy;
if (Math.max(x, y) > SQRT_SAFE_UPPER) {
- // Overflow
- sx = x * 0x1.0p-4;
- sy = y * 0x1.0p-4;
- rescale = 0x1.0p2;
+ // Overflow. Scale down by 16 and rescale by sqrt(16).
+ sx = x / 16;
+ sy = y / 16;
+ rescale = 4;
} else {
// Sub-normal numbers. Make them normal by scaling by 2^54,
- // i.e. more than the mantissa digits and rescale by half.
+ // i.e. more than the mantissa digits, and rescale by sqrt(2^54) = 2^27.
sx = x * 0x1.0p54;
sy = y * 0x1.0p54;
rescale = 0x1.0p-27;
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 ebd8a87..ac20b48 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
@@ -562,41 +562,64 @@ public class ComplexEdgeCaseTest {
final String name = "sqrt";
final UnaryOperator<Complex> operation = Complex::sqrt;
- // Computed in polar coordinates:
- // real = (x^2 + y^2)^0.25 * cos(0.5 * atan2(y, x))
- // imag = (x^2 + y^2)^0.25 * sin(0.5 * atan2(y, x))
- // ---
- // Note:
- // If x is positive and y is +/-0.0 atan2 returns +/-0.
- // If x is negative and y is +/-0.0 atan2 returns +/-PI.
- // This causes problems as
- // cos(0.5 * PI) = 6.123233995736766e-17
- // assert: Math.cos(Math.acos(0)) != 0.0
- // Thus an unchecked polar computation will produce incorrect output when
- // there is no imaginary component and real is negative.
- // The computation should be done for real only numbers separately.
- // This condition is tested in the reference test against known data
- // so not repeated here.
- // ---
+ // Test real/imaginary only numbers satisfy the definition using polar coordinates:
+ // real = sqrt(abs()) * Math.cos(arg() / 2)
+ // imag = sqrt(abs()) * Math.sin(arg() / 2)
+ // However direct use of sin/cos will result in incorrect results due floating-point error.
+ // This test asserts the on the closest result to the exact answer which is possible
+ // if not using a simple polar computation.
+ // Note: If this test fails in the set-up assertions it is due to a change in the
+ // precision of java.util.Math.
+
+ // For positive real-only the argument is +/-0.
+ // For negative real-only the argument is +/-pi.
+ Assertions.assertEquals(0, Math.atan2(0, 1));
+ Assertions.assertEquals(Math.PI, Math.atan2(0, -1));
+ // In both cases the trigonomic functions should be exact but
+ // cos(pi/2) cannot be as pi/2 is not exact.
+ final double cosArgRe = 1.0;
+ final double sinArgRe = 0.0;
+ Assertions.assertNotEquals(0.0, Math.cos(Math.PI / 2), "Expected cos(pi/2) to be non-zero");
+ Assertions.assertEquals(0.0, Math.cos(Math.PI / 2), 6.123233995736766e-17);
+ // For imaginary-only the argument is Math.atan2(y, 0) = +/- pi/2.
+ Assertions.assertEquals(Math.PI / 2, Math.atan2(1, 0));
+ Assertions.assertEquals(-Math.PI / 2, Math.atan2(-1, 0));
+ // There is 1 ULP difference in the result of cos/sin of pi/4.
+ // It should be sqrt(2) / 2 for both.
+ final double cosArgIm = Math.cos(Math.PI / 4);
+ final double sinArgIm = Math.sin(Math.PI / 4);
+ final double root2over2 = Math.sqrt(2) / 2;
+ final double ulp = Math.ulp(cosArgIm);
+ Assertions.assertNotEquals(cosArgIm, sinArgIm, "Expected cos(pi/4) to not exactly equal sin(pi/4)");
+ Assertions.assertEquals(root2over2, cosArgIm, 0, "Expected cos(pi/4) to be sqrt(2) / 2");
+ Assertions.assertEquals(root2over2, sinArgIm, ulp, "Expected sin(pi/4) to be 1 ulp from sqrt(2) / 2");
+ for (final double a : new double[] {0.5, 1.0, 1.2322, 345345.234523}) {
+ final double rootA = Math.sqrt(a);
+ assertComplex(a, 0, name, operation, rootA * cosArgRe, rootA * sinArgRe, 0);
+ // This should be exact. It will fail if using the polar computation
+ // real = sqrt(abs()) * Math.cos(arg() / 2) as cos(pi/2) is not 0.0 but 6.123233995736766e-17
+ assertComplex(-a, 0, name, operation, rootA * sinArgRe, rootA * cosArgRe, 0);
+ // This should be exact. It won't be if Complex is using polar computation
+ // with sin/cos which does not output the same result for angle pi/4.
+ assertComplex(0, a, name, operation, rootA * root2over2, rootA * root2over2, 0);
+ assertComplex(0, -a, name, operation, rootA * root2over2, -rootA * root2over2, 0);
+ }
// Check overflow safe.
double a = Double.MAX_VALUE;
final double b = a / 4;
Assertions.assertEquals(inf, Complex.ofCartesian(a, b).abs(), "Expected overflow");
- // The expected absolute value has been compute using BigDecimal on Java 9
+ // The expected absolute value has been computed using BigDecimal on Java 9
//final double newAbs = new BigDecimal(a).multiply(new BigDecimal(a)).add(
// new BigDecimal(b).multiply(new BigDecimal(b)))
// .sqrt(MathContext.DECIMAL128).sqrt(MathContext.DECIMAL128).doubleValue()
final double newAbs = 1.3612566508088272E154;
-// assertComplex(a, b, name, operation, newAbs * Math.cos(0.5 * Math.atan2(b, a)),
-// newAbs * Math.sin(0.5 * Math.atan2(b, a)), 3);
-// assertComplex(b, a, name, operation, newAbs * Math.cos(0.5 * Math.atan2(a, b)),
-// newAbs * Math.sin(0.5 * Math.atan2(a, b)), 2);
+ assertComplex(a, b, name, operation, newAbs * Math.cos(0.5 * Math.atan2(b, a)),
+ newAbs * Math.sin(0.5 * Math.atan2(b, a)), 3);
+ assertComplex(b, a, name, operation, newAbs * Math.cos(0.5 * Math.atan2(a, b)),
+ newAbs * Math.sin(0.5 * Math.atan2(a, b)), 2);
- // In polar coords:
- // real = sqrt(abs()) * Math.cos(arg() / 2)
- // imag = sqrt(abs()) * Math.sin(arg() / 2)
- // This is possible if abs() does not overflow.
+ // Note that the computation is possible in polar coords if abs() does not overflow.
a = Double.MAX_VALUE / 2;
assertComplex(-a, a, name, operation, 4.3145940638864765e+153, 1.0416351505169177e+154, 2);
assertComplex(a, a, name, operation, 1.0416351505169177e+154, 4.3145940638864758e+153);
@@ -604,7 +627,7 @@ public class ComplexEdgeCaseTest {
assertComplex(a, -a, name, operation, 1.0416351505169177e+154, -4.3145940638864758e+153);
// Check minimum normal value conditions
- // Computing in Polar coords produces a very different result with
+ // Computing in polar coords produces a very different result with
// MIN_VALUE so use MIN_NORMAL
a = Double.MIN_NORMAL;
assertComplex(-a, a, name, operation, 6.7884304867749663e-155, 1.6388720948399111e-154);
[commons-numbers] 01/03: Complex: log() and sqrt() should use the
same magnitude as abs().
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 b2622effae96912f154fa7f629bb064a647886e8
Author: aherbert <ah...@apache.org>
AuthorDate: Sun Feb 2 14:29:18 2020 +0000
Complex: log() and sqrt() should use the same magnitude as abs().
Added tests to check the functions return what is documented as the
formula in the javadoc.
Fixed sqrt and log to use the hypot function to compute the absolute of
the complex number.
---
.../apache/commons/numbers/complex/Complex.java | 177 +++++++++++++--------
.../numbers/complex/ComplexEdgeCaseTest.java | 19 ++-
.../commons/numbers/complex/ComplexTest.java | 120 ++++++++++++++
3 files changed, 243 insertions(+), 73 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 04c7a6e..fb6174e 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
@@ -82,8 +82,6 @@ public final class Complex implements Serializable {
private static final double PI_OVER_2 = 0.5 * Math.PI;
/** π/4. */
private static final double PI_OVER_4 = 0.25 * Math.PI;
- /** Mask an integer number to even by discarding the lowest bit. */
- 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). */
@@ -94,8 +92,6 @@ public final class Complex implements Serializable {
private static final double HALF = 0.5;
/** {@code sqrt(2)}. */
private static final double ROOT2 = Math.sqrt(2);
- /** 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. */
@@ -150,6 +146,8 @@ public final class Complex implements Serializable {
* Equal to sqrt(u) * 2 with u the smallest normalised floating-point value.
*/
private static final double SAFE_LOWER = Math.sqrt(Double.MIN_NORMAL) * 2;
+ /** The safe maximum double value {@code x} to avoid overflow in sqrt. */
+ private static final double SQRT_SAFE_UPPER = Double.MAX_VALUE / 8;
/**
* A safe maximum double value {@code m} where {@code e^m} is not infinite.
* This can be used when functions require approximations of sinh(x) or cosh(x)
@@ -481,22 +479,57 @@ public final class Complex implements Serializable {
*
* <p>\[ \text{abs}(x + i y) = \sqrt{(x^2 + y^2)} \]
*
- * <p>If either component is infinite then the result is positive infinity. If either
- * component is NaN and this is not {@link #isInfinite() infinite} then the result is NaN.
+ * <p>Special cases:
+ *
+ * <ul>
+ * <li>{@code \text{abs}(x + iy) == \text{abs}(y + ix) == \text{abs}(x - iy)}.
+ * <li>If {@code z} is ±∞ + iy for any y, returns +∞.
+ * <li>If {@code z} is x + iNaN for non-infinite x, returns NaN.
+ * <li>If {@code z} is x + i0, returns |x|.
+ * </ul>
+ *
+ * <p>The cases ensure that if either component is infinite then the result is positive
+ * infinity. If either component is NaN and this is not {@link #isInfinite() infinite} then
+ * the result is NaN.
*
* <p>This code follows the
* <a href="http://www.iso-9899.info/wiki/The_Standard">ISO C Standard</a>, Annex G,
- * in calculating the returned value using the {@code hypot(x, y)} method for complex
- * \( x + i y \).
+ * in calculating the returned value without intermediate overflow or underflow.
*
* @return The absolute value.
* @see #isInfinite()
* @see #isNaN()
- * @see Math#hypot(double, double)
* @see <a href="http://mathworld.wolfram.com/ComplexModulus.html">Complex modulus</a>
*/
public double abs() {
- // Delegate
+ return abs(real, imaginary);
+ }
+
+ /**
+ * Returns the absolute value of the complex number.
+ * <pre>abs(x + i y) = sqrt(x^2 + y^2)</pre>
+ *
+ * <p>This should satisfy the special cases of the hypot function in ISO C99 F.9.4.3:
+ * "The hypot functions compute the square root of the sum of the squares of x and y,
+ * without undue overflow or underflow."
+ *
+ * <ul>
+ * <li>hypot(x, y), hypot(y, x), and hypot(x, −y) are equivalent.
+ * <li>hypot(x, ±0) is equivalent to |x|.
+ * <li>hypot(±∞, y) returns +∞, even if y is a NaN.
+ * </ul>
+ *
+ * <p>This method is called by all methods that require the absolute value of the complex
+ * number, e.g. abs(), sqrt() and log().
+ *
+ * @param real Real part.
+ * @param imaginary Imaginary part.
+ * @return The absolute value.
+ */
+ private static double abs(double real, double imaginary) {
+ // Delegate.
+ // Note that Math.hypot is slow in JDK 8. Using JDK 9+ improves performance 7-fold.
+ // This function is a candidate for optimisation.
return Math.hypot(real, imaginary);
}
@@ -2292,50 +2325,47 @@ public final class Complex implements Serializable {
double re;
+ // This alters the implementation of Hull et al (1994) which used a standard
+ // precision representation of |z|: sqrt(x*x + y*y).
+ // This formula should use the same definition of the magnitude returned
+ // by Complex.abs() which is a high precision computation with scaling.
+ // The checks for overflow thus only require ensuring the output of |z|
+ // will not overflow or underflow.
+
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 || y < SAFE_MIN) {
- // Over/underflow of sqrt(x^2+y^2)
- // Note: Since y<x no check for y > SAFE_MAX or x < SAFE_MIN.
- if (isPosInfinite(x)) {
- // Handle infinity
- re = x;
- } else {
- // Do scaling
- final int expx = Math.getExponent(x);
- final int expy = Math.getExponent(y);
- // 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 {
- // 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;
- }
- final double sx = Math.scalb(x, -scale);
- final double sy = Math.scalb(y, -scale);
- re = scale * logOf2 + 0.5 * log.apply(sx * sx + sy * sy);
+ } else {
+ // Check for over/underflow in |z|
+ // When scaling:
+ // log(a / b) = log(a) - log(b)
+ // So initialise the result with the log of the scale factor.
+ re = 0;
+ if (x > Double.MAX_VALUE / 2) {
+ // Potential overflow.
+ if (isPosInfinite(x)) {
+ // Handle infinity
+ return constructor.create(x, arg());
}
+ // Scale down.
+ x /= 2;
+ y /= 2;
+ // log(2)
+ re = logOf2;
+ } else if (y < Double.MIN_NORMAL) {
+ // Potential underflow.
+ if (y == 0) {
+ // Handle real only number
+ return constructor.create(log.apply(x), arg());
+ }
+ // Scale up sub-normal numbers to make them normal by scaling by 2^54,
+ // i.e. more than the mantissa digits.
+ x *= 0x1.0p54;
+ y *= 0x1.0p54;
+ // log(2^-54) = -54 * log(2)
+ re = -54 * logOf2;
}
- } else {
- // Safe region that avoids under/overflow
- re = 0.5 * log.apply(x * x + y * y);
+ re += log.apply(abs(x, y));
}
// All ISO C99 edge cases for the imaginary are satisfied by the Math library.
@@ -2794,15 +2824,22 @@ public final class Complex implements Serializable {
}
// Compute with positive values and determine sign at the end
- final double x = Math.abs(real);
- final double y = Math.abs(imaginary);
+ double x = Math.abs(real);
+ double y = Math.abs(imaginary);
// Compute
double t;
- if (inRegion(x, y, SAFE_MIN, SAFE_MAX)) {
- // No over/underflow of x^2 + y^2
- t = Math.sqrt(2 * (Math.sqrt(x * x + y * y) + x));
+ // This alters the implementation of Hull et al (1994) which used a standard
+ // precision representation of |z|: sqrt(x*x + y*y).
+ // This formula should use the same definition of the magnitude returned
+ // by Complex.abs() which is a high precision computation with scaling.
+ // Worry about overflow if 2 * (|z| + |x|) will overflow.
+ // Worry about underflow if |z| or |x| are sub-normal components.
+
+ if (inRegion(x, y, Double.MIN_NORMAL, SQRT_SAFE_UPPER)) {
+ // No over/underflow
+ t = Math.sqrt(2 * (abs(x, y) + x));
} else {
// Potential over/underflow. First check infinites and real/imaginary only.
@@ -2826,15 +2863,29 @@ public final class Complex implements Serializable {
final double sqrtAbs = Math.sqrt(y) / ROOT2;
return new Complex(sqrtAbs, Math.copySign(sqrtAbs, imaginary));
} else {
- // Over/underflow
- // scale so that abs(x) is near 1, with even exponent.
- 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));
- // Rescale. This works if exponent is even:
- // st * sqrt(2^scale) = st * (2^scale)^0.5 = st * 2^(scale*0.5)
- t = Math.scalb(st, scale / 2);
+ // Over/underflow.
+ // Full scaling is not required as this is done in the hypotenuse function.
+ // Keep the number as big as possible for maximum precision in the second sqrt.
+ // Note if we scale by an even power of 2, we can re-scale by sqrt of the number.
+ // a = sqrt(b)
+ // a = sqrt(b/4) * sqrt(4)
+
+ double rescale;
+ double sx;
+ double sy;
+ if (Math.max(x, y) > SQRT_SAFE_UPPER) {
+ // Overflow
+ sx = x * 0x1.0p-4;
+ sy = y * 0x1.0p-4;
+ rescale = 0x1.0p2;
+ } else {
+ // Sub-normal numbers. Make them normal by scaling by 2^54,
+ // i.e. more than the mantissa digits and rescale by half.
+ sx = x * 0x1.0p54;
+ sy = y * 0x1.0p54;
+ rescale = 0x1.0p-27;
+ }
+ t = rescale * Math.sqrt(2 * (abs(sx, sy) + sx));
}
}
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 5e4e1f4..ebd8a87 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
@@ -583,16 +583,15 @@ public class ComplexEdgeCaseTest {
double a = Double.MAX_VALUE;
final double b = a / 4;
Assertions.assertEquals(inf, Complex.ofCartesian(a, b).abs(), "Expected overflow");
- // Compute the expected new magnitude by expressing b as a scale factor of a:
- // (x^2 + y^2)^0.25
- // = sqrt(sqrt(a^2 + (b/a)^2 * a^2))
- // = sqrt(sqrt((1+(b/a)^2) * a^2))
- // = sqrt(sqrt((1+(b/a)^2))) * sqrt(a)
- final double newAbs = Math.sqrt(Math.sqrt(1 + (b / a) * (b / a))) * Math.sqrt(a);
- assertComplex(a, b, name, operation, newAbs * Math.cos(0.5 * Math.atan2(b, a)),
- newAbs * Math.sin(0.5 * Math.atan2(b, a)), 3);
- assertComplex(b, a, name, operation, newAbs * Math.cos(0.5 * Math.atan2(a, b)),
- newAbs * Math.sin(0.5 * Math.atan2(a, b)), 2);
+ // The expected absolute value has been compute using BigDecimal on Java 9
+ //final double newAbs = new BigDecimal(a).multiply(new BigDecimal(a)).add(
+ // new BigDecimal(b).multiply(new BigDecimal(b)))
+ // .sqrt(MathContext.DECIMAL128).sqrt(MathContext.DECIMAL128).doubleValue()
+ final double newAbs = 1.3612566508088272E154;
+// assertComplex(a, b, name, operation, newAbs * Math.cos(0.5 * Math.atan2(b, a)),
+// newAbs * Math.sin(0.5 * Math.atan2(b, a)), 3);
+// assertComplex(b, a, name, operation, newAbs * Math.cos(0.5 * Math.atan2(a, b)),
+// newAbs * Math.sin(0.5 * Math.atan2(a, b)), 2);
// In polar coords:
// real = sqrt(abs()) * Math.cos(arg() / 2)
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 a98b65d..08b6118 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
@@ -22,6 +22,8 @@ import java.util.Arrays;
import java.util.List;
import java.util.function.BiFunction;
import java.util.function.DoubleFunction;
+import java.util.function.Supplier;
+
import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.rng.simple.RandomSource;
import org.junit.jupiter.api.Assertions;
@@ -2037,4 +2039,122 @@ public class ComplexTest {
Assertions.assertEquals(Math.cosh(Double.MIN_NORMAL), 1.0);
Assertions.assertEquals(Math.sinh(Double.MIN_NORMAL), Double.MIN_NORMAL);
}
+
+ /**
+ * Test the abs and sqrt functions are consistent. The definition of sqrt uses abs
+ * and the result should be computed using the same representation of the
+ * complex number's magnitude (abs). If the sqrt function uses a simple representation
+ * {@code sqrt(x^2 + y^2)} then this may have a 1 ulp or more difference from the high
+ * accuracy result computed by abs. This will propagate to create differences in sqrt.
+ *
+ * <p>Note: This test is separated from the similar test for log to allow testing
+ * different numbers.
+ */
+ @Test
+ public void testAbsVsSqrt() {
+ final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_128_PP);
+ // Note: All methods implement scaling to ensure the magnitude can be computed.
+ // Try very large or small numbers that will over/underflow to test that the scaling
+ // is consistent. Note that:
+ // - sqrt will reduce the size of the real and imaginary
+ // components when |z|>1 and increase them when |z|<1.
+
+ // Each sample fails approximately 3% of the time if using a standard x^2+y^2 in sqrt()
+ // and high accuracy representation in abs().
+ // Use 1000 samples to ensure the behaviour is OK.
+ // Do not use data which will over/underflow so we can use a simple computation in the test
+ assertAbsVsSqrt(1000, () -> Complex.ofCartesian(createFixedExponentNumber(rng, 1000),
+ createFixedExponentNumber(rng, 1000)));
+ assertAbsVsSqrt(1000, () -> Complex.ofCartesian(createFixedExponentNumber(rng, -1000),
+ createFixedExponentNumber(rng, -1000)));
+ }
+
+ private static void assertAbsVsSqrt(int samples, Supplier<Complex> supplier) {
+ // Note: All methods implement scaling to ensure the magnitude can be computed.
+ // Try very large or small numbers that will over/underflow to test that the scaling
+ // is consistent.
+ for (int i = 0; i < samples; i++) {
+ final Complex z = supplier.get();
+ final double abs = z.abs();
+ final double x = Math.abs(z.getReal());
+ final double y = Math.abs(z.getImaginary());
+
+ // Target the formula provided in the documentation for sqrt:
+ // sqrt(x + iy)
+ // t = sqrt( 2 (|x| + |x + iy|) )
+ // if x >= 0: (t/2, y/t)
+ // else : (|y| / t, t/2 * sgn(y))
+ // Note this is not the definitional polar computation using absolute and argument:
+ // real = sqrt(|z|) * cos(0.5 * arg(z))
+ // imag = sqrt(|z|) * sin(0.5 * arg(z))
+ final Complex c = z.sqrt();
+ final double t = Math.sqrt(2 * (x + abs));
+ if (z.getReal() >= 0) {
+ Assertions.assertEquals(t / 2, c.getReal());
+ Assertions.assertEquals(z.getImaginary() / t, c.getImaginary());
+ } else {
+ Assertions.assertEquals(y / t, c.getReal());
+ Assertions.assertEquals(Math.copySign(t / 2, z.getImaginary()), c.getImaginary());
+ }
+ }
+ }
+
+ /**
+ * Test the abs and log functions are consistent. The definition of log uses abs
+ * and the result should be computed using the same representation of the
+ * complex number's magnitude (abs). If the log function uses a simple representation
+ * {@code sqrt(x^2 + y^2)} then this may have a 1 ulp or more difference from the high
+ * accuracy result computed by abs. This will propagate to create differences in log.
+ *
+ * <p>Note: This test is separated from the similar test for sqrt to allow testing
+ * different numbers.
+ */
+ @Test
+ public void testAbsVsLog() {
+ final UniformRandomProvider rng = RandomSource.create(RandomSource.XO_RO_SHI_RO_128_PP);
+ // Note: All methods implement scaling to ensure the magnitude can be computed.
+ // Try very large or small numbers that will over/underflow to test that the scaling
+ // is consistent. Note that:
+ // - log will set the real component using log(|z|). This will massively reduce
+ // the magnitude when |z| >> 1. Highest accuracy will be when |z| is as large
+ // as possible before computing the log.
+
+ // No test around |z| == 1 as a high accuracy computation is required: Math.log1p(x*x+y*y-1)
+
+ // Each sample fails approximately 25% of the time if using a standard x^2+y^2 in log()
+ // and high accuracy representation in abs(). Use 100 samples to ensure the behaviour is OK.
+ assertAbsVsLog(100, () -> Complex.ofCartesian(createFixedExponentNumber(rng, 1022),
+ createFixedExponentNumber(rng, 1022)));
+ assertAbsVsLog(100, () -> Complex.ofCartesian(createFixedExponentNumber(rng, -1022),
+ createFixedExponentNumber(rng, -1022)));
+ }
+
+ private static void assertAbsVsLog(int samples, Supplier<Complex> supplier) {
+ // Note: All methods implement scaling to ensure the magnitude can be computed.
+ // Try very large or small numbers that will over/underflow to test that the scaling
+ // is consistent.
+ for (int i = 0; i < samples; i++) {
+ final Complex z = supplier.get();
+ final double abs = z.abs();
+ final double x = Math.abs(z.getReal());
+ final double y = Math.abs(z.getImaginary());
+
+ // log(x + iy) = log(|x + i y|) + i arg(x + i y)
+ // Only test the real component
+ final Complex c = z.log();
+ Assertions.assertEquals(Math.log(abs), c.getReal());
+ }
+ }
+
+ /**
+ * Creates a number in the range {@code [1, 2)} with up to 52-bits in the mantissa.
+ * Then modifies the exponent by the given amount.
+ *
+ * @param rng Source of randomness
+ * @param exponent Amount to change the exponent (in range [-1023, 1023])
+ * @return the number
+ */
+ private static double createFixedExponentNumber(UniformRandomProvider rng, int exponent) {
+ return Double.longBitsToDouble((rng.nextLong() >>> 12) | ((1023L + exponent) << 52));
+ }
}
[commons-numbers] 03/03: Change tanh to use the identities for all
computations.
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 8af1a5dd8d1efaf2a29ed403e98cbf216212220b
Author: aherbert <ah...@apache.org>
AuthorDate: Mon Feb 3 16:45:05 2020 +0000
Change tanh to use the identities for all computations.
This removes the requirement to assert that the switch from using
sin(2x) and cos(2x) and suitable identities is smooth. The test to
assert this failed on Java 11 but passed on Java 8. This may be due to a
change in the trigonomic functions between JDKs. The test assertion has
been removed as it is no longer valid for the new computation. The
functions to compute the identities for overflow are now redundant and
have been removed. The function is now overflow safe for finite
imaginary with no computation switch as y approaches infinity.
This lowers the ulp from the c reference test from 3 to 2 for tanh. Thus
the computation accuracy is similar and the c reference data may be
using the same computation.
---
.../apache/commons/numbers/complex/Complex.java | 85 +++++++---------------
.../commons/numbers/complex/CReferenceTest.java | 2 +-
.../commons/numbers/complex/ComplexTest.java | 23 ++----
3 files changed, 37 insertions(+), 73 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 4498530..5a322c3 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
@@ -2965,10 +2965,13 @@ public final class Complex implements Serializable {
* <a href="http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1892.htm#dr_471">
* DR 471: Complex math functions cacosh and ctanh</a>.
*
- * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
+ * <p>This is defined using real \( x \) and imaginary \( y \) parts:
*
* <p>\[ \tan(x + iy) = \frac{\sinh(2x)}{\cosh(2x)+\cos(2y)} + i \frac{\sin(2y)}{\cosh(2x)+\cos(2y)} \]
*
+ * <p>The implementation uses double-angle identities to avoid overflow of {@code 2x}
+ * and {@code 2y}.
+ *
* @return The hyperbolic tangent of this complex number.
* @see <a href="http://functions.wolfram.com/ElementaryFunctions/Tanh/">Tanh</a>
*/
@@ -2996,8 +2999,15 @@ public final class Complex implements Serializable {
if (!(x <= Double.MAX_VALUE) || !Double.isFinite(imaginary)) {
if (isPosInfinite(x)) {
if (Double.isFinite(imaginary)) {
+ // The sign is copied from sin(2y)
+ // The identity sin(2a) = 2 sin(a) cos(a) is used for consistency
+ // with the computation below. Only the magnitude is important
+ // so drop the 2. When |y| is small sign(sin(2y)) = sign(y).
+ final double sign = Math.abs(imaginary) < PI_OVER_2 ?
+ imaginary :
+ Math.sin(imaginary) * Math.cos(imaginary);
return constructor.create(Math.copySign(1, real),
- Math.copySign(0, sin2(imaginary)));
+ Math.copySign(0, sign));
}
// imaginary is infinite or NaN
return constructor.create(Math.copySign(1, real), Math.copySign(0, imaginary));
@@ -3033,7 +3043,8 @@ public final class Complex implements Serializable {
// cosh(2x) = 2 sinh^2(x) + 1
// cos(2y) = 2 cos^2(y) - 1
// tanh(x+iy) = (sinh(x)cosh(x) + i sin(y)cos(y)) / (sinh^2(x) + cos^2(y))
- // tanh(x+iy) = (sinh(x)cosh(x) + i 0.5 sin(2y)) / (sinh^2(x) + cos^2(y))
+ // To avoid a junction when swapping between the double angles and the identities
+ // the identities are used in all cases.
if (x > SAFE_EXP / 2) {
// Potential overflow in sinh/cosh(2x).
@@ -3045,74 +3056,34 @@ public final class Complex implements Serializable {
// sinh(x) -> sign(x) * e^|x| / 2 when x is large.
// sinh^2(x) -> e^2|x| / 4 when x is large.
// imag = sin(2y) / 2 (e^2|x| / 4) = 2 sin(2y) / e^2|x|
+ // = 4 * sin(y) cos(y) / e^2|x|
// Underflow safe divide as e^2|x| may overflow:
- // imag = 2 sin(2y) / e^m / e^(2|x| - m)
- double im = sin2(imaginary);
+ // imag = 4 * sin(y) cos(y) / e^m / e^(2|x| - m)
+ // (|im| is a maximum of 2)
+ double im = Math.sin(imaginary) * Math.cos(imaginary);
if (x > SAFE_EXP) {
// e^2|x| > e^m * e^m
// This will underflow 2.0 / e^m / e^m
im = Math.copySign(0.0, im);
} else {
// e^2|x| = e^m * e^(2|x| - m)
- im = 2 * im / EXP_M / Math.exp(2 * x - SAFE_EXP);
+ im = 4 * im / EXP_M / Math.exp(2 * x - SAFE_EXP);
}
return constructor.create(re, im);
}
// No overflow of sinh(2x) and cosh(2x)
- // Note: This does not use the double angle identities and returns the
- // definitional formula when 2y is finite. The transition when 2y overflows
- // and cos(2y) and sin(2y) switch to identities is monotonic at the junction
- // but the function is not smooth due to the sampling
- // of the 2 pi period at very large jumps of x. Thus this returns a value
- // but the usefulness as y -> inf may be limited.
-
- final double real2 = 2 * real;
- final double divisor = Math.cosh(real2) + cos2(imaginary);
- return constructor.create(Math.sinh(real2) / divisor,
- sin2(imaginary) / divisor);
- }
-
- /**
- * Safely compute {@code cos(2*a)} when {@code a} is finite.
- * Note that {@link Math#cos(double)} returns NaN when the input is infinite.
- * If {@code 2*a} is finite use {@code Math.cos(2*a)}; otherwise use the identity:
- *
- * <pre>
- * <code>cos(2a) = 2 cos<sup>2</sup>(a) - 1</code></pre>
- *
- * @param a Angle a.
- * @return The cosine of 2a.
- * @see Math#cos(double)
- */
- private static double cos2(double a) {
- final double twoA = 2 * a;
- if (Double.isFinite(twoA)) {
- return Math.cos(twoA);
- }
- final double cosA = Math.cos(a);
- return 2 * cosA * cosA - 1;
- }
+ // Note: This does not use the definitional formula but uses the identity:
+ // tanh(x+iy) = (sinh(x)cosh(x) + i sin(y)cos(y)) / (sinh^2(x) + cos^2(y))
- /**
- * Safely compute {@code sin(2*a)} when {@code a} is finite.
- * Note that {@link Math#sin(double)} returns NaN when the input is infinite.
- * If {@code 2*a} is finite use {@code Math.sin(2*a)}; otherwise use the identity:
- *
- * <pre>
- * <code>sin(2a) = 2 sin(a) cos(a)</code></pre>
- *
- * @param a Angle a.
- * @return The sine of 2a.
- * @see Math#sin(double)
- */
- private static double sin2(double a) {
- final double twoA = 2 * a;
- if (Double.isFinite(twoA)) {
- return Math.sin(twoA);
- }
- return 2 * Math.sin(a) * Math.cos(a);
+ final double sinhx = Math.sinh(real);
+ final double coshx = Math.cosh(real);
+ final double siny = Math.sin(imaginary);
+ final double cosy = Math.cos(imaginary);
+ final double divisor = sinhx * sinhx + cosy * cosy;
+ return constructor.create(sinhx * coshx / divisor,
+ siny * cosy / divisor);
}
/**
diff --git a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java
index 3b02907..6eab130 100644
--- a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java
+++ b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java
@@ -302,7 +302,7 @@ public class CReferenceTest {
@Test
public void testTanh() {
// Odd function: negative real cases defined by positive real cases
- assertOperation("tanh", Complex::tanh, 3);
+ assertOperation("tanh", Complex::tanh, 2);
}
@Test
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 08b6118..b66b872 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
@@ -1992,9 +1992,11 @@ public class ComplexTest {
Assertions.assertEquals(Math.sinh(-safeExpMax), -big / 2);
// Can we assume sinh(x/2) * cosh(x/2) is finite
+ // Can we assume sinh(x/2)^2 is finite
Assertions.assertTrue(Double.isFinite(Math.sinh(safeExpMax / 2) * Math.cosh(safeExpMax / 2)));
+ Assertions.assertTrue(Double.isFinite(Math.sinh(safeExpMax / 2) * Math.sinh(safeExpMax / 2)));
- // Will 2.0 / (2 * (e^|x|)) underflow
+ // Will 2.0 / e^|x| / e^|x| underflow
Assertions.assertNotEquals(0.0, 2.0 / big);
Assertions.assertEquals(0.0, 2.0 / big / big);
@@ -2003,20 +2005,11 @@ public class ComplexTest {
Assertions.assertTrue(Double.isFinite(0.5 * big * Double.MIN_VALUE * big));
Assertions.assertTrue(Double.isInfinite(0.5 * big * Double.MIN_VALUE * big * big));
- // Is overflow point monotonic:
- // The values at the switch are exact. Previous values are close.
- final double x = Double.MAX_VALUE / 2;
- final double x1 = Math.nextDown(x);
- final double x2 = Math.nextDown(x1);
- // sin(2x) = 2 sin(x) cos(x) (within 1 ulp)
- Assertions.assertEquals(Double.POSITIVE_INFINITY, 2 * Math.nextUp(x));
- Assertions.assertEquals(Math.sin(2 * x), 2 * Math.sin(x) * Math.cos(x));
- Assertions.assertEquals(Math.sin(2 * x1), 2 * Math.sin(x1) * Math.cos(x1), Math.ulp(Math.sin(2 * x1)));
- Assertions.assertEquals(Math.sin(2 * x2), 2 * Math.sin(x2) * Math.cos(x2), Math.ulp(Math.sin(2 * x2)));
- // cos(2x) = 2 * cos^2(x) - 1 (within 4 ulp)
- Assertions.assertEquals(Math.cos(2 * x), 2 * Math.cos(x) * Math.cos(x) - 1);
- Assertions.assertEquals(Math.cos(2 * x1), 2 * Math.cos(x1) * Math.cos(x1) - 1, 4 * Math.ulp(Math.cos(2 * x1)));
- Assertions.assertEquals(Math.cos(2 * x2), 2 * Math.cos(x2) * Math.cos(x2) - 1, 2 * Math.ulp(Math.cos(2 * x2)));
+ // Assume the sign of sin(2y) = sin(y) * cos(y) when |y| < pi/2
+ for (final double y : new double[] {Math.PI / 2, Math.PI / 4, 1.0, 0.5, 0.0}) {
+ Assertions.assertEquals(Math.signum(Math.sin(2 * y)), Math.signum(Math.sin(y) * Math.cos(y)));
+ Assertions.assertEquals(Math.signum(Math.sin(2 * -y)), Math.signum(Math.sin(-y) * Math.cos(-y)));
+ }
// tanh: 2.0 / Double.MAX_VALUE does not underflow.
// Thus 2 sin(2y) / e^2|x| can be computed when e^2|x| only just overflows