Approximation of sinh/cosh using exp can be overflow/underflow safe.
Updated the sinh(), cosh() and tanh() methods for large absolute real to
avoid overflow of exp(x) to infinite.
sinh()/cosh() now use a common function to compute the result as only
the sign of the result components differs. The computation of e^x can be
made in stages to allow computation of results even when e^x / 2
overflows.
tanh() now computes e^x in stages.
Added edge cases where exp(x) overflows but the result can be computed.
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 3b7eff0..81c3ac5 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
@@ -147,17 +147,37 @@
*/
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.
+ * 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)
- * when x is large:
+ * when x is large using exp(x):
* <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.783. However it is set to an integer (709)
- * to provide headroom.
+ * <p>This value can be used to approximate e^x using a product:
+ *
+ * <pre>
+ * e^x = product_n (e^m) * e^(x-nm)
+ * n = (int) x/m
+ * e.g. e^2000 = e^m * e^m * e^(2000 - 2m) </pre>
+ *
+ * <p>The value should be below ln(max_value) ~ 709.783.
+ * The value m is set to an integer for less error when subtracting m and chosen as
+ * even (m=708) as it is used as a threshold in tanh with m/2.
+ *
+ * <p>The value is used to compute e^x multiplied by a small number avoiding
+ * overflow (sinh/cosh) or a small number divided by e^x without underflow due to
+ * infinite e^x (tanh). The following conditions are used:
+ * <pre>
+ * 0.5 * e^m * Double.MIN_VALUE * e^m * e^m = Infinity
+ * 2.0 / e^m / e^m = 0.0 </pre>
*/
- private static final double SAFE_EXP_MAX = 709;
+ private static final double SAFE_EXP = 708;
+ /**
+ * The value of Math.exp(SAFE_EXP): e^708.
+ * To be used in overflow/underflow safe products of e^m to approximate e^x where x > m.
+ */
+ private static final double EXP_M = Math.exp(SAFE_EXP);
/** Serializable version identifier. */
private static final long serialVersionUID = 20180201L;
@@ -2022,16 +2042,76 @@
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));
+ if (x > SAFE_EXP) {
+ // Approximate sinh/cosh(x) using exp^|x| / 2
+ return coshsinh(x, real, imaginary, false, constructor);
}
+ // No overflow of sinh/cosh
return constructor.create(Math.cosh(real) * Math.cos(imaginary),
Math.sinh(real) * Math.sin(imaginary));
}
/**
+ * Compute cosh or sinh when the absolute real component |x| is large. In this case
+ * cosh(x) and sinh(x) can be approximated by exp(|x|) / 2:
+ *
+ * <pre>
+ * cosh(x+iy) real = (e^|x| / 2) * cos(y)
+ * cosh(x+iy) imag = (e^|x| / 2) * sin(y) * sign(x)
+ * sinh(x+iy) real = (e^|x| / 2) * cos(y) * sign(x)
+ * sinh(x+iy) imag = (e^|x| / 2) * sin(y)
+ * </pre>
+ *
+ * @param x Absolute real component |x|.
+ * @param real Real part (x).
+ * @param imaginary Imaginary part (y).
+ * @param sinh Set to true to compute sinh, otherwise cosh.
+ * @param constructor Constructor.
+ * @return The hyperbolic sine/cosine of the complex number.
+ */
+ private static Complex coshsinh(double x, double real, double imaginary, boolean sinh,
+ ComplexConstructor constructor) {
+ // Always require the cos and sin.
+ double re = Math.cos(imaginary);
+ double im = Math.sin(imaginary);
+ // Compute the correct function
+ if (sinh) {
+ re = changeSign(re, real);
+ } else {
+ im = changeSign(im, real);
+ }
+ // Multiply by (e^|x| / 2).
+ // Overflow safe computation since sin/cos can be very small allowing a result
+ // when e^x overflows: e^x / 2 = (e^m / 2) * e^m * e^(x-2m)
+ if (x > SAFE_EXP * 3) {
+ // e^x > e^m * e^m * e^m
+ // y * (e^m / 2) * e^m * e^m will overflow when starting with Double.MIN_VALUE.
+ // Note: Do not multiply by +inf to safeguard against sin(y)=0.0 which
+ // will create 0 * inf = nan.
+ re *= Double.MAX_VALUE * Double.MAX_VALUE * Double.MAX_VALUE;
+ im *= Double.MAX_VALUE * Double.MAX_VALUE * Double.MAX_VALUE;
+ } else {
+ // Initial part of (e^x / 2) using (e^m / 2)
+ re *= EXP_M / 2;
+ im *= EXP_M / 2;
+ double xm;
+ if (x > SAFE_EXP * 2) {
+ // e^x = e^m * e^m * e^(x-2m)
+ re *= EXP_M;
+ im *= EXP_M;
+ xm = x - SAFE_EXP * 2;
+ } else {
+ // e^x = e^m * e^(x-m)
+ xm = x - SAFE_EXP;
+ }
+ final double exp = Math.exp(xm);
+ re *= exp;
+ im *= exp;
+ }
+ return constructor.create(re, im);
+ }
+
+ /**
* Returns the
* <a href="http://mathworld.wolfram.com/ExponentialFunction.html">
* exponential function</a> of this complex number.
@@ -2095,7 +2175,7 @@
return new Complex(zeroOrInf * Math.cos(imaginary),
zeroOrInf * Math.sin(imaginary));
} else if (Double.isNaN(real)) {
- // (NaN + i0) returns (NaN + i0);
+ // (NaN + i0) returns (NaN + i0)
// (NaN + iy) returns (NaN + iNaN) and optionally raises the invalid floating-point exception
// (NaN + iNaN) returns (NaN + iNaN)
return imaginary == 0 ? this : NAN;
@@ -2617,11 +2697,11 @@
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));
+ if (x > SAFE_EXP) {
+ // Approximate sinh/cosh(x) using exp^|x| / 2
+ return coshsinh(x, real, imaginary, true, constructor);
}
+ // No overflow of sinh/cosh
return constructor.create(Math.sinh(real) * Math.cos(imaginary),
Math.cosh(real) * Math.sin(imaginary));
}
@@ -2837,7 +2917,8 @@
// Cache the absolute real value
final double x = Math.abs(real);
- // Handle inf or nan
+ // Handle inf or nan.
+ // Deliberate logic inversion using x to match !Double.isFinite(x) knowing x is absolute.
if (!(x <= Double.MAX_VALUE) || !Double.isFinite(imaginary)) {
if (isPosInfinite(x)) {
if (Double.isFinite(imaginary)) {
@@ -2878,19 +2959,30 @@
// 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))
- // 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.
+ if (x > SAFE_EXP / 2) {
+ // Potential overflow in sinh/cosh(2x).
+ // Approximate sinh/cosh using exp^x.
+ // Ignore cos^2(y) in the divisor as it is insignificant.
// 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));
+ final double re = Math.copySign(1, real);
+ // imag = sin(2y) / 2 sinh^2(x)
+ // 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|
+ // Underflow safe divide as e^2|x| may overflow:
+ // imag = 2 sin(2y) / e^m / e^(2|x| - m)
+ double im = sin2(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);
+ }
+ return constructor.create(re, im);
}
// No overflow of sinh(2x) and cosh(2x)
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 606a9b1..400b843 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
@@ -247,11 +247,53 @@
assertComplex(small, big, name, operation, -0.99998768942655991, 1.1040715888508271e-310);
assertComplex(small, medium, name, operation, -0.41614683654714241, 2.0232539340376892e-308);
assertComplex(small, small, name, operation, 1, 0);
+
+ // Overflow test.
+ // Based on MATH-901 discussion of FastMath functionality.
+ // https://issues.apache.org/jira/browse/MATH-901#comment-13500669
+ // sinh(x)/cosh(x) can be approximated by exp(x) but must be overflow safe.
+
+ // sinh(x) = sign(x) * e^|x| / 2 when x is large.
+ // cosh(x) = e^|x| / 2 when x is large.
+ // Thus e^|x| can overflow but e^|x| / 2 may not.
+ // (e^|x| / 2) * sin/cos will always be smaller.
+ final double tiny = Double.MIN_VALUE;
+ final double x = 709.783;
+ Assertions.assertEquals(inf, Math.exp(x));
+ // As computed by GNU g++
+ assertComplex(x, 0, name, operation, 8.9910466927705402e+307, 0.0);
+ assertComplex(-x, 0, name, operation, 8.9910466927705402e+307, -0.0);
+ // sub-normal number x:
+ // cos(x) = 1 => real = (e^|x| / 2)
+ // sin(x) = x => imaginary = x * (e^|x| / 2)
+ assertComplex(x, small, name, operation, 8.9910466927705402e+307, 2.0005742956701358);
+ assertComplex(-x, small, name, operation, 8.9910466927705402e+307, -2.0005742956701358);
+ assertComplex(x, tiny, name, operation, 8.9910466927705402e+307, 4.4421672910524807e-16);
+ assertComplex(-x, tiny, name, operation, 8.9910466927705402e+307, -4.4421672910524807e-16);
+ // Should not overflow imaginary.
+ assertComplex(2 * x, tiny, name, operation, inf, 7.9879467061901743e+292);
+ assertComplex(-2 * x, tiny, name, operation, inf, -7.9879467061901743e+292);
+ // Test when large enough to overflow any non-zero value to infinity. Result should be
+ // as if x was infinite and y was finite.
+ assertComplex(3 * x, tiny, name, operation, inf, inf);
+ assertComplex(-3 * x, tiny, name, operation, inf, -inf);
+ // pi / 2 x:
+ // cos(x) = ~0 => real = x * (e^|x| / 2)
+ // sin(x) = ~1 => imaginary = (e^|x| / 2)
+ final double pi2 = Math.PI / 2;
+ assertComplex(x, pi2, name, operation, 5.5054282766429199e+291, 8.9910466927705402e+307);
+ assertComplex(-x, pi2, name, operation, 5.5054282766429199e+291, -8.9910466927705402e+307);
+ assertComplex(2 * x, pi2, name, operation, inf, inf);
+ assertComplex(-2 * x, pi2, name, operation, inf, -inf);
+ // Test when large enough to overflow any non-zero value to infinity. Result should be
+ // as if x was infinite and y was finite.
+ assertComplex(3 * x, pi2, name, operation, inf, inf);
+ assertComplex(-3 * x, pi2, name, operation, inf, -inf);
}
@Test
public void testSinh() {
- // sinh(a + b i) = sinh(a)cos(b)) + i cosh(a)sin(b)
+ // sinh(a + b i) = sinh(a)cos(b) + i cosh(a)sin(b)
// Odd function: negative real cases defined by positive real cases
final String name = "sinh";
final UnaryOperator<Complex> operation = Complex::sinh;
@@ -270,6 +312,46 @@
assertComplex(small, big, name, operation, -2.2250464665720564e-308, 0.004961954789184062);
assertComplex(small, medium, name, operation, -9.2595744730151568e-309, 0.90929742682568171);
assertComplex(small, small, name, operation, 2.2250738585072014e-308, 2.2250738585072014e-308);
+
+ // Overflow test.
+ // As per cosh with sign changes to real and imaginary
+
+ // sinh(x) = sign(x) * e^|x| / 2 when x is large.
+ // cosh(x) = e^|x| / 2 when x is large.
+ // Thus e^|x| can overflow but e^|x| / 2 may not.
+ // sinh(x) * sin/cos will always be smaller.
+ final double tiny = Double.MIN_VALUE;
+ final double x = 709.783;
+ Assertions.assertEquals(inf, Math.exp(x));
+ // As computed by GNU g++
+ assertComplex(x, 0, name, operation, 8.9910466927705402e+307, 0.0);
+ assertComplex(-x, 0, name, operation, -8.9910466927705402e+307, 0.0);
+ // sub-normal number:
+ // cos(x) = 1 => real = (e^|x| / 2)
+ // sin(x) = x => imaginary = x * (e^|x| / 2)
+ assertComplex(x, small, name, operation, 8.9910466927705402e+307, 2.0005742956701358);
+ assertComplex(-x, small, name, operation, -8.9910466927705402e+307, 2.0005742956701358);
+ assertComplex(x, tiny, name, operation, 8.9910466927705402e+307, 4.4421672910524807e-16);
+ assertComplex(-x, tiny, name, operation, -8.9910466927705402e+307, 4.4421672910524807e-16);
+ // Should not overflow imaginary.
+ assertComplex(2 * x, tiny, name, operation, inf, 7.9879467061901743e+292);
+ assertComplex(-2 * x, tiny, name, operation, -inf, 7.9879467061901743e+292);
+ // Test when large enough to overflow any non-zero value to infinity. Result should be
+ // as if x was infinite and y was finite.
+ assertComplex(3 * x, tiny, name, operation, inf, inf);
+ assertComplex(-3 * x, tiny, name, operation, -inf, inf);
+ // pi / 2 x:
+ // cos(x) = ~0 => real = x * (e^|x| / 2)
+ // sin(x) = ~1 => imaginary = (e^|x| / 2)
+ final double pi2 = Math.PI / 2;
+ assertComplex(x, pi2, name, operation, 5.5054282766429199e+291, 8.9910466927705402e+307);
+ assertComplex(-x, pi2, name, operation, -5.5054282766429199e+291, 8.9910466927705402e+307);
+ assertComplex(2 * x, pi2, name, operation, inf, inf);
+ assertComplex(-2 * x, pi2, name, operation, -inf, inf);
+ // Test when large enough to overflow any non-zero value to infinity. Result should be
+ // as if x was infinite and y was finite.
+ assertComplex(3 * x, pi2, name, operation, inf, inf);
+ assertComplex(-3 * x, pi2, name, operation, -inf, inf);
}
@Test
@@ -305,6 +387,18 @@
// cosh(2a) -> 0
assertComplex(Double.MIN_NORMAL, 1, name, operation, 7.6220323800193346e-308, 1.5574077246549021);
assertComplex(Double.MIN_VALUE, 1, name, operation, 1.4821969375237396e-323, 1.5574077246549021);
+
+ // Underflow test.
+ // sinh(x) can be approximated by exp(x) but must be overflow safe.
+ // im = 2 sin(2y) / e^2|x|
+ // This can be computed when e^2|x| only just overflows.
+ // Set a case where e^2|x| overflows but the imaginary can be computed
+ double x = 709.783 / 2;
+ double y = Math.PI / 4;
+ Assertions.assertEquals(1.0, Math.sin(2 * y), 1e-16);
+ Assertions.assertEquals(Double.POSITIVE_INFINITY, Math.exp(2 * x));
+ // As computed by GNU g++
+ assertComplex(x, y, name, operation, 1, 1.1122175583895849e-308);
}
@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 6a910f4..a98b65d 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
@@ -1967,16 +1967,16 @@
}
@Test
- public void testTanhAssumptions() {
- // Use the same constants used by tanh
- final double safeExpMax = 709;
+ public void testCoshSinhTanhAssumptions() {
+ // Use the same constants used to approximate cosh/sinh with e^|x| / 2
+ final double safeExpMax = 708;
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)));
+ Assertions.assertTrue(Double.isInfinite(Math.exp(safeExpMax + 2)));
// Can we assume cosh(x) = (e^x + e^-x) / 2 = e^|x| / 2
Assertions.assertEquals(big + small, big);
@@ -1992,6 +1992,15 @@
// Can we assume sinh(x/2) * cosh(x/2) is finite
Assertions.assertTrue(Double.isFinite(Math.sinh(safeExpMax / 2) * Math.cosh(safeExpMax / 2)));
+ // Will 2.0 / (2 * (e^|x|)) underflow
+ Assertions.assertNotEquals(0.0, 2.0 / big);
+ Assertions.assertEquals(0.0, 2.0 / big / big);
+
+ // This is an assumption used in sinh/cosh.
+ // Will 3 * (e^|x|/2) * y overflow for any positive y
+ 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;
@@ -2006,6 +2015,10 @@
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)));
+
+ // 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
+ Assertions.assertTrue(2.0 / Double.MAX_VALUE > 0);
}
/**
