Fix ctan, ctanh overflow (bug 11521).

ChangeLog

@@ -1,3 +1,19 @@
+2012-04-09  Joseph Myers  <joseph@codesourcery.com>
+
+	[BZ #11521]
+	* math/s_ctan.c: Include <float.h>.
+	(__ctan): Avoid internal overflow or cancellation in calculating
+	denominator.
+	* math/s_ctanf.c: Likewise.
+	* math/s_ctanl.c: Likewise.
+	* math/s_ctanh.c: Likewise.
+	* math/s_ctanhf.c: Likewise.
+	* math/s_ctanhl.c: Likewise.
+	* math/libm-test.inc (ctan_test): Add more tests.
+	(ctanh_test): Likewise.
+	* sysdeps/i386/fpu/libm-test-ulps: Update.
+	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
+
 2012-04-09  Andreas Jaeger  <aj@suse.de>
 
 	[BZ #6894]

NEWS

@@ -14,13 +14,13 @@ Version 2.16
   4596, 4822, 5077, 5461, 5805, 5993, 6471, 6486, 6578, 6649, 6730, 6770,
   6884, 6890, 6894, 6895, 6907, 6911, 9739, 9902, 10110, 10135, 10140,
   10153, 10210, 10254, 10346, 10545, 10716, 11174, 11322, 11365, 11451,
-  11494, 12047, 12340, 13058, 13525, 13526, 13527, 13528, 13529, 13530,
+  11494, 11521, 12047, 12340, 13058, 13525, 13526, 13527, 13528, 13529,
-  13531, 13532, 13533, 13547, 13551, 13552, 13553, 13555, 13559, 13566,
+  13530, 13531, 13532, 13533, 13547, 13551, 13552, 13553, 13555, 13559,
-  13583, 13592, 13618, 13637, 13656, 13658, 13673, 13691, 13695, 13704,
+  13566, 13583, 13592, 13618, 13637, 13656, 13658, 13673, 13691, 13695,
-  13705, 13706, 13726, 13738, 13760, 13761, 13786, 13792, 13806, 13824,
+  13704, 13705, 13706, 13726, 13738, 13760, 13761, 13786, 13792, 13806,
-  13840, 13841, 13844, 13846, 13851, 13852, 13854, 13871, 13873, 13879,
+  13824, 13840, 13841, 13844, 13846, 13851, 13852, 13854, 13871, 13873,
-  13883, 13892, 13895, 13908, 13910, 13911, 13912, 13913, 13915, 13916,
+  13879, 13883, 13892, 13895, 13908, 13910, 13911, 13912, 13913, 13915,
-  13917, 13918, 13919, 13920, 13921, 13926, 13928, 13938, 13963
+  13916, 13917, 13918, 13919, 13920, 13921, 13926, 13928, 13938, 13963
 
 * ISO C11 support:
 

math/libm-test.inc

@@ -2840,6 +2840,36 @@ ctan_test (void)
   TEST_c_c (ctan, 0.75L, 1.25L, 0.160807785916206426725166058173438663L, 0.975363285031235646193581759755216379L);
   TEST_c_c (ctan, -2, -3, 0.376402564150424829275122113032269084e-2L, -1.00323862735360980144635859782192726L);
 
+  TEST_c_c (ctan, 1, 45, 1.490158918874345552942703234806348520895e-39L, 1.000000000000000000000000000000000000001L);
+  TEST_c_c (ctan, 1, 47, 2.729321264492904590777293425576722354636e-41L, 1.0);
+
+#ifndef TEST_FLOAT
+  TEST_c_c (ctan, 1, 355, 8.140551093483276762350406321792653551513e-309L, 1.0);
+  TEST_c_c (ctan, 1, 365, 1.677892637497921890115075995898773550884e-317L, 1.0);
+#endif
+
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+  TEST_c_c (ctan, 1, 5680, 4.725214596136812019616700920476949798307e-4934L, 1.0);
+  TEST_c_c (ctan, 1, 5690, 9.739393181626937151720816611272607059057e-4943L, 1.0);
+#endif
+
+  TEST_c_c (ctan, 0x3.243f6cp-1, 0, -2.287733242885645987394874673945769518150e7L, 0.0);
+
+  TEST_c_c (ctan, 0x1p127, 1, 0.2446359391192790896381501310437708987204L, 0.9101334047676183761532873794426475906201L);
+
+#ifndef TEST_FLOAT
+  TEST_c_c (ctan, 0x1p1023, 1, -0.2254627924997545057926782581695274244229L, 0.8786063118883068695462540226219865087189L);
+#endif
+
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+  TEST_c_c (ctan, 0x1p16383L, 1, 0.1608598776370396607204448234354670036772L, 0.8133818522051542536316746743877629761488L);
+#endif
+
+  TEST_c_c (ctan, 50000, 50000, plus_zero, 1.0);
+  TEST_c_c (ctan, 50000, -50000, plus_zero, -1.0);
+  TEST_c_c (ctan, -50000, 50000, minus_zero, 1.0);
+  TEST_c_c (ctan, -50000, -50000, minus_zero, -1.0);
+
   END (ctan, complex);
 }
 
@@ -2899,6 +2929,36 @@ ctanh_test (void)
   TEST_c_c (ctanh, 0.75L, 1.25L, 1.37260757053378320258048606571226857L, 0.385795952609750664177596760720790220L);
   TEST_c_c (ctanh, -2, -3, -0.965385879022133124278480269394560686L, 0.988437503832249372031403430350121098e-2L);
 
+  TEST_c_c (ctanh, 45, 1, 1.000000000000000000000000000000000000001L, 1.490158918874345552942703234806348520895e-39L);
+  TEST_c_c (ctanh, 47, 1, 1.0, 2.729321264492904590777293425576722354636e-41L);
+
+#ifndef TEST_FLOAT
+  TEST_c_c (ctanh, 355, 1, 1.0, 8.140551093483276762350406321792653551513e-309L);
+  TEST_c_c (ctanh, 365, 1, 1.0, 1.677892637497921890115075995898773550884e-317L);
+#endif
+
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+  TEST_c_c (ctanh, 5680, 1, 1.0, 4.725214596136812019616700920476949798307e-4934L);
+  TEST_c_c (ctanh, 5690, 1, 1.0, 9.739393181626937151720816611272607059057e-4943L);
+#endif
+
+  TEST_c_c (ctanh, 0, 0x3.243f6cp-1, 0.0, -2.287733242885645987394874673945769518150e7L);
+
+  TEST_c_c (ctanh, 1, 0x1p127, 0.9101334047676183761532873794426475906201L, 0.2446359391192790896381501310437708987204L);
+
+#ifndef TEST_FLOAT
+  TEST_c_c (ctanh, 1, 0x1p1023, 0.8786063118883068695462540226219865087189L, -0.2254627924997545057926782581695274244229L);
+#endif
+
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+  TEST_c_c (ctanh, 1, 0x1p16383L, 0.8133818522051542536316746743877629761488L, 0.1608598776370396607204448234354670036772L);
+#endif
+
+  TEST_c_c (ctanh, 50000, 50000, 1.0, plus_zero);
+  TEST_c_c (ctanh, 50000, -50000, 1.0, minus_zero);
+  TEST_c_c (ctanh, -50000, 50000, -1.0, plus_zero);
+  TEST_c_c (ctanh, -50000, -50000, -1.0, minus_zero);
+
   END (ctanh, complex);
 }
 

