Fix catan, catanh missing underflows (bug 15406).

7 files changed, 128 insertions, 6 deletions

diff --git a/math/libm-test.inc b/math/libm-test.inc
index a3d7731..a62ae9a 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -4336,6 +4336,34 @@ catan_test (void)
   TEST_c_c (catan, 0.0L, -0x1p5000L, 1.570796326794896619231321691639751442099L, -7.079811261048172892385615158694057552948e-1506L);
   TEST_c_c (catan, -0.0L, -0x1p5000L, -1.570796326794896619231321691639751442099L, -7.079811261048172892385615158694057552948e-1506L);
 #endif
+  TEST_c_c (catan, 0x1p63L, 0.5L, 1.570796326794896619122901474391200998698L, 5.877471754111437539843682686111228389007e-39L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catan, 0x1p63L, -0.5L, 1.570796326794896619122901474391200998698L, -5.877471754111437539843682686111228389007e-39L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catan, -0x1p63L, 0.5L, -1.570796326794896619122901474391200998698L, 5.877471754111437539843682686111228389007e-39L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catan, -0x1p63L, -0.5L, -1.570796326794896619122901474391200998698L, -5.877471754111437539843682686111228389007e-39L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catan, 0.5L, 0x1p63L, 1.570796326794896619231321691639751442093L, 1.084202172485504434007452800869941711427e-19L);
+  TEST_c_c (catan, -0.5L, 0x1p63L, -1.570796326794896619231321691639751442093L, 1.084202172485504434007452800869941711427e-19L);
+  TEST_c_c (catan, 0.5L, -0x1p63L, 1.570796326794896619231321691639751442093L, -1.084202172485504434007452800869941711427e-19L);
+  TEST_c_c (catan, -0.5L, -0x1p63L, -1.570796326794896619231321691639751442093L, -1.084202172485504434007452800869941711427e-19L);
+#ifndef TEST_FLOAT
+  TEST_c_c (catan, 0x1p511L, 0.5L, 1.570796326794896619231321691639751442099L, 1.112536929253600691545116358666202032110e-308L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catan, 0x1p511L, -0.5L, 1.570796326794896619231321691639751442099L, -1.112536929253600691545116358666202032110e-308L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catan, -0x1p511L, 0.5L, -1.570796326794896619231321691639751442099L, 1.112536929253600691545116358666202032110e-308L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catan, -0x1p511L, -0.5L, -1.570796326794896619231321691639751442099L, -1.112536929253600691545116358666202032110e-308L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catan, 0.5L, 0x1p511L, 1.570796326794896619231321691639751442099L, 1.491668146240041348658193063092586767475e-154L);
+  TEST_c_c (catan, -0.5L, 0x1p511L, -1.570796326794896619231321691639751442099L, 1.491668146240041348658193063092586767475e-154L);
+  TEST_c_c (catan, 0.5L, -0x1p511L, 1.570796326794896619231321691639751442099L, -1.491668146240041348658193063092586767475e-154L);
+  TEST_c_c (catan, -0.5L, -0x1p511L, -1.570796326794896619231321691639751442099L, -1.491668146240041348658193063092586767475e-154L);
+#endif
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+  TEST_c_c (catan, 0x1p8191L, 0.5L, 1.570796326794896619231321691639751442099L, 1.681051571556046753131338908660876301299e-4932L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catan, 0x1p8191L, -0.5L, 1.570796326794896619231321691639751442099L, -1.681051571556046753131338908660876301299e-4932L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catan, -0x1p8191L, 0.5L, -1.570796326794896619231321691639751442099L, 1.681051571556046753131338908660876301299e-4932L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catan, -0x1p8191L, -0.5L, -1.570796326794896619231321691639751442099L, -1.681051571556046753131338908660876301299e-4932L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catan, 0.5L, 0x1p8191L, 1.570796326794896619231321691639751442099L, 1.833603867554847165621412392048483165956e-2466L);
+  TEST_c_c (catan, -0.5L, 0x1p8191L, -1.570796326794896619231321691639751442099L, 1.833603867554847165621412392048483165956e-2466L);
+  TEST_c_c (catan, 0.5L, -0x1p8191L, 1.570796326794896619231321691639751442099L, -1.833603867554847165621412392048483165956e-2466L);
+  TEST_c_c (catan, -0.5L, -0x1p8191L, -1.570796326794896619231321691639751442099L, -1.833603867554847165621412392048483165956e-2466L);
+#endif
 
   TEST_c_c (catan, 0.75L, 1.25L, 1.10714871779409050301706546017853704L, 0.549306144334054845697622618461262852L);
   TEST_c_c (catan, -2, -3, -1.4099210495965755225306193844604208L, -0.22907268296853876629588180294200276L);
