Fix wrong ussage of sincos for subnormal arguments

18 files changed, 274 insertions, 30 deletions

diff --git a/math/s_ccosh.c b/math/s_ccosh.c
index 44c9944466..91477ee2d0 100644
--- a/math/s_ccosh.c
+++ b/math/s_ccosh.c
@@ -39,7 +39,15 @@ __ccosh (__complex__ double x)
 	  const int t = (int) ((DBL_MAX_EXP - 1) * M_LN2);
 	  double sinix, cosix;
 
-	  __sincos (__imag__ x, &sinix, &cosix);
+	  if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+	    {
+	      __sincos (__imag__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __imag__ x;
+	      cosix = 1.0;
+	    }
 
 	  if (fabs (__real__ x) > t)
 	    {
@@ -92,7 +100,15 @@ __ccosh (__complex__ double x)
 	  /* Imaginary part is finite.  */
 	  double sinix, cosix;
 
-	  __sincos (__imag__ x, &sinix, &cosix);
+	  if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+	    {
+	      __sincos (__imag__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __imag__ x;
+	      cosix = 1.0;
+	    }
 
 	  __real__ retval = __copysign (HUGE_VAL, cosix);
 	  __imag__ retval = (__copysign (HUGE_VAL, sinix)
diff --git a/math/s_ccoshf.c b/math/s_ccoshf.c
index d6f811049e..9ca57b25b9 100644
--- a/math/s_ccoshf.c
+++ b/math/s_ccoshf.c
@@ -39,7 +39,15 @@ __ccoshf (__complex__ float x)
 	  const int t = (int) ((FLT_MAX_EXP - 1) * M_LN2);
 	  float sinix, cosix;
 
-	  __sincosf (__imag__ x, &sinix, &cosix);
+	  if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+	    {
+	      __sincosf (__imag__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __imag__ x;
+	      cosix = 1.0f;
+	    }
 
 	  if (fabsf (__real__ x) > t)
 	    {
@@ -92,7 +100,15 @@ __ccoshf (__complex__ float x)
 	  /* Imaginary part is finite.  */
 	  float sinix, cosix;
 
-	  __sincosf (__imag__ x, &sinix, &cosix);
+	  if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+	    {
+	      __sincosf (__imag__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __imag__ x;
+	      cosix = 1.0f;
+	    }
 
 	  __real__ retval = __copysignf (HUGE_VALF, cosix);
 	  __imag__ retval = (__copysignf (HUGE_VALF, sinix)
diff --git a/math/s_ccoshl.c b/math/s_ccoshl.c
index 77a9ae3869..1d561b4e3f 100644
--- a/math/s_ccoshl.c
+++ b/math/s_ccoshl.c
@@ -39,7 +39,15 @@ __ccoshl (__complex__ long double x)
 	  const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l);
 	  long double sinix, cosix;
 
-	  __sincosl (__imag__ x, &sinix, &cosix);
+	  if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+	    {
+	      __sincosl (__imag__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __imag__ x;
+	      cosix = 1.0;
+	    }
 
 	  if (fabsl (__real__ x) > t)
 	    {
@@ -92,7 +100,15 @@ __ccoshl (__complex__ long double x)
 	  /* Imaginary part is finite.  */
 	  long double sinix, cosix;
 
-	  __sincosl (__imag__ x, &sinix, &cosix);
+	  if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+	    {
+	      __sincosl (__imag__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __imag__ x;
+	      cosix = 1.0;
+	    }
 
 	  __real__ retval = __copysignl (HUGE_VALL, cosix);
 	  __imag__ retval = (__copysignl (HUGE_VALL, sinix)
diff --git a/math/s_cexp.c b/math/s_cexp.c
index 1d7a5a2c40..d0a95ba79c 100644
--- a/math/s_cexp.c
+++ b/math/s_cexp.c
@@ -39,7 +39,15 @@ __cexp (__complex__ double x)
 	  const int t = (int) ((DBL_MAX_EXP - 1) * M_LN2);
 	  double sinix, cosix;
 
-	  __sincos (__imag__ x, &sinix, &cosix);
+	  if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+	    {
+	      __sincos (__imag__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __imag__ x;
+	      cosix = 1.0;
+	    }
 
 	  if (__real__ x > t)
 	    {
@@ -95,7 +103,15 @@ __cexp (__complex__ double x)
 	    {
 	      double sinix, cosix;
 
-	      __sincos (__imag__ x, &sinix, &cosix);
+	      if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+		{
+		  __sincos (__imag__ x, &sinix, &cosix);
+		}
+	      else
+		{
+		  sinix = __imag__ x;
+		  cosix = 1.0;
+		}
 
 	      __real__ retval = __copysign (value, cosix);
 	      __imag__ retval = __copysign (value, sinix);
diff --git a/math/s_cexpf.c b/math/s_cexpf.c
index 4aa9765818..41fcea51db 100644
--- a/math/s_cexpf.c
+++ b/math/s_cexpf.c
@@ -39,7 +39,15 @@ __cexpf (__complex__ float x)
 	  const int t = (int) ((FLT_MAX_EXP - 1) * M_LN2);
 	  float sinix, cosix;
 
-	  __sincosf (__imag__ x, &sinix, &cosix);
+	  if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+	    {
+	      __sincosf (__imag__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __imag__ x;
+	      cosix = 1.0f;
+	    }
 
 	  if (__real__ x > t)
 	    {
@@ -95,7 +103,15 @@ __cexpf (__complex__ float x)
 	    {
 	      float sinix, cosix;
 
-	      __sincosf (__imag__ x, &sinix, &cosix);
+	      if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+		{
+		  __sincosf (__imag__ x, &sinix, &cosix);
+		}
+	      else
+		{
+		  sinix = __imag__ x;
+		  cosix = 1.0f;
+		}
 
 	      __real__ retval = __copysignf (value, cosix);
 	      __imag__ retval = __copysignf (value, sinix);
diff --git a/math/s_cexpl.c b/math/s_cexpl.c
index 256824924f..d67cc5e5eb 100644
--- a/math/s_cexpl.c
+++ b/math/s_cexpl.c
@@ -39,7 +39,15 @@ __cexpl (__complex__ long double x)
 	  const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l);
 	  long double sinix, cosix;
 
-	  __sincosl (__imag__ x, &sinix, &cosix);
+	  if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+	    {
+	      __sincosl (__imag__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __imag__ x;
+	      cosix = 1.0;
+	    }
 
 	  if (__real__ x > t)
 	    {
@@ -95,7 +103,15 @@ __cexpl (__complex__ long double x)
 	    {
 	      long double sinix, cosix;
 
-	      __sincosl (__imag__ x, &sinix, &cosix);
+	      if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+	        {
+		  __sincosl (__imag__ x, &sinix, &cosix);
+	        }
+	      else
+		{
+		  sinix = __imag__ x;
+		  cosix = 1.0;
+		}
 
 	      __real__ retval = __copysignl (value, cosix);
 	      __imag__ retval = __copysignl (value, sinix);
diff --git a/math/s_csin.c b/math/s_csin.c
index 602c14e537..725989c9c7 100644
--- a/math/s_csin.c
+++ b/math/s_csin.c
@@ -42,7 +42,15 @@ __csin (__complex__ double x)
 	  const int t = (int) ((DBL_MAX_EXP - 1) * M_LN2);
 	  double sinix, cosix;
 
-	  __sincos (__real__ x, &sinix, &cosix);
+	  if (__builtin_expect (rcls != FP_SUBNORMAL, 1))
+	    {
+	      __sincos (__real__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __real__ x;
+	      cosix = 1.0;
+	    }
 
 	  if (fabs (__imag__ x) > t)
 	    {
@@ -115,7 +123,15 @@ __csin (__complex__ double x)
 	  /* Real part is finite.  */
 	  double sinix, cosix;
 
-	  __sincos (__real__ x, &sinix, &cosix);
+	  if (__builtin_expect (rcls != FP_SUBNORMAL, 1))
+	    {
+	      __sincos (__real__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __real__ x;
+	      cosix = 1.0;
+	    }
 
 	  __real__ retval = __copysign (HUGE_VAL, sinix);
 	  __imag__ retval = __copysign (HUGE_VAL, cosix);
diff --git a/math/s_csinf.c b/math/s_csinf.c
index c1d6a4f28e..d53f943fa7 100644
--- a/math/s_csinf.c
+++ b/math/s_csinf.c
@@ -42,7 +42,15 @@ __csinf (__complex__ float x)
 	  const int t = (int) ((FLT_MAX_EXP - 1) * M_LN2);
 	  float sinix, cosix;
 
-	  __sincosf (__real__ x, &sinix, &cosix);
+	  if (__builtin_expect (rcls != FP_SUBNORMAL, 1))
+	    {
+	      __sincosf (__real__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __real__ x;
+	      cosix = 1.0f;
+	    }
 
 	  if (fabsf (__imag__ x) > t)
 	    {
@@ -115,7 +123,15 @@ __csinf (__complex__ float x)
 	  /* Real part is finite.  */
 	  float sinix, cosix;
 
-	  __sincosf (__real__ x, &sinix, &cosix);
+	  if (__builtin_expect (rcls != FP_SUBNORMAL, 1))
+	    {
+	      __sincosf (__real__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __real__ x;
+	      cosix = 1.0f;
+	    }
 
 	  __real__ retval = __copysignf (HUGE_VALF, sinix);
 	  __imag__ retval = __copysignf (HUGE_VALF, cosix);
diff --git a/math/s_csinh.c b/math/s_csinh.c
index 56bb256239..20edbd13f0 100644
--- a/math/s_csinh.c
+++ b/math/s_csinh.c
@@ -42,7 +42,15 @@ __csinh (__complex__ double x)
 	  const int t = (int) ((DBL_MAX_EXP - 1) * M_LN2);
 	  double sinix, cosix;
 
-	  __sincos (__imag__ x, &sinix, &cosix);
+	  if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+	    {
+	      __sincos (__imag__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __imag__ x;
+	      cosix = 1.0;
+	    }
 
 	  if (fabs (__real__ x) > t)
 	    {
@@ -109,7 +117,15 @@ __csinh (__complex__ double x)
 	  /* Imaginary part is finite.  */
 	  double sinix, cosix;
 
-	  __sincos (__imag__ x, &sinix, &cosix);
+	  if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+	    {
+	      __sincos (__imag__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __imag__ x;
+	      cosix = 1.0;
+	    }
 
 	  __real__ retval = __copysign (HUGE_VAL, cosix);
 	  __imag__ retval = __copysign (HUGE_VAL, sinix);
diff --git a/math/s_csinhf.c b/math/s_csinhf.c
index ba85e79b8a..4b019a0ed8 100644
--- a/math/s_csinhf.c
+++ b/math/s_csinhf.c
@@ -42,7 +42,15 @@ __csinhf (__complex__ float x)
 	  const int t = (int) ((FLT_MAX_EXP - 1) * M_LN2);
 	  float sinix, cosix;
 
-	  __sincosf (__imag__ x, &sinix, &cosix);
+	  if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+	    {
+	      __sincosf (__imag__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __imag__ x;
+	      cosix = 1.0f;
+	    }
 
 	  if (fabsf (__real__ x) > t)
 	    {
@@ -109,7 +117,15 @@ __csinhf (__complex__ float x)
 	  /* Imaginary part is finite.  */
 	  float sinix, cosix;
 
-	  __sincosf (__imag__ x, &sinix, &cosix);
+	  if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+	    {
+	      __sincosf (__imag__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __imag__ x;
+	      cosix = 1.0f;
+	    }
 
 	  __real__ retval = __copysignf (HUGE_VALF, cosix);
 	  __imag__ retval = __copysignf (HUGE_VALF, sinix);
diff --git a/math/s_csinhl.c b/math/s_csinhl.c
index e482e3a62e..d9a928a401 100644
--- a/math/s_csinhl.c
+++ b/math/s_csinhl.c
@@ -42,7 +42,15 @@ __csinhl (__complex__ long double x)
 	  const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l);
 	  long double sinix, cosix;
 
-	  __sincosl (__imag__ x, &sinix, &cosix);
+	  if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+	    {
+	      __sincosl (__imag__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __imag__ x;
+	      cosix = 1.0;
+	    }
 
 	  if (fabsl (__real__ x) > t)
 	    {
@@ -109,7 +117,15 @@ __csinhl (__complex__ long double x)
 	  /* Imaginary part is finite.  */
 	  long double sinix, cosix;
 
-	  __sincosl (__imag__ x, &sinix, &cosix);
+	  if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+	    {
+	      __sincosl (__imag__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __imag__ x;
+	      cosix = 1.0;
+	    }
 
 	  __real__ retval = __copysignl (HUGE_VALL, cosix);
 	  __imag__ retval = __copysignl (HUGE_VALL, sinix);
diff --git a/math/s_csinl.c b/math/s_csinl.c
index ff43256e54..9812bddd5c 100644
--- a/math/s_csinl.c
+++ b/math/s_csinl.c
@@ -42,7 +42,15 @@ __csinl (__complex__ long double x)
 	  const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l);
 	  long double sinix, cosix;
 
-	  __sincosl (__real__ x, &sinix, &cosix);
+	  if (__builtin_expect (rcls != FP_SUBNORMAL, 1))
+	    {
+	      __sincosl (__real__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __real__ x;
+	      cosix = 1.0;
+	    }
 
 	  if (fabsl (__imag__ x) > t)
 	    {
@@ -115,7 +123,15 @@ __csinl (__complex__ long double x)
 	  /* Real part is finite.  */
 	  long double sinix, cosix;
 
-	  __sincosl (__real__ x, &sinix, &cosix);
+	  if (__builtin_expect (rcls != FP_SUBNORMAL, 1))
+	    {
+	      __sincosl (__real__ x, &sinix, &cosix);
+	    }
+	  else
+	    {
+	      sinix = __real__ x;
+	      cosix = 1.0;
+	    }
 
 	  __real__ retval = __copysignl (HUGE_VALL, sinix);
 	  __imag__ retval = __copysignl (HUGE_VALL, cosix);
diff --git a/math/s_ctan.c b/math/s_ctan.c
index 89c0fef91e..fdba847447 100644
--- a/math/s_ctan.c
+++ b/math/s_ctan.c
@@ -53,11 +53,20 @@ __ctan (__complex__ double x)
       double sinrx, cosrx;
       double den;
       const int t = (int) ((DBL_MAX_EXP - 1) * M_LN2 / 2);
+      int rcls = fpclassify (__real__ x);
 
       /* tan(x+iy) = (sin(2x) + i*sinh(2y))/(cos(2x) + cosh(2y))
 	 = (sin(x)*cos(x) + i*sinh(y)*cosh(y)/(cos(x)^2 + sinh(y)^2). */
 
-      __sincos (__real__ x, &sinrx, &cosrx);
+      if (__builtin_expect (rcls != FP_SUBNORMAL, 1))
+	{
+	  __sincos (__real__ x, &sinrx, &cosrx);
+	}
+      else
+	{
+	  sinrx = __real__ x;
+	  cosrx = 1.0;
+	}
 
       if (fabs (__imag__ x) > t)
 	{
diff --git a/math/s_ctanf.c b/math/s_ctanf.c
index 2559f83f84..fd2b797882 100644
--- a/math/s_ctanf.c
+++ b/math/s_ctanf.c
@@ -57,7 +57,15 @@ __ctanf (__complex__ float x)
       /* tan(x+iy) = (sin(2x) + i*sinh(2y))/(cos(2x) + cosh(2y))
 	 = (sin(x)*cos(x) + i*sinh(y)*cosh(y)/(cos(x)^2 + sinh(y)^2). */
 
-      __sincosf (__real__ x, &sinrx, &cosrx);
+      if (__builtin_expect (fpclassify(__real__ x) != FP_SUBNORMAL, 1))
+	{
+	  __sincosf (__real__ x, &sinrx, &cosrx);
+	}
+      else
+	{
+	  sinrx = __real__ x;
+	  cosrx = 1.0f;
+	}
 
       if (fabsf (__imag__ x) > t)
 	{
diff --git a/math/s_ctanh.c b/math/s_ctanh.c
index d288b7d168..fee1910551 100644
--- a/math/s_ctanh.c
+++ b/math/s_ctanh.c
@@ -53,11 +53,20 @@ __ctanh (__complex__ double x)
       double sinix, cosix;
       double den;
       const int t = (int) ((DBL_MAX_EXP - 1) * M_LN2 / 2);
+      int icls = fpclassify (__imag__ x);
 
       /* tanh(x+iy) = (sinh(2x) + i*sin(2y))/(cosh(2x) + cos(2y))
 	 = (sinh(x)*cosh(x) + i*sin(y)*cos(y))/(sinh(x)^2 + cos(y)^2).  */
 
-      __sincos (__imag__ x, &sinix, &cosix);
+      if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+	{
+	  __sincos (__imag__ x, &sinix, &cosix);
+	}
+      else
+	{
+	  sinix = __imag__ x;
+	  cosix = 1.0;
+	}
 
       if (fabs (__real__ x) > t)
 	{
diff --git a/math/s_ctanhf.c b/math/s_ctanhf.c
index ca36a83bfb..862845f8c2 100644
--- a/math/s_ctanhf.c
+++ b/math/s_ctanhf.c
@@ -57,7 +57,15 @@ __ctanhf (__complex__ float x)
       /* tanh(x+iy) = (sinh(2x) + i*sin(2y))/(cosh(2x) + cos(2y))
 	 = (sinh(x)*cosh(x) + i*sin(y)*cos(y))/(sinh(x)^2 + cos(y)^2).  */
 
-      __sincosf (__imag__ x, &sinix, &cosix);
+      if (__builtin_expect (fpclassify(__imag__ x) != FP_SUBNORMAL, 1))
+	{
+	  __sincosf (__imag__ x, &sinix, &cosix);
+	}
+      else
+	{
+	  sinix = __imag__ x;
+	  cosix = 1.0f;
+	}
 
       if (fabsf (__real__ x) > t)
 	{
diff --git a/math/s_ctanhl.c b/math/s_ctanhl.c
index dbf1612707..c4fc1d3d04 100644
--- a/math/s_ctanhl.c
+++ b/math/s_ctanhl.c
@@ -53,11 +53,20 @@ __ctanhl (__complex__ long double x)
       long double sinix, cosix;
       long double den;
       const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l / 2);
+      int icls = fpclassify (__imag__ x);
 
       /* tanh(x+iy) = (sinh(2x) + i*sin(2y))/(cosh(2x) + cos(2y))
 	 = (sinh(x)*cosh(x) + i*sin(y)*cos(y))/(sinh(x)^2 + cos(y)^2).  */
 
-      __sincosl (__imag__ x, &sinix, &cosix);
+      if (__builtin_expect (icls != FP_SUBNORMAL, 1))
+	{
+	  __sincosl (__imag__ x, &sinix, &cosix);
+	}
+      else
+	{
+	  sinix = __imag__ x;
+	  cosix = 1.0;
+	}
 
       if (fabsl (__real__ x) > t)
 	{
diff --git a/math/s_ctanl.c b/math/s_ctanl.c
index 4fe26119c8..0affe93211 100644
--- a/math/s_ctanl.c
+++ b/math/s_ctanl.c
@@ -53,11 +53,20 @@ __ctanl (__complex__ long double x)
       long double sinrx, cosrx;
       long double den;
       const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l / 2);
+      int rcls = fpclassify (__real__ x);
 
       /* tan(x+iy) = (sin(2x) + i*sinh(2y))/(cos(2x) + cosh(2y))
 	 = (sin(x)*cos(x) + i*sinh(y)*cosh(y)/(cos(x)^2 + sinh(y)^2). */
 
-      __sincosl (__real__ x, &sinrx, &cosrx);
+      if (__builtin_expect (rcls != FP_SUBNORMAL, 1))
+	{
+	  __sincosl (__real__ x, &sinrx, &cosrx);
+	}
+      else
+	{
+	  sinrx = __real__ x;
+	  cosrx = 1.0;
+	}
 
       if (fabsl (__imag__ x) > t)
 	{