math/s_ctan.c

@@ -1,5 +1,5 @@
 /* Complex tangent function for double.
-   Copyright (C) 1997, 2005, 2011 Free Software Foundation, Inc.
+   Copyright (C) 1997-2012 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
@@ -20,9 +20,8 @@
 #include <complex.h>
 #include <fenv.h>
 #include <math.h>
-
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ double
 __ctan (__complex__ double x)
@@ -51,24 +50,45 @@ __ctan (__complex__ double x)
     }
   else
     {
-      double sin2rx, cos2rx;
+      double sinrx, cosrx;
       double den;
+      const int t = (int) ((DBL_MAX_EXP - 1) * M_LN2 / 2);
 
-      __sincos (2.0 * __real__ x, &sin2rx, &cos2rx);
+      /* 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). */
 
-      den = cos2rx + __ieee754_cosh (2.0 * __imag__ x);
+      __sincos (__real__ x, &sinrx, &cosrx);
 
-      if (den == 0.0)
+      if (fabs (__imag__ x) > t)
 	{
-	  __complex__ double ez = __cexp (1.0i * x);
+	  /* Avoid intermediate overflow when the real part of the
-	  __complex__ double emz = __cexp (-1.0i * x);
+	     result may be subnormal.  Ignoring negligible terms, the
+	     imaginary part is +/- 1, the real part is
+	     sin(x)*cos(x)/sinh(y)^2 = 4*sin(x)*cos(x)/exp(2y).  */
+	  double exp_2t = __ieee754_exp (2 * t);
 
-	  res = (ez - emz) / (ez + emz) * -1.0i;
+	  __imag__ res = __copysign (1.0, __imag__ x);
+	  __real__ res = 4 * sinrx * cosrx;
+	  __imag__ x = fabs (__imag__ x);
+	  __imag__ x -= t;
+	  __real__ res /= exp_2t;
+	  if (__imag__ x > t)
+	    {
+	      /* Underflow (original imaginary part of x has absolute
+		 value > 2t).  */
+	      __real__ res /= exp_2t;
+	    }
+	  else
+	    __real__ res /= __ieee754_exp (2 * __imag__ x);
 	}
       else
 	{
-	  __real__ res = sin2rx / den;
+	  double sinhix = __ieee754_sinh (__imag__ x);
-	  __imag__ res = __ieee754_sinh (2.0 * __imag__ x) / den;
+	  double coshix = __ieee754_cosh (__imag__ x);
+
+	  den = cosrx * cosrx + sinhix * sinhix;
+	  __real__ res = sinrx * cosrx / den;
+	  __imag__ res = sinhix * coshix / den;
 	}
     }
 

