Fix ctan, ctanh overflow (bug 11521).

diff --git a/ChangeLog b/ChangeLog
index 109b6a100af6694f16e13faa4e4cd9bb3e623506..3686491a2d6dfbfb48ce5027674cbaa1afc73fc0 100644 (file)
--- a/ChangeLog
+++ b/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]

diff --git a/NEWS b/NEWS
index 2636c22864d774cc42dae3e479b2399dbc47f4c7..bb303f8b9320b2887db51340e3260b9d7ebfc122 100644 (file)
--- a/NEWS
+++ b/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,
-  13531, 13532, 13533, 13547, 13551, 13552, 13553, 13555, 13559, 13566,
-  13583, 13592, 13618, 13637, 13656, 13658, 13673, 13691, 13695, 13704,
-  13705, 13706, 13726, 13738, 13760, 13761, 13786, 13792, 13806, 13824,
-  13840, 13841, 13844, 13846, 13851, 13852, 13854, 13871, 13873, 13879,
-  13883, 13892, 13895, 13908, 13910, 13911, 13912, 13913, 13915, 13916,
-  13917, 13918, 13919, 13920, 13921, 13926, 13928, 13938, 13963
+  11494, 11521, 12047, 12340, 13058, 13525, 13526, 13527, 13528, 13529,
+  13530, 13531, 13532, 13533, 13547, 13551, 13552, 13553, 13555, 13559,
+  13566, 13583, 13592, 13618, 13637, 13656, 13658, 13673, 13691, 13695,
+  13704, 13705, 13706, 13726, 13738, 13760, 13761, 13786, 13792, 13806,
+  13824, 13840, 13841, 13844, 13846, 13851, 13852, 13854, 13871, 13873,
+  13879, 13883, 13892, 13895, 13908, 13910, 13911, 13912, 13913, 13915,
+  13916, 13917, 13918, 13919, 13920, 13921, 13926, 13928, 13938, 13963
 
 * ISO C11 support:
 

diff --git a/math/libm-test.inc b/math/libm-test.inc
index 05334837492a597f128d626c1930c2f7ae5a7f0c..a551b9f3f4626ddb54ecabfe61c20da3df623456 100644 (file)
--- a/math/libm-test.inc
+++ b/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);
 }
 

