From fdfbcf37db0cb9fdeadb2be6f77c83c92c901f37 Mon Sep 17 00:00:00 2001
From: diluculo <diluculo@gmail.com>
Date: Wed, 17 Jan 2018 18:45:16 +0900
Subject: [PATCH] Set huge = 22.0 Apply approximation to Sinh(z) and Cosh(z).

---
 src/Numerics/Trigonometry.cs      | 38 ++++++++++++++++++++++---------
 src/UnitTests/TrigonometryTest.cs |  2 ++
 2 files changed, 29 insertions(+), 11 deletions(-)

diff --git a/src/Numerics/Trigonometry.cs b/src/Numerics/Trigonometry.cs
index 45b4b8fa..fc75847d 100644
--- a/src/Numerics/Trigonometry.cs
+++ b/src/Numerics/Trigonometry.cs
@@ -447,9 +447,20 @@ namespace MathNet.Numerics
                 return new Complex(Sinh(value.Real), 0.0);
             }
 
+            // sinh(x + j y) = sinh(x)*cos(y) + j*cosh(x)*sin(y)
+            // if x > huge, sinh(x + jy) = sign(x)*exp(|x|)/2*cos(y) + j*exp(|x|)/2*sin(y)
+
+            if (Math.Abs(value.Real) >= 22.0) // Taken from the msun library in FreeBSD
+            {
+                double h = Math.Exp(Math.Abs(value.Real)) * 0.5;
+                return new Complex(
+                    Math.Sign(value.Real)*h*Cos(value.Imaginary),
+                    h*Sin(value.Imaginary));
+            }
+
             return new Complex(
                 Sinh(value.Real) * Cos(value.Imaginary),
-                Cosh(value.Real) * Sin(value.Imaginary));
+                Cosh(value.Real) * Sin(value.Imaginary));            
         }
 
         /// <summary>
@@ -474,6 +485,17 @@ namespace MathNet.Numerics
                 return new Complex(Cosh(value.Real), 0.0);
             }
 
+            // cosh(x + j*y) = cosh(x)*cos(y) + j*sinh(x)*sin(y)
+            // if x > huge, cosh(x + j*y) = exp(|x|)/2*cos(y) + j*sign(x)*exp(|x|)/2*sin(y)
+
+            if (Math.Abs(value.Real) >= 22.0) // Taken from the msun library in FreeBSD
+            {
+                double h = Math.Exp(Math.Abs(value.Real)) * 0.5;
+                return new Complex(
+                    h * Cos(value.Imaginary),
+                    Math.Sign(value.Real) * h * Sin(value.Imaginary));
+            }
+
             return new Complex(
                 Cosh(value.Real) * Cos(value.Imaginary),
                 Sinh(value.Real) * Sin(value.Imaginary));
@@ -519,10 +541,8 @@ namespace MathNet.Numerics
             // if tan(y) = +/- oo or 1/cos^2(y) = 1 + tan^2(y) = oo, tanh(z) = cosh(x)/sinh(x)
             //
             // The algorithm is based on Kahan.
-
-            double upperLimit = 177.6189650; // Asinh(double.MaxValue) / 4.0;
-
-            if (Math.Abs(value.Real) > upperLimit)
+            
+            if (Math.Abs(value.Real) >= 22.0) // Taken from the msun library in FreeBSD
             {
                 double e = Math.Exp(-Math.Abs(value.Real));
                 return e == 0.0
@@ -610,15 +630,13 @@ namespace MathNet.Numerics
             //
             // The algorithm is based on Kahan.
 
-            double upperLimit = 177.6189650; // Asinh(double.MaxValue) / 4.0;
-
             double tani = Tan(value.Imaginary);
             double cosi = Cos(value.Imaginary);
             double beta = 1.0 + tani * tani;
             double sinhr = Math.Sinh(value.Real);
             double coshr = Math.Cosh(value.Real);
             
-            if (Math.Abs(value.Real) > upperLimit)
+            if (Math.Abs(value.Real) >= 22.0) // Taken from the msun library in FreeBSD
             {
                 double e = Math.Exp(-Math.Abs(value.Real));
                 return e == 0.0
@@ -664,15 +682,13 @@ namespace MathNet.Numerics
             //
             // The algorithm is based on Kahan.
 
-            double upperLimit = 177.6189650; // Asinh(double.MaxValue) / 4.0;
-
             double coti = Cot(value.Imaginary);
             double sini = Sin(value.Imaginary);
             double beta = 1 + coti * coti;
             double sinhr = Sinh(value.Real);
             double coshr = Cosh(value.Real);
             
-            if (Math.Abs(value.Real) > upperLimit)
+            if (Math.Abs(value.Real) >= 22.0) // Taken from the msun library in FreeBSD
             {
                 double e = Math.Exp(-Math.Abs(value.Real));
                 return e == 0.0
diff --git a/src/UnitTests/TrigonometryTest.cs b/src/UnitTests/TrigonometryTest.cs
index a7c4bc82..c55e0d27 100644
--- a/src/UnitTests/TrigonometryTest.cs
+++ b/src/UnitTests/TrigonometryTest.cs
@@ -660,6 +660,7 @@ namespace MathNet.Numerics.UnitTests
         [TestCase(-512.0, -32.0, -9.52855589474962752845e221, -6.29843301304522728211e221)]
         [TestCase(-32.0, 512.0, +3.93564576751566e13, +3.13950783826274e12)]
         [TestCase(32.0, -512.0, -3.93564576751566e13, -3.13950783826274e12)]
+        [TestCase(709.0, 709.0, +2.22042071181169647963e307, -3.45764185010438140812e307)]
         public void CanComputeComplexHyperbolicSine(double real, double imag, double expectedReal, double expectedImag)
         {
             var actual = new Complex(real, imag).Sinh();
@@ -689,6 +690,7 @@ namespace MathNet.Numerics.UnitTests
         [TestCase(-512.0, -32.0, 9.52855589474962752845e221, 6.29843301304522728211e221)]
         [TestCase(-32.0, 512.0, -3.93564576751566e13, -3.13950783826274e12)]
         [TestCase(32.0, -512.0, -3.93564576751566e13, -3.13950783826274e12)]
+        [TestCase(709.0, 709.0, +2.22042071181169647963e307, -3.45764185010438140812e307)]
         public void CanComputeComplexHyperbolicCosine(double real, double imag, double expectedReal, double expectedImag)
         {
             var actual = new Complex(real, imag).Cosh();