math/s_ctanf.c

@@ -1,5 +1,5 @@
 /* Complex tangent function for float.
-   Copyright (C) 1997, 2005, 2011 Free Software Foundation, Inc.
+   Copyright (C) 1997-2012 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
@@ -21,7 +21,7 @@
 #include <fenv.h>
 #include <math.h>
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ float
 __ctanf (__complex__ float x)
@@ -50,25 +50,45 @@ __ctanf (__complex__ float x)
     }
   else
     {
-      float sin2rx, cos2rx;
+      float sinrx, cosrx;
       float den;
+      const int t = (int) ((FLT_MAX_EXP - 1) * M_LN2 / 2);
 
-      __sincosf (2.0 * __real__ x, &sin2rx, &cos2rx);
+      /* 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). */
 
-      den = cos2rx + __ieee754_coshf (2.0 * __imag__ x);
+      __sincosf (__real__ x, &sinrx, &cosrx);
 
-
+      if (fabsf (__imag__ x) > t)
-      if (den == 0.0)
 	{
-	  __complex__ float ez = __cexpf (1.0i * x);
+	  /* Avoid intermediate overflow when the real part of the
-	  __complex__ float emz = __cexpf (-1.0i * x);
+	     result may be subnormal.  Ignoring negligible terms, the
+	     imaginary part is +/- 1, the real part is
+	     sin(x)*cos(x)/sinh(y)^2 = 4*sin(x)*cos(x)/exp(2y).  */
+	  float exp_2t = __ieee754_expf (2 * t);
 
-	  res = (ez - emz) / (ez + emz) * -1.0i;
+	  __imag__ res = __copysignf (1.0, __imag__ x);
+	  __real__ res = 4 * sinrx * cosrx;
+	  __imag__ x = fabsf (__imag__ x);
+	  __imag__ x -= t;
+	  __real__ res /= exp_2t;
+	  if (__imag__ x > t)
+	    {
+	      /* Underflow (original imaginary part of x has absolute
+		 value > 2t).  */
+	      __real__ res /= exp_2t;
+	    }
+	  else
+	    __real__ res /= __ieee754_expf (2 * __imag__ x);
 	}
       else
 	{
-	  __real__ res = sin2rx / den;
+	  float sinhix = __ieee754_sinhf (__imag__ x);
-	  __imag__ res = __ieee754_sinhf (2.0 * __imag__ x) / den;
+	  float coshix = __ieee754_coshf (__imag__ x);
+
+	  den = cosrx * cosrx + sinhix * sinhix;
+	  __real__ res = sinrx * cosrx / den;
+	  __imag__ res = sinhix * coshix / den;
 	}
     }
 

