diff --git a/math/s_ctanh.c b/math/s_ctanh.c
index fe38dae..b87855f 100644
--- a/math/s_ctanh.c
+++ b/math/s_ctanh.c
@@ -54,25 +54,58 @@ __ctanh (__complex__ double x)
     }
   else
     {
-      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)
+	 |a| >= 0.5*arccosh(1/eps) ~= 18.368
+
+	 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.
+      */
+      if (fabs (__real__ x) < 22.0)
 	{
-	  __complex__ double ez = __cexp (x);
-	  __complex__ double emz = __cexp (-x);
+	  double den;
+
+	  /* avoid cancellation from cos(2y)<0
+	     cosh(2x) + cos(2y) = 2*(sinh(x)^2 + cos(x)^2) */
+	  /* note: actually compute 0.5*den and adjust res */
+	  den = (__ieee754_pow (__ieee754_sinh (__real__ x), 2.0) + __ieee754_pow (__cos (__imag__ x), 2.0));
 
-	  res = (ez - emz) / (ez + emz);
+	  /* sinh(2x) = 2*sinh(x)*cosh(x) */
+	  __real__ res = __ieee754_sinh (__real__ x) * __ieee754_cosh (__real__ x) / den;
+	  __imag__ res = 0.5 * __sin (2.0 * __imag__ x) / den;
 	}
       else
-	{
-	  __real__ res = __ieee754_sinh (2.0 * __real__ x) / den;
-	  __imag__ res = sin2ix / den;
-	}
+        {
+	  // 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;