diff --git a/math/s_ctanh.c b/math/s_ctanh.c
index fe38dae..b6a9bc1 100644
--- a/math/s_ctanh.c
+++ b/math/s_ctanh.c
@@ -2,6 +2,7 @@
    Copyright (C) 1997, 2005 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <address@hidden>, 1997.
+   Improved by Judd Storrs and Jaroslav Hajek, 2009.
 
    The GNU C Library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
@@ -52,27 +53,64 @@ __ctanh (__complex__ double x)
 #endif
 	}
     }
-  else
+  else if (fabs (__real__ x) < 22.0)
     {
-      double sin2ix, cos2ix;
-      double den;
+      /* tanh(a+bi) = (sinh(2a) + i*sin(2b)) / (cosh(2a) + cos(2b)) */
+      
+      /* The denominator becomes cosh(2a) for all b to numerical
+	 precision when:
 
-      __sincos (2.0 * __imag__ x, &sin2ix, &cos2ix);
+	 cosh(2a) +/- 1 < cosh(2a)*(1 +/- eps)
 
-      den = (__ieee754_cosh (2.0 * __real__ x) + cos2ix);
+	 This happens for doubles approximately when:
 
-      if (den == 0.0)
-	{
-	  __complex__ double ez = __cexp (x);
-	  __complex__ double emz = __cexp (-x);
+	 |a| >= 0.5*arccosh(1/eps) ~= 18.368
 
-	  res = (ez - emz) / (ez + emz);
-	}
+	 Therefore, when |a| > 18.368
+	 
+	 tanh(a+bi) ~= 1/tanh(a) + i*sin(2b)/cosh(2a)
+
+	 This can be approximated further. Now note that the real
+	 component:
+
+	 1/tanh(18.368) ~= 1 - eps
+
+	 and the imaginary component:
+
+	 |sin(2b)|/cosh(2*18.368) <= 0+/-1/cosh(2*18.368) < 0 +/- eps
+
+	 So, for doubles when a > 18.368 this can be approximated by
+
+	 tanh(a+bi) ~= sign(a) + i*0*sign(sin(2b))
+	 
+	 and maintain eps precision on the real and imaginary
+	 components.  Choose a > 22 for the cutoff only because fdlibm
+	 used 22 as the limiting cutoff for the real-valued tanh.
+      */
+
+      double sinix, cosix;
+      double sinhx, coshx;
+      double den;
+
+      __sincos (__imag__ x, &sinix, &cosix);
+      sinhx = __ieee754_sinh (__real__ x);
+      coshx = __ieee754_cosh (__real__ x);
+
+      den = (sinhx*sinhx + cosix*cosix);
+      if (den != 0.0)
+        {
+          __real__ res = sinhx * coshx / den;
+          __imag__ res = sinix * cosix / den;
+        }
       else
-	{
-	  __real__ res = __ieee754_sinh (2.0 * __real__ x) / den;
-	  __imag__ res = sin2ix / den;
-	}
+        res = x;
+    }
+  else
+    {
+      // x >= 22 return limits
+      // cosh(2x) +/- 1 == cosh(2x)
+      __real__ res = __copysign (1.0, __real__ x);
+      __imag__ res = __copysign (0.0, __sin (2.0 * __imag__ x));
     }
 
   return res;