math/s_ctanh.c

@@ -1,5 +1,5 @@
 /* Complex hyperbole tangent for double.
-   Copyright (C) 1997, 2005, 2011 Free Software Foundation, Inc.
+   Copyright (C) 1997-2012 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
@@ -21,7 +21,7 @@
 #include <fenv.h>
 #include <math.h>
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ double
 __ctanh (__complex__ double x)
@@ -50,24 +50,45 @@ __ctanh (__complex__ double x)
     }
   else
     {
-      double sin2ix, cos2ix;
+      double sinix, cosix;
       double den;
+      const int t = (int) ((DBL_MAX_EXP - 1) * M_LN2 / 2);
 
-      __sincos (2.0 * __imag__ x, &sin2ix, &cos2ix);
+      /* 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).  */
 
-      den = (__ieee754_cosh (2.0 * __real__ x) + cos2ix);
+      __sincos (__imag__ x, &sinix, &cosix);
 
-      if (den == 0.0)
+      if (fabs (__real__ x) > t)
 	{
-	  __complex__ double ez = __cexp (x);
+	  /* Avoid intermediate overflow when the imaginary part of
-	  __complex__ double emz = __cexp (-x);
+	     the result may be subnormal.  Ignoring negligible terms,
+	     the real part is +/- 1, the imaginary part is
+	     sin(y)*cos(y)/sinh(x)^2 = 4*sin(y)*cos(y)/exp(2x).  */
+	  double exp_2t = __ieee754_exp (2 * t);
 
-	  res = (ez - emz) / (ez + emz);
+	  __real__ res = __copysign (1.0, __real__ x);
+	  __imag__ res = 4 * sinix * cosix;
+	  __real__ x = fabs (__real__ x);
+	  __real__ x -= t;
+	  __imag__ res /= exp_2t;
+	  if (__real__ x > t)
+	    {
+	      /* Underflow (original real part of x has absolute value
+		 > 2t).  */
+	      __imag__ res /= exp_2t;
+	    }
+	  else
+	    __imag__ res /= __ieee754_exp (2 * __real__ x);
 	}
       else
 	{
-	  __real__ res = __ieee754_sinh (2.0 * __real__ x) / den;
+	  double sinhrx = __ieee754_sinh (__real__ x);
-	  __imag__ res = sin2ix / den;
+	  double coshrx = __ieee754_cosh (__real__ x);
+
+	  den = sinhrx * sinhrx + cosix * cosix;
+	  __real__ res = sinhrx * coshrx / den;
+	  __imag__ res = sinix * cosix / den;
 	}
     }
 

math/s_ctanhf.c

@@ -1,5 +1,5 @@
 /* Complex hyperbole tangent for float.
-   Copyright (C) 1997, 2005, 2011 Free Software Foundation, Inc.
+   Copyright (C) 1997-2012 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
@@ -21,7 +21,7 @@
 #include <fenv.h>
 #include <math.h>
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ float
 __ctanhf (__complex__ float x)
@@ -50,24 +50,45 @@ __ctanhf (__complex__ float x)
     }
   else
     {
-      float sin2ix, cos2ix;
+      float sinix, cosix;
       float den;
+      const int t = (int) ((FLT_MAX_EXP - 1) * M_LN2 / 2);
 
-      __sincosf (2.0 * __imag__ x, &sin2ix, &cos2ix);
+      /* 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).  */
 
-      den = (__ieee754_coshf (2.0 * __real__ x) + cos2ix);
+      __sincosf (__imag__ x, &sinix, &cosix);
 