@@ -4440,6 +4468,34 @@ catanh_test (void)
   TEST_c_c (catanh, 0.0L, -0x1p5000L, 0.0L, -1.570796326794896619231321691639751442099L);
   TEST_c_c (catanh, -0.0L, -0x1p5000L, -0.0L, -1.570796326794896619231321691639751442099L);
 #endif
+  TEST_c_c (catanh, 0x1p63L, 0.5L, 1.084202172485504434007452800869941711427e-19L, 1.570796326794896619231321691639751442093L);
+  TEST_c_c (catanh, 0x1p63L, -0.5L, 1.084202172485504434007452800869941711427e-19L, -1.570796326794896619231321691639751442093L);
+  TEST_c_c (catanh, -0x1p63L, 0.5L, -1.084202172485504434007452800869941711427e-19L, 1.570796326794896619231321691639751442093L);
+  TEST_c_c (catanh, -0x1p63L, -0.5L, -1.084202172485504434007452800869941711427e-19L, -1.570796326794896619231321691639751442093L);
+  TEST_c_c (catanh, 0.5L, 0x1p63L, 5.877471754111437539843682686111228389007e-39L, 1.570796326794896619122901474391200998698L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catanh, -0.5L, 0x1p63L, -5.877471754111437539843682686111228389007e-39L, 1.570796326794896619122901474391200998698L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catanh, 0.5L, -0x1p63L, 5.877471754111437539843682686111228389007e-39L, -1.570796326794896619122901474391200998698L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catanh, -0.5L, -0x1p63L, -5.877471754111437539843682686111228389007e-39L, -1.570796326794896619122901474391200998698L, UNDERFLOW_EXCEPTION_FLOAT);
+#ifndef TEST_FLOAT
+  TEST_c_c (catanh, 0x1p511L, 0.5L, 1.491668146240041348658193063092586767475e-154L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0x1p511L, -0.5L, 1.491668146240041348658193063092586767475e-154L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0x1p511L, 0.5L, -1.491668146240041348658193063092586767475e-154L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0x1p511L, -0.5L, -1.491668146240041348658193063092586767475e-154L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0.5L, 0x1p511L, 1.112536929253600691545116358666202032110e-308L, 1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catanh, -0.5L, 0x1p511L, -1.112536929253600691545116358666202032110e-308L, 1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catanh, 0.5L, -0x1p511L, 1.112536929253600691545116358666202032110e-308L, -1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catanh, -0.5L, -0x1p511L, -1.112536929253600691545116358666202032110e-308L, -1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION_DOUBLE);
+#endif
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+  TEST_c_c (catanh, 0x1p8191L, 0.5L, 1.833603867554847165621412392048483165956e-2466L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0x1p8191L, -0.5L, 1.833603867554847165621412392048483165956e-2466L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0x1p8191L, 0.5L, -1.833603867554847165621412392048483165956e-2466L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0x1p8191L, -0.5L, -1.833603867554847165621412392048483165956e-2466L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0.5L, 0x1p8191L, 1.681051571556046753131338908660876301299e-4932L, 1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catanh, -0.5L, 0x1p8191L, -1.681051571556046753131338908660876301299e-4932L, 1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catanh, 0.5L, -0x1p8191L, 1.681051571556046753131338908660876301299e-4932L, -1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catanh, -0.5L, -0x1p8191L, -1.681051571556046753131338908660876301299e-4932L, -1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION);
+#endif
 
   TEST_c_c (catanh, 0.75L, 1.25L, 0.261492138795671927078652057366532140L, 0.996825126463918666098902241310446708L);
   TEST_c_c (catanh, -2, -3, -0.14694666622552975204743278515471595L, -1.3389725222944935611241935759091443L);
diff --git a/math/s_catan.c b/math/s_catan.c
index 783941a..6a1016e 100644
--- a/math/s_catan.c
+++ b/math/s_catan.c
@@ -20,7 +20,7 @@
 #include <complex.h>
 #include <math.h>
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ double
 __catan (__complex__ double x)
@@ -83,6 +83,17 @@ __catan (__complex__ double x)
 	  num = 4.0 * __imag__ x;
 	  __imag__ res = 0.25 * __log1p (num / den);
 	}
+
+      if (fabs (__real__ res) < DBL_MIN)
+	{
+	  volatile double force_underflow = __real__ res * __real__ res;
+	  (void) force_underflow;
+	}
+      if (fabs (__imag__ res) < DBL_MIN)
+	{
+	  volatile double force_underflow = __imag__ res * __imag__ res;
+	  (void) force_underflow;
+	}
     }
 
   return res;
