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 44c9944..91477ee 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 d6f8110..9ca57b2 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 77a9ae3..1d561b4 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 1d7a5a2..d0a95ba 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 4aa9765..41fcea5 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 2568249..d67cc5e 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 602c14e..725989c 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 c1d6a4f..d53f943 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 56bb256..20edbd1 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 ba85e79..4b019a0 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 e482e3a..d9a928a 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 ff43256..9812bdd 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 89c0fef..fdba847 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 2559f83..fd2b797 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 d288b7d..fee1910 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 ca36a83..862845f 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 dbf1612..c4fc1d3 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 4fe2611..0affe93 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)
 	{