-      if (den == 0.0f)
+      if (fabsf (__real__ x) > t)
 	{
-	  __complex__ float ez = __cexpf (x);
+	  /* Avoid intermediate overflow when the imaginary part of
-	  __complex__ float emz = __cexpf (-x);
+	     the result may be subnormal.  Ignoring negligible terms,
+	     the real part is +/- 1, the imaginary part is
+	     sin(y)*cos(y)/sinh(x)^2 = 4*sin(y)*cos(y)/exp(2x).  */
+	  float exp_2t = __ieee754_expf (2 * t);
 
-	  res = (ez - emz) / (ez + emz);
+	  __real__ res = __copysignf (1.0, __real__ x);
+	  __imag__ res = 4 * sinix * cosix;
+	  __real__ x = fabsf (__real__ x);
+	  __real__ x -= t;
+	  __imag__ res /= exp_2t;
+	  if (__real__ x > t)
+	    {
+	      /* Underflow (original real part of x has absolute value
+		 > 2t).  */
+	      __imag__ res /= exp_2t;
+	    }
+	  else
+	    __imag__ res /= __ieee754_expf (2 * __real__ x);
 	}
       else
 	{
-	  __real__ res = __ieee754_sinhf (2.0 * __real__ x) / den;
+	  float sinhrx = __ieee754_sinhf (__real__ x);
-	  __imag__ res = sin2ix / den;
+	  float coshrx = __ieee754_coshf (__real__ x);
+
+	  den = sinhrx * sinhrx + cosix * cosix;
+	  __real__ res = sinhrx * coshrx / den;
+	  __imag__ res = sinix * cosix / den;
 	}
     }
 

math/s_ctanhl.c

@@ -1,5 +1,5 @@
 /* Complex hyperbole tangent for long double.
-   Copyright (C) 1997, 2005, 2011 Free Software Foundation, Inc.
+   Copyright (C) 1997-2012 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
@@ -21,7 +21,7 @@
 #include <fenv.h>
 #include <math.h>
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ long double
 __ctanhl (__complex__ long double x)
@@ -50,24 +50,45 @@ __ctanhl (__complex__ long double x)
     }
   else
     {
-      long double sin2ix, cos2ix;
+      long double sinix, cosix;
       long double den;
+      const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l / 2);
 
-      __sincosl (2.0 * __imag__ x, &sin2ix, &cos2ix);
+      /* 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).  */
 
-      den = (__ieee754_coshl (2.0 * __real__ x) + cos2ix);
+      __sincosl (__imag__ x, &sinix, &cosix);
 
-      if (den == 0.0L)
+      if (fabsl (__real__ x) > t)
 	{
-	  __complex__ long double ez = __cexpl (x);
+	  /* Avoid intermediate overflow when the imaginary part of
-	  __complex__ long double emz = __cexpl (-x);
+	     the result may be subnormal.  Ignoring negligible terms,
+	     the real part is +/- 1, the imaginary part is
+	     sin(y)*cos(y)/sinh(x)^2 = 4*sin(y)*cos(y)/exp(2x).  */
+	  long double exp_2t = __ieee754_expl (2 * t);
 
-	  res = (ez - emz) / (ez + emz);
+	  __real__ res = __copysignl (1.0, __real__ x);
+	  __imag__ res = 4 * sinix * cosix;
+	  __real__ x = fabsl (__real__ x);
+	  __real__ x -= t;
+	  __imag__ res /= exp_2t;
+	  if (__real__ x > t)
+	    {
+	      /* Underflow (original real part of x has absolute value
+		 > 2t).  */
+	      __imag__ res /= exp_2t;
+	    }
+	  else
+	    __imag__ res /= __ieee754_expl (2 * __real__ x);
 	}
       else
 	{
-	  __real__ res = __ieee754_sinhl (2.0 * __real__ x) / den;
+	  long double sinhrx = __ieee754_sinhl (__real__ x);
-	  __imag__ res = sin2ix / den;
+	  long double coshrx = __ieee754_coshl (__real__ x);
+
+	  den = sinhrx * sinhrx + cosix * cosix;
+	  __real__ res = sinhrx * coshrx / den;
+	  __imag__ res = sinix * cosix / den;
 	}
     }
 

math/s_ctanl.c

