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);
     }
 
     /**