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/14 02:42:28 UTC
[commons-numbers] 02/03: Better overflow handling of sinh/cosh.
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 642741cc893505a93351bbdfe7232dd8aa895ce6
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Tue Jan 14 01:00:17 2020 +0000
Better overflow handling of sinh/cosh.
At large x sinh and cosh can be approximated by exp(|x|) / 2.
---
.../apache/commons/numbers/complex/Complex.java | 109 ++++++++++++++++-----
.../commons/numbers/complex/ComplexTest.java | 59 +++++++++++
2 files changed, 141 insertions(+), 27 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 952c3b7..27737c3 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
@@ -146,6 +146,18 @@ 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;
+ /**
+ * A safe maximum double value {@code x} where {@code e^x} is not infinite.
+ * This can be used when functions require approximations of sinh(x) or cosh(x)
+ * when x is large:
+ * <pre>
+ * sinh(x) = (e^x - e^-x) / 2 = sign(x) * e^|x| / 2
+ * cosh(x) = (e^x + e^-x) / 2 = e^|x| / 2 </pre>
+ *
+ * <p>The value should be ln(max_value) ~ 709.08. However it is set to an integer (709)
+ * to provide headroom.
+ */
+ private static final double SAFE_EXP_MAX = 709;
/** Serializable version identifier. */
private static final long serialVersionUID = 20180201L;
@@ -2009,6 +2021,12 @@ public final class Complex implements Serializable {
// If real is NaN the sign of the imaginary part is unspecified.
return constructor.create(Math.cosh(real), changeSign(imaginary, real));
}
+ final double x = Math.abs(real);
+ if (x > SAFE_EXP_MAX) {
+ final double exp2 = Math.exp(x) / 2;
+ return constructor.create(exp2 * Math.cos(imaginary),
+ Math.copySign(exp2, real) * Math.sin(imaginary));
+ }
return constructor.create(Math.cosh(real) * Math.cos(imaginary),
Math.sinh(real) * Math.sin(imaginary));
}
@@ -2598,6 +2616,12 @@ public final class Complex implements Serializable {
// Real-only sinh(x).
return constructor.create(Math.sinh(real), imaginary);
}
+ final double x = Math.abs(real);
+ if (x > SAFE_EXP_MAX) {
+ final double exp2 = Math.exp(x) / 2;
+ return constructor.create(Math.copySign(exp2, real) * Math.cos(imaginary),
+ exp2 * Math.sin(imaginary));
+ }
return constructor.create(Math.sinh(real) * Math.cos(imaginary),
Math.cosh(real) * Math.sin(imaginary));
}
@@ -2810,46 +2834,77 @@ public final class Complex implements Serializable {
* @return The hyperbolic tangent of the complex number.
*/
private static Complex tanh(double real, double imaginary, ComplexConstructor constructor) {
- if (Double.isInfinite(real)) {
- if (Double.isFinite(imaginary)) {
- return constructor.create(Math.copySign(1, real), Math.copySign(0, sin2(imaginary)));
+ // Cache the absolute real value
+ final double x = Math.abs(real);
+
+ // Handle inf or nan
+ if (!(x <= Double.MAX_VALUE) || !Double.isFinite(imaginary)) {
+ if (isPosInfinite(x)) {
+ if (Double.isFinite(imaginary)) {
+ return constructor.create(Math.copySign(1, real),
+ Math.copySign(0, sin2(imaginary)));
+ }
+ // imaginary is infinite or NaN
+ return constructor.create(Math.copySign(1, real), Math.copySign(0, imaginary));
}
- // imaginary is infinite or NaN
- return constructor.create(Math.copySign(1, real), Math.copySign(0, imaginary));
+ // Remaining cases:
+ // (0 + i inf), returns (0 + i NaN)
+ // (0 + i NaN), returns (0 + i NaN)
+ // (x + i inf), returns (NaN + i NaN) for non-zero x (including infinite)
+ // (x + i NaN), returns (NaN + i NaN) for non-zero x (including infinite)
+ // (NaN + i 0), returns (NaN + i 0)
+ // (NaN + i y), returns (NaN + i NaN) for non-zero y (including infinite)
+ // (NaN + i NaN), returns (NaN + i NaN)
+ return constructor.create(real == 0 ? real : Double.NaN,
+ imaginary == 0 ? imaginary : Double.NaN);
}
+ // Finite components
+ // tanh(x+iy) = (sinh(2x) + i sin(2y)) / (cosh(2x) + cos(2y))
+
if (real == 0) {
// Imaginary-only tanh(iy) = i tan(y)
- if (Double.isFinite(imaginary)) {
- // Identity: sin x / (1 + cos x) = tan(x/2)
- return constructor.create(real, Math.tan(imaginary));
- }
- return constructor.create(real, Double.NaN);
+ // Identity: sin 2y / (1 + cos 2y) = tan(y)
+ return constructor.create(real, Math.tan(imaginary));
}
if (imaginary == 0) {
- if (Double.isNaN(real)) {
- return constructor.create(Double.NaN, imaginary);
- }
- // Identity: sinh x / (1 + cosh x) = tanh(x/2)
+ // Identity: sinh 2x / (1 + cosh 2x) = tanh(x)
return constructor.create(Math.tanh(real), imaginary);
}
- final double real2 = 2 * real;
+ // The double angles can be avoided using the identities:
+ // sinh(2x) = 2 sinh(x) cosh(x)
+ // sin(2y) = 2 sin(y) cos(y)
+ // 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))
+
+ // Check for overflow in sinh/cosh:
+ // sinh = (e^x - e^-x) / 2
+ // cosh = (e^x + e^-x) / 2
+ // sinh(x) -> sign(x) * e^|x| / 2 when x is large.
+ // cosh(x) -> e^|x| / 2 when x is large.
+ // sinh^2(x) -> e^2|x| / 4 when x is large.
+ if (x > SAFE_EXP_MAX / 2) {
+ // Ignore cos^2(y) in the divisor.
+ // real = sinh(x)cosh(x) / sinh^2(x) = +/-1
+ // imag = sin(2y) / (e^2|x| / 4) = 2 sin(2y) / e^2|x|
+ return constructor.create(Math.copySign(1, real),
+ 2 * sin2(imaginary) / Math.exp(2 * x));
+ }
- // Math.cosh returns positive infinity for infinity.
- // cosh -> inf
- final double divisor = Math.cosh(real2) + cos2(imaginary);
+ // No overflow of sinh(2x) and cosh(2x)
- // Math.sinh returns the input infinity for infinity.
- // sinh -> inf for positive x; else -inf
- final double sinhRe2 = Math.sinh(real2);
+ // 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.
- // Avoid inf / inf
- if (Double.isInfinite(sinhRe2) && Double.isInfinite(divisor)) {
- // Handle as if real was infinite
- return constructor.create(Math.copySign(1, real), Math.copySign(0, imaginary));
- }
- return constructor.create(sinhRe2 / divisor,
+ final double real2 = 2 * real;
+ final double divisor = Math.cosh(real2) + cos2(imaginary);
+ return constructor.create(Math.sinh(real2) / divisor,
sin2(imaginary) / divisor);
}
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 b791c9b..6a910f4 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
@@ -1965,4 +1965,63 @@ public class ComplexTest {
Assertions.assertNotEquals(0, 1 - Math.nextUp(1));
Assertions.assertNotEquals(0, 1 - Math.nextDown(1));
}
+
+ @Test
+ public void testTanhAssumptions() {
+ // Use the same constants used by tanh
+ final double safeExpMax = 709;
+
+ final double big = Math.exp(safeExpMax);
+ final double small = Math.exp(-safeExpMax);
+
+ // Overflow assumptions
+ Assertions.assertTrue(Double.isFinite(big));
+ Assertions.assertTrue(Double.isInfinite(Math.exp(safeExpMax + 1)));
+
+ // Can we assume cosh(x) = (e^x + e^-x) / 2 = e^|x| / 2
+ Assertions.assertEquals(big + small, big);
+ Assertions.assertEquals(Math.cosh(safeExpMax), big / 2);
+ Assertions.assertEquals(Math.cosh(-safeExpMax), big / 2);
+
+ // Can we assume sinh(x) = (e^x - e^-x) / 2 = sign(x) * e^|x| / 2
+ Assertions.assertEquals(big - small, big);
+ Assertions.assertEquals(small - big, -big);
+ Assertions.assertEquals(Math.sinh(safeExpMax), big / 2);
+ Assertions.assertEquals(Math.sinh(-safeExpMax), -big / 2);
+
+ // Can we assume sinh(x/2) * cosh(x/2) is finite
+ Assertions.assertTrue(Double.isFinite(Math.sinh(safeExpMax / 2) * Math.cosh(safeExpMax / 2)));
+
+ // 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)));
+ }
+
+ /**
+ * Test that sin and cos are linear around zero. This can be used for fast computation
+ * of sin and cos together when |x| is small.
+ */
+ @Test
+ public void testSinCosLinearAssumptions() {
+ // Are cos and sin linear around zero?
+ // If cos is still 1 then since d(sin) dx = cos then sin is linear.
+ Assertions.assertEquals(Math.cos(Double.MIN_NORMAL), 1.0);
+ Assertions.assertEquals(Math.sin(Double.MIN_NORMAL), Double.MIN_NORMAL);
+
+ // Are cosh and sinh linear around zero?
+ // If cosh is still 1 then since d(sinh) dx = cosh then sinh is linear.
+ Assertions.assertEquals(Math.cosh(Double.MIN_NORMAL), 1.0);
+ Assertions.assertEquals(Math.sinh(Double.MIN_NORMAL), Double.MIN_NORMAL);
+ }
}