@@ -1,5 +1,5 @@
 /* Complex tangent function for long double.
-   Copyright (C) 1997, 2005, 2011 Free Software Foundation, Inc.
+   Copyright (C) 1997-2012 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
@@ -20,9 +20,8 @@
 #include <complex.h>
 #include <fenv.h>
 #include <math.h>
-
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ long double
 __ctanl (__complex__ long double x)
@@ -51,25 +50,45 @@ __ctanl (__complex__ long double x)
     }
   else
     {
-      long double sin2rx, cos2rx;
+      long double sinrx, cosrx;
       long double den;
+      const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l / 2);
 
-      __sincosl (2.0 * __real__ x, &sin2rx, &cos2rx);
+      /* 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). */
 
-      den = cos2rx + __ieee754_coshl (2.0 * __imag__ x);
+      __sincosl (__real__ x, &sinrx, &cosrx);
 
-
+      if (fabsl (__imag__ x) > t)
-      if (den == 0.0)
 	{
-	  __complex__ long double ez = __cexpl (1.0i * x);
+	  /* Avoid intermediate overflow when the real part of the
-	  __complex__ long double emz = __cexpl (-1.0i * x);
+	     result may be subnormal.  Ignoring negligible terms, the
+	     imaginary part is +/- 1, the real part is
+	     sin(x)*cos(x)/sinh(y)^2 = 4*sin(x)*cos(x)/exp(2y).  */
+	  long double exp_2t = __ieee754_expl (2 * t);
 
-	  res = (ez - emz) / (ez + emz) * -1.0i;
+	  __imag__ res = __copysignl (1.0, __imag__ x);
+	  __real__ res = 4 * sinrx * cosrx;
+	  __imag__ x = fabsl (__imag__ x);
+	  __imag__ x -= t;
+	  __real__ res /= exp_2t;
+	  if (__imag__ x > t)
+	    {
+	      /* Underflow (original imaginary part of x has absolute
+		 value > 2t).  */
+	      __real__ res /= exp_2t;
+	    }
+	  else
+	    __real__ res /= __ieee754_expl (2 * __imag__ x);
 	}
       else
 	{
-	  __real__ res = sin2rx / den;
+	  long double sinhix = __ieee754_sinhl (__imag__ x);
-	  __imag__ res = __ieee754_sinhl (2.0 * __imag__ x) / den;
+	  long double coshix = __ieee754_coshl (__imag__ x);
+
+	  den = cosrx * cosrx + sinhix * sinhix;
+	  __real__ res = sinrx * cosrx / den;
+	  __imag__ res = sinhix * coshix / den;
 	}
     }
 

sysdeps/i386/fpu/libm-test-ulps

@@ -1009,9 +1009,33 @@ idouble: 1
 ifloat: 1
 ildouble: 3
 ldouble: 3
+Test "Real part of: ctan (0x1p1023 + 1 i) == -0.2254627924997545057926782581695274244229 + 0.8786063118883068695462540226219865087189 i":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: ctan (0x1p1023 + 1 i) == -0.2254627924997545057926782581695274244229 + 0.8786063118883068695462540226219865087189 i":
+ildouble: 1
+ldouble: 1
+Test "Real part of: ctan (0x1p127 + 1 i) == 0.2446359391192790896381501310437708987204 + 0.9101334047676183761532873794426475906201 i":
+float: 1
+ifloat: 1
+Test "Real part of: ctan (0x3.243f6cp-1 + 0 i) == -2.287733242885645987394874673945769518150e7 + 0.0 i":
+float: 1
+ifloat: 1
+Test "Real part of: ctan (1 + 45 i) == 1.490158918874345552942703234806348520895e-39 + 1.000000000000000000000000000000000000001 i":
+ildouble: 1
+ldouble: 1
+Test "Real part of: ctan (1 + 47 i) == 2.729321264492904590777293425576722354636e-41 + 1.0 i":
+double: 1
+idouble: 1
+ildouble: 2
+ldouble: 2
 
 # ctanh
 Test "Real part of: ctanh (-2 - 3 i) == -0.965385879022133124278480269394560686 + 0.988437503832249372031403430350121098e-2 i":