diff --git a/math/s_catanf.c b/math/s_catanf.c
index 3ffc6db..d0d188c 100644
--- a/math/s_catanf.c
+++ b/math/s_catanf.c
@@ -20,7 +20,7 @@
 #include <complex.h>
 #include <math.h>
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ float
 __catanf (__complex__ float x)
@@ -83,6 +83,17 @@ __catanf (__complex__ float x)
 	  num = 4.0f * __imag__ x;
 	  __imag__ res = 0.25f * __log1pf (num / den);
 	}
+
+      if (fabsf (__real__ res) < FLT_MIN)
+	{
+	  volatile float force_underflow = __real__ res * __real__ res;
+	  (void) force_underflow;
+	}
+      if (fabsf (__imag__ res) < FLT_MIN)
+	{
+	  volatile float force_underflow = __imag__ res * __imag__ res;
+	  (void) force_underflow;
+	}
     }
 
   return res;
diff --git a/math/s_catanh.c b/math/s_catanh.c
index 0ee8c64..5d18cd6 100644
--- a/math/s_catanh.c
+++ b/math/s_catanh.c
@@ -20,7 +20,7 @@
 #include <complex.h>
 #include <math.h>
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ double
 __catanh (__complex__ double x)
@@ -76,6 +76,17 @@ __catanh (__complex__ double x)
       den = 1 - __real__ x * __real__ x - i2;
 
       __imag__ res = 0.5 * __ieee754_atan2 (2.0 * __imag__ x, den);
+
+      if (fabs (__real__ res) < DBL_MIN)
+	{
+	  volatile double force_underflow = __real__ res * __real__ res;
+	  (void) force_underflow;
+	}
+      if (fabs (__imag__ res) < DBL_MIN)
+	{
+	  volatile double force_underflow = __imag__ res * __imag__ res;
+	  (void) force_underflow;
+	}
     }
 
   return res;
diff --git a/math/s_catanhf.c b/math/s_catanhf.c
index 9d90a00..6376a2c 100644
--- a/math/s_catanhf.c
+++ b/math/s_catanhf.c
@@ -20,7 +20,7 @@
 #include <complex.h>
 #include <math.h>
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ float
 __catanhf (__complex__ float x)
@@ -76,6 +76,17 @@ __catanhf (__complex__ float x)
       den = 1 - __real__ x * __real__ x - i2;
 
       __imag__ res = 0.5f * __ieee754_atan2f (2.0f * __imag__ x, den);
+
+      if (fabsf (__real__ res) < FLT_MIN)
+	{
+	  volatile float force_underflow = __real__ res * __real__ res;
+	  (void) force_underflow;
+	}
+      if (fabsf (__imag__ res) < FLT_MIN)
+	{
+	  volatile float force_underflow = __imag__ res * __imag__ res;
+	  (void) force_underflow;
+	}
     }
 
   return res;
diff --git a/math/s_catanhl.c b/math/s_catanhl.c
index 7e2b894..30fd277 100644
--- a/math/s_catanhl.c
+++ b/math/s_catanhl.c
@@ -20,7 +20,7 @@
 #include <complex.h>
 #include <math.h>
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ long double
 __catanhl (__complex__ long double x)
@@ -76,6 +76,17 @@ __catanhl (__complex__ long double x)
       den = 1 - __real__ x * __real__ x - i2;
 
       __imag__ res = 0.5L * __ieee754_atan2l (2.0L * __imag__ x, den);
+
+      if (fabsl (__real__ res) < LDBL_MIN)
+	{
+	  volatile long double force_underflow = __real__ res * __real__ res;
+	  (void) force_underflow;
+	}
+      if (fabsl (__imag__ res) < LDBL_MIN)
+	{
+	  volatile long double force_underflow = __imag__ res * __imag__ res;
+	  (void) force_underflow;
+	}
     }
 
   return res;
diff --git a/math/s_catanl.c b/math/s_catanl.c
index 67b6b52..2440e6d 100644
--- a/math/s_catanl.c
+++ b/math/s_catanl.c
@@ -20,7 +20,7 @@
 #include <complex.h>
 #include <math.h>
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ long double
 __catanl (__complex__ long double x)
@@ -83,6 +83,17 @@ __catanl (__complex__ long double x)
 	  num = 4.0L * __imag__ x;
 	  __imag__ res = 0.25L * __log1pl (num / den);
 	}
+
+      if (fabsl (__real__ res) < LDBL_MIN)
+	{
+	  volatile long double force_underflow = __real__ res * __real__ res;
+	  (void) force_underflow;
+	}
+      if (fabsl (__imag__ res) < LDBL_MIN)
+	{
+	  volatile long double force_underflow = __imag__ res * __imag__ res;
+	  (void) force_underflow;
+	}
     }
 
   return res;