diff --git a/math/s_ctan.c b/math/s_ctan.c
index c838fadebb05bbff5a25b9e9dd46c9fc41f85274..78117b31031d81ba672683cb60bcfb78e8953908 100644 (file)
--- a/math/s_ctan.c
+++ b/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);
-         __complex__ double emz = __cexp (-1.0i * x);
+         /* Avoid intermediate overflow when the real part of the
+            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;
-         __imag__ res = __ieee754_sinh (2.0 * __imag__ x) / den;
+         double sinhix = __ieee754_sinh (__imag__ x);
+         double coshix = __ieee754_cosh (__imag__ x);
+
+         den = cosrx * cosrx + sinhix * sinhix;
+         __real__ res = sinrx * cosrx / den;
+         __imag__ res = sinhix * coshix / den;
        }
     }
 

diff --git a/math/s_ctanf.c b/math/s_ctanf.c
index 5f7f28ad07daf48a5f2a285efb3e3d866f4359dc..4cba559a44a99cb612f019d5c9a94650f874b309 100644 (file)
--- a/math/s_ctanf.c
+++ b/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);
-
-      den = cos2rx + __ieee754_coshf (2.0 * __imag__ 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 (den == 0.0)
+      if (fabsf (__imag__ x) > t)
        {
-         __complex__ float ez = __cexpf (1.0i * x);
-         __complex__ float emz = __cexpf (-1.0i * x);
+         /* Avoid intermediate overflow when the real part of the
+            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;
-         __imag__ res = __ieee754_sinhf (2.0 * __imag__ x) / den;
+         float sinhix = __ieee754_sinhf (__imag__ x);
+         float coshix = __ieee754_coshf (__imag__ x);
+
+         den = cosrx * cosrx + sinhix * sinhix;
+         __real__ res = sinrx * cosrx / den;
+         __imag__ res = sinhix * coshix / den;
        }
     }
 

diff --git a/math/s_ctanh.c b/math/s_ctanh.c
index 9cecb8bdb77dcc696812eb75e94b7ceb77c850fa..201871e7ec5e54a9cb65845de5f629831563494d 100644 (file)
--- a/math/s_ctanh.c
+++ b/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);
-         __complex__ double emz = __cexp (-x);
+         /* Avoid intermediate overflow when the imaginary part of
+            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;
-         __imag__ res = sin2ix / den;
+         double sinhrx = __ieee754_sinh (__real__ x);
+         double coshrx = __ieee754_cosh (__real__ x);
+
+         den = sinhrx * sinhrx + cosix * cosix;
+         __real__ res = sinhrx * coshrx / den;
+         __imag__ res = sinix * cosix / den;
        }
     }
 

diff --git a/math/s_ctanhf.c b/math/s_ctanhf.c
index fce5aaf290041077ad2e9da12c0c8b972af78ce5..e505155774dd85d96add2ac8aae180a067d3b283 100644 (file)
--- a/math/s_ctanhf.c
+++ b/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);
-         __complex__ float emz = __cexpf (-x);
+         /* Avoid intermediate overflow when the imaginary part of
+            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;
-         __imag__ res = sin2ix / den;
+         float sinhrx = __ieee754_sinhf (__real__ x);
+         float coshrx = __ieee754_coshf (__real__ x);
+
+         den = sinhrx * sinhrx + cosix * cosix;
+         __real__ res = sinhrx * coshrx / den;
+         __imag__ res = sinix * cosix / den;
        }
     }
 

diff --git a/math/s_ctanhl.c b/math/s_ctanhl.c
index d22e13a975a8c0aa75777697d6fc40135c8590fc..e5d677903f26e42735343d191cb71a23d4fbae77 100644 (file)
--- a/math/s_ctanhl.c
+++ b/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);
-         __complex__ long double emz = __cexpl (-x);
+         /* Avoid intermediate overflow when the imaginary part of
+            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;
-         __imag__ res = sin2ix / den;
+         long double sinhrx = __ieee754_sinhl (__real__ x);
+         long double coshrx = __ieee754_coshl (__real__ x);
+
+         den = sinhrx * sinhrx + cosix * cosix;
+         __real__ res = sinhrx * coshrx / den;
+         __imag__ res = sinix * cosix / den;
        }
     }
 

diff --git a/math/s_ctanl.c b/math/s_ctanl.c
index 112dd723d8bc6241820bf34d2ccaf871766b22d3..12d700cad9bc6f0153331682b593175373e5d815 100644 (file)
--- a/math/s_ctanl.c
+++ b/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);
-
-      den = cos2rx + __ieee754_coshl (2.0 * __imag__ 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 (den == 0.0)
+      if (fabsl (__imag__ x) > t)
        {
-         __complex__ long double ez = __cexpl (1.0i * x);
-         __complex__ long double emz = __cexpl (-1.0i * x);
+         /* Avoid intermediate overflow when the real part of the
+            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;
-         __imag__ res = __ieee754_sinhl (2.0 * __imag__ x) / den;
+         long double sinhix = __ieee754_sinhl (__imag__ x);
+         long double coshix = __ieee754_coshl (__imag__ x);
+
+         den = cosrx * cosrx + sinhix * sinhix;
+         __real__ res = sinrx * cosrx / den;
+         __imag__ res = sinhix * coshix / den;
        }
     }
 

diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 1c791405abde61b910a1772b9798ddd76c15e99a..c3a3ce0da2742e3cde687124c1a0e0d6b3878656 100644 (file)
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/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
-ldouble: 1
+ifloat: 1
+ildouble: 2
+ldouble: 2
 
 Function: Imaginary part of "ctan":
 double: 1
 float: 1
 idouble: 1
 ifloat: 1
-ildouble: 3
-ldouble: 3
+ildouble: 1
+ldouble: 1
 
 Function: Real part of "ctanh":
 double: 1
 float: 1
 idouble: 1
 ifloat: 1
-ildouble: 3
-ldouble: 3
+ildouble: 1
+ldouble: 1
 
 Function: Imaginary part of "ctanh":
 double: 1
 float: 1
 idouble: 1
 ifloat: 1
-ildouble: 1
-ldouble: 1
+ildouble: 2
+ldouble: 2
 
 Function: "erf":
 double: 1

diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index dce60cfb5b12f7c29a38818fadb34e6a9df57997..46d858f2981d15ee273beebacffce5edcc6e7b9e 100644 (file)
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/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 + 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
-ldouble: 1
+ifloat: 1
+ildouble: 2
+ldouble: 2
 
 Function: Imaginary part of "ctan":
 double: 1
 float: 1
 idouble: 1
 ifloat: 1
-ildouble: 3
-ldouble: 3
+ildouble: 1
+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