+float: 1
+ifloat: 1
 ildouble: 3
 ldouble: 3
 Test "Imaginary part of: ctanh (-2 - 3 i) == -0.965385879022133124278480269394560686 + 0.988437503832249372031403430350121098e-2 i":
@@ -1019,9 +1043,14 @@ float: 1
 ifloat: 1
 ildouble: 1
 ldouble: 1
+Test "Imaginary part of: ctanh (0 + 0x3.243f6cp-1 i) == 0.0 - 2.287733242885645987394874673945769518150e7 i":
+float: 1
+ifloat: 1
 Test "Imaginary part of: ctanh (0 + pi/4 i) == 0.0 + 1.0 i":
 float: 1
 ifloat: 1
+ildouble: 1
+ldouble: 1
 Test "Real part of: ctanh (0.75 + 1.25 i) == 1.37260757053378320258048606571226857 + 0.385795952609750664177596760720790220 i":
 double: 1
 float: 1
@@ -1034,6 +1063,25 @@ idouble: 1
 ifloat: 1
 ildouble: 1
 ldouble: 1
+Test "Real part of: ctanh (1 + 0x1p1023 i) == 0.8786063118883068695462540226219865087189 - 0.2254627924997545057926782581695274244229 i":
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: ctanh (1 + 0x1p1023 i) == 0.8786063118883068695462540226219865087189 - 0.2254627924997545057926782581695274244229 i":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: ctanh (1 + 0x1p127 i) == 0.9101334047676183761532873794426475906201 + 0.2446359391192790896381501310437708987204 i":
+float: 1
+ifloat: 1
+Test "Imaginary part of: ctanh (45 + 1 i) == 1.000000000000000000000000000000000000001 + 1.490158918874345552942703234806348520895e-39 i":
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: ctanh (47 + 1 i) == 1.0 + 2.729321264492904590777293425576722354636e-41 i":
+double: 1
+idouble: 1
+ildouble: 2
+ldouble: 2
 
 # erf
 Test "erf (1.25) == 0.922900128256458230136523481197281140":
@@ -2342,33 +2390,35 @@ ldouble: 1
 
 Function: Real part of "ctan":
 double: 1
+float: 1
 idouble: 1
-ildouble: 1
+ifloat: 1
-ldouble: 1
+ildouble: 2
+ldouble: 2
 
 Function: Imaginary part of "ctan":
 double: 1
 float: 1
 idouble: 1
 ifloat: 1
-ildouble: 3
+ildouble: 1
-ldouble: 3
+ldouble: 1
 
 Function: Real part of "ctanh":
 double: 1
 float: 1
 idouble: 1
 ifloat: 1
-ildouble: 3
+ildouble: 1
-ldouble: 3
+ldouble: 1
 
 Function: Imaginary part of "ctanh":
 double: 1
 float: 1
 idouble: 1
 ifloat: 1
-ildouble: 1
+ildouble: 2
-ldouble: 1
+ldouble: 2
 
 Function: "erf":
 double: 1

sysdeps/x86_64/fpu/libm-test-ulps

@@ -1011,11 +1011,15 @@ ldouble: 1
 # ctan
 Test "Real part of: ctan (-2 - 3 i) == 0.376402564150424829275122113032269084e-2 - 1.00323862735360980144635859782192726 i":
 double: 1
+float: 1
 idouble: 1
+ifloat: 1
 ildouble: 1
 ldouble: 1
 Test "Imaginary part of: ctan (-2 - 3 i) == 0.376402564150424829275122113032269084e-2 - 1.00323862735360980144635859782192726 i":
+double: 1
 float: 1
+idouble: 1
 ifloat: 1
 ildouble: 2
 ldouble: 2
@@ -1029,6 +1033,26 @@ idouble: 1
 ifloat: 1
 ildouble: 3
 ldouble: 3
+Test "Real part of: ctan (0x1p1023 + 1 i) == -0.2254627924997545057926782581695274244229 + 0.8786063118883068695462540226219865087189 i":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: ctan (0x1p1023 + 1 i) == -0.2254627924997545057926782581695274244229 + 0.8786063118883068695462540226219865087189 i":
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: ctan (0x1p127 + 1 i) == 0.2446359391192790896381501310437708987204 + 0.9101334047676183761532873794426475906201 i":
+double: 1
+idouble: 1
+Test "Real part of: ctan (0x3.243f6cp-1 + 0 i) == -2.287733242885645987394874673945769518150e7 + 0.0 i":
+float: 1
+ifloat: 1
+Test "Real part of: ctan (1 + 45 i) == 1.490158918874345552942703234806348520895e-39 + 1.000000000000000000000000000000000000001 i":
+ildouble: 1
+ldouble: 1
+Test "Real part of: ctan (1 + 47 i) == 2.729321264492904590777293425576722354636e-41 + 1.0 i":
+ildouble: 2
+ldouble: 2
 
 # ctanh
 Test "Real part of: ctanh (-2 - 3 i) == -0.965385879022133124278480269394560686 + 0.988437503832249372031403430350121098e-2 i":
@@ -1039,19 +1063,51 @@ ifloat: 2
 ildouble: 3
 ldouble: 3
 Test "Imaginary part of: ctanh (-2 - 3 i) == -0.965385879022133124278480269394560686 + 0.988437503832249372031403430350121098e-2 i":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
 ildouble: 1
 ldouble: 1
-Test "Imaginary part of: ctanh (0 + pi/4 i) == 0.0 + 1.0 i":
+Test "Imaginary part of: ctanh (0 + 0x3.243f6cp-1 i) == 0.0 - 2.287733242885645987394874673945769518150e7 i":
 float: 1
 ifloat: 1
+Test "Imaginary part of: ctanh (0 + pi/4 i) == 0.0 + 1.0 i":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
 Test "Real part of: ctanh (0.75 + 1.25 i) == 1.37260757053378320258048606571226857 + 0.385795952609750664177596760720790220 i":
 double: 1
+float: 1
 idouble: 1
+ifloat: 1
 Test "Imaginary part of: ctanh (0.75 + 1.25 i) == 1.37260757053378320258048606571226857 + 0.385795952609750664177596760720790220 i":
 double: 1
+float: 2
+idouble: 1
+ifloat: 2
+ildouble: 1
+ldouble: 1
+Test "Real part of: ctanh (1 + 0x1p1023 i) == 0.8786063118883068695462540226219865087189 - 0.2254627924997545057926782581695274244229 i":
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: ctanh (1 + 0x1p1023 i) == 0.8786063118883068695462540226219865087189 - 0.2254627924997545057926782581695274244229 i":
+double: 1
 idouble: 1
 ildouble: 1
 ldouble: 1
+Test "Real part of: ctanh (1 + 0x1p127 i) == 0.9101334047676183761532873794426475906201 + 0.2446359391192790896381501310437708987204 i":
+double: 1
+idouble: 1
+Test "Imaginary part of: ctanh (45 + 1 i) == 1.000000000000000000000000000000000000001 + 1.490158918874345552942703234806348520895e-39 i":
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: ctanh (47 + 1 i) == 1.0 + 2.729321264492904590777293425576722354636e-41 i":
+ildouble: 2
+ldouble: 2
 
 # erf
 Test "erf (1.25) == 0.922900128256458230136523481197281140":
@@ -2300,34 +2356,36 @@ ldouble: 1
 
 Function: Real part of "ctan":
 double: 1
+float: 1
 idouble: 1
-ildouble: 1
+ifloat: 1
-ldouble: 1
+ildouble: 2
+ldouble: 2
 
 Function: Imaginary part of "ctan":
 double: 1
 float: 1
 idouble: 1
 ifloat: 1
-ildouble: 3
+ildouble: 1
-ldouble: 3
+ldouble: 1
 
 Function: Real part of "ctanh":
 double: 1
-float: 2
-idouble: 1
-ifloat: 2
-ildouble: 3
-ldouble: 3
-
-Function: Imaginary part of "ctanh":
-double: 1
 float: 1
 idouble: 1
 ifloat: 1
 ildouble: 1
 ldouble: 1
 
+Function: Imaginary part of "ctanh":
+double: 1
+float: 2
+idouble: 1
+ifloat: 2
+ildouble: 2
+ldouble: 2
+
 Function: "erf":
 double: 1
 idouble: 1