Index: configure
===================================================================
--- configure	(revision 59184)
+++ configure	(working copy)
@@ -1,5 +1,5 @@
 #! /bin/sh
-# From configure.in Revision: 58653 .
+# From configure.in Revision: 58784 .
 # Guess values for system-dependent variables and create Makefiles.
 # Generated by GNU Autoconf 2.61 for python 2.6.
 #
@@ -15226,11 +15226,14 @@
 
 
 
-for ac_func in alarm bind_textdomain_codeset chflags chown clock confstr \
- ctermid execv fork fpathconf ftime ftruncate \
+
+
+
+for ac_func in alarm asinh bind_textdomain_codeset chflags chown clock \
+ confstr copysign ctermid execv fork fpathconf ftime ftruncate \
  gai_strerror getgroups getlogin getloadavg getpeername getpgid getpid \
  getpriority getpwent getspnam getspent getsid getwd \
- kill killpg lchflags lchown lstat mkfifo mknod mktime \
+ kill killpg lchflags lchown log1p lstat mkfifo mknod mktime \
  mremap nice pathconf pause plock poll pthread_init \
  putenv readlink realpath \
  select setegid seteuid setgid \
Index: configure.in
===================================================================
--- configure.in	(revision 59184)
+++ configure.in	(working copy)
@@ -2303,11 +2303,11 @@
 AC_MSG_RESULT(MACHDEP_OBJS)
 
 # checks for library functions
-AC_CHECK_FUNCS(alarm bind_textdomain_codeset chflags chown clock confstr \
- ctermid execv fork fpathconf ftime ftruncate \
+AC_CHECK_FUNCS(alarm asinh bind_textdomain_codeset chflags chown clock \
+ confstr copysign ctermid execv fork fpathconf ftime ftruncate \
  gai_strerror getgroups getlogin getloadavg getpeername getpgid getpid \
  getpriority getpwent getspnam getspent getsid getwd \
- kill killpg lchflags lchown lstat mkfifo mknod mktime \
+ kill killpg lchflags lchown log1p lstat mkfifo mknod mktime \
  mremap nice pathconf pause plock poll pthread_init \
  putenv readlink realpath \
  select setegid seteuid setgid \
Index: Doc/library/cmath.rst
===================================================================
--- Doc/library/cmath.rst	(revision 59184)
+++ Doc/library/cmath.rst	(working copy)
@@ -14,6 +14,15 @@
 floating-point number, respectively, and the function is then applied to the
 result of the conversion.
 
+.. note::
+
+   On platforms with hardware and system-level support for signed
+   zeros, functions involving branch cuts are continuous on *both*
+   sides of the branch cut: the sign of the zero distinguishes one
+   side of the branch cut from the other.  On platforms that do not
+   support signed zeros the continuity is as specified below.
+
+
 The functions are:
 
 
@@ -37,32 +46,37 @@
 
 .. function:: asinh(x)
 
-   Return the hyperbolic arc sine of *x*. There are two branch cuts, extending
-   left from ``±1j`` to ``±∞j``, both continuous from above. These branch cuts
-   should be considered a bug to be corrected in a future release. The correct
-   branch cuts should extend along the imaginary axis, one from ``1j`` up to
-   ``∞j`` and continuous from the right, and one from ``-1j`` down to ``-∞j``
-   and continuous from the left.
+   Return the hyperbolic arc sine of *x*. There are two branch cuts:
+   One extends from ``1j`` along the imaginary axis to ``∞j``,
+   continuous from the right.  The other extends from ``-1j`` along
+   the imaginary axis to ``-∞j``, continuous from the left.
 
+   .. versionchanged:: 2.6
+      branch cuts moved to match those recommended by the C99 standard
 
+
 .. function:: atan(x)
 
    Return the arc tangent of *x*. There are two branch cuts: One extends from
-   ``1j`` along the imaginary axis to ``∞j``, continuous from the left. The
+   ``1j`` along the imaginary axis to ``∞j``, continuous from the right. The
    other extends from ``-1j`` along the imaginary axis to ``-∞j``, continuous
-   from the left. (This should probably be changed so the upper cut becomes
-   continuous from the other side.)
+   from the left.
 
+   .. versionchanged:: 2.6
+      direction of continuity of upper cut reversed
 
+
 .. function:: atanh(x)
 
    Return the hyperbolic arc tangent of *x*. There are two branch cuts: One
-   extends from ``1`` along the real axis to ``∞``, continuous from above. The
+   extends from ``1`` along the real axis to ``∞``, continuous from below. The
    other extends from ``-1`` along the real axis to ``-∞``, continuous from
-   above. (This should probably be changed so the right cut becomes continuous
-   from the other side.)
+   above.
 
+   .. versionchanged:: 2.6
+      direction of continuity of right cut reversed
 
+
 .. function:: cos(x)
 
    Return the cosine of *x*.
Index: Modules/cmathmodule.c
===================================================================
--- Modules/cmathmodule.c	(revision 59184)
+++ Modules/cmathmodule.c	(working copy)
@@ -5,29 +5,172 @@
 #include "Python.h"
 
 #ifndef M_PI
-#define M_PI (3.141592653589793239)
+#define M_PI (3.141592653589793238)
 #endif
 
-/* First, the C functions that do the real work */
+#ifndef M_PI_2
+#define M_PI_2 (1.570796326794896619)
+#endif
 
-/* constants */
-static Py_complex c_one = {1., 0.};
-static Py_complex c_half = {0.5, 0.};
-static Py_complex c_i = {0., 1.};
-static Py_complex c_halfi = {0., 0.5};
+#ifndef M_E
+#define M_E (2.718281828459045235)
+#endif
 
+#ifndef M_LN2
+#define M_LN2 (0.6931471805599453094) /* natural log of 2 */
+#endif
+
+#ifndef M_LN4
+#define M_LN4 (1.386294361119890619) /* natural log of 4 */
+#endif
+
+#ifndef M_LN10
+#define M_LN10 (2.302585092994045684) /* natural log of 10 */
+#endif
+
+
+/* Machine-dependent constants used as cutoff points for some of the
+   algorithms below.  The values chosen below should work for DEC VAX, IEEE
+   and IBM S/390 double-precision formats.
+
+   LARGE_DOUBLE is a large positive number used to avoid spurious overflow in
+   the sqrt, log, inverse trig and inverse hyperbolic trig functions.  It should be large
+   enough that 1/sqrt(LARGE_DOUBLE) is significantly smaller than the machine
+   epsilon (around 1e-16 on IEEE machines), and small enough that
+   4*LARGE_DOUBLE does not overflow.  SQRT_LARGE_DOUBLE is (approximately) the
+   square root of LARGE_DOUBLE; R_LARGE_DOUBLE is (approximately) the
+   reciprocal of LARGE_DOUBLE.
+
+   EXP_CUTOFF is used in the evaluation of exp, cos, cosh, sin, sinh, tan,
+   tanh to avoid unecessary overflow.  It should be small enough that
+   exp(EXP_CUTOFF) is representable.
+
+   LOG1P_TAYLOR_CUTOFF is used to decide when to use a Taylor-series expansion
+   for the log1p function.  It should be chosen so that LOG1P_TAYLOR_CUTOFF is
+   larger than the machine epsilon, but its square is smaller than the machine
+   epsilon.
+ */
+
+#define LARGE_DOUBLE (1e300)
+#define R_LARGE_DOUBLE (1./LARGE_DOUBLE)
+#define SQRT_LARGE_DOUBLE (1e150)
+#define R_SQRT_LARGE_DOUBLE (1./SQRT_LARGE_DOUBLE)
+#define EXP_CUTOFF (700.)
+#define LOG1P_TAYLOR_CUTOFF (1e-10)
+
+
 /* forward declarations */
-static Py_complex c_log(Py_complex);
-static Py_complex c_prodi(Py_complex);
+static Py_complex c_asinh(Py_complex);
+static Py_complex c_atanh(Py_complex);
+static Py_complex c_cosh(Py_complex);
+static Py_complex c_sinh(Py_complex);
 static Py_complex c_sqrt(Py_complex);
+static Py_complex c_tanh(Py_complex);
 static PyObject * math_error(void);
 
+/* the copysign, log1p and asinh functions are all part of the C99 standard,
+   but not part of the C89 standard.  Here are substitutes for these three
+   functions, to be used if the math library doesn't already have them. */
 
+#ifndef HAVE_COPYSIGN
+static double
+copysign(double x, double y)
+{
+	/* use atan2 to distinguish -0. from 0. */
+	if (y > 0. || (y == 0. && atan2(y, -1.) > 0.)) {
+		return fabs(x);
+	} else {
+		return -fabs(x);
+	}
+}
+#endif
+
+#ifndef HAVE_LOG1P
+static double
+log1p(double x)
+{
+	double y;
+	if (fabs(x) < LOG1P_TAYLOR_CUTOFF) {
+		/* for tiny arguments, use two terms of the Taylor series. */
+		return x*(1.-x/2.);
+	} else if (-0.5 <= x && x <= 1.) {
+		/* WARNING: it's possible than an errant compiler will
+		   incorrectly optimize the following two lines to the
+		   equivalent of "return log(1.+x)". If this happens, then
+		   results from log1p will be inaccurate for small x. */
+		y = 1.+x;
+		return log(y)-((y-1.)-x)/y;
+	} else {
+		/* NaNs and infinities should end up here */
+		return log(1.+x);
+	}
+}
+#endif
+
+#ifndef HAVE_ASINH
+static double
+asinh(double x)
+{
+	/* asinh is an odd function, so we can use
+
+	     asinh(x) = copysign(asinh(abs(x)), x)
+
+	   to reduce to the case when x is nonnegative.  Here, the usual
+	   formula is:
+
+	     asinh(x) = log(x + hypot(1, x))
+
+	   For x small, the argument of the log is close to 1, resulting in
+	   a loss of precision; in this case we use the formula
+
+	     asinh(x) = log1p(x + x*x/(hypot(1, x)+1))
+
+	   For x very large there's a danger of overflow in computing x + hypot(1, x).
+	   In this case we use the approximation
+
+	     asinh(x) ~ log(2*x) = log(x) + log(2).
+
+	*/
+
+	double ax;
+	ax = fabs(x);
+	if (ax > LARGE_DOUBLE) {
+		return copysign(log(ax) + M_LN2, x);
+	} else if (ax < 1.) {
+		return copysign(log1p(ax*(1.+ax/(hypot(1., ax)+1.))), x);
+	} else {
+		return copysign(log(ax + hypot(1., ax)), x);
+	}
+}
+#endif
+
+/* First, the C functions that do the real work */
+
 static Py_complex
-c_acos(Py_complex x)
+c_acos(Py_complex z)
 {
-	return c_neg(c_prodi(c_log(c_sum(x,c_prod(c_i,
-		    c_sqrt(c_diff(c_one,c_prod(x,x))))))));
+	Py_complex s1, s2, r;
+        if (fabs(z.real) > LARGE_DOUBLE || fabs(z.imag) > LARGE_DOUBLE) {
+		/* avoid unnecessary overflow for large arguments */
+		r.real = atan2(fabs(z.imag), z.real);
+		/* split into cases to make sure that the branch cut has the
+		   correct continuity on systems with unsigned zeros */
+		if (z.real < 0.) {
+			r.imag = -copysign(log(hypot(z.real/2., z.imag/2.)) + M_LN4, z.imag);
+		} else {
+			r.imag = copysign(log(hypot(z.real/2., z.imag/2.)) + M_LN4, -z.imag);
+		}
+	} else {
+		s1.real = 1.-z.real;
+		s1.imag = -z.imag;
+		s1 = c_sqrt(s1);
+		s2.real = 1.+z.real;
+		s2.imag = z.imag;
+		s2 = c_sqrt(s2);
+		r.real = 2.*atan2(s1.real, s2.real);
+		r.imag = asinh(s2.real*s1.imag - s2.imag*s1.real);
+	}
+	return r;
 }
 
 PyDoc_STRVAR(c_acos_doc,
@@ -37,13 +180,25 @@
 
 
 static Py_complex
-c_acosh(Py_complex x)
+c_acosh(Py_complex z)
 {
-	Py_complex z;
-	z = c_sqrt(c_half);
-	z = c_log(c_prod(z, c_sum(c_sqrt(c_sum(x,c_one)),
-				  c_sqrt(c_diff(x,c_one)))));
-	return c_sum(z, z);
+	Py_complex s1, s2, r;
+
+        if (fabs(z.real) > LARGE_DOUBLE || fabs(z.imag) > LARGE_DOUBLE) {
+		/* avoid unnecessary overflow for large arguments */
+		r.real = log(hypot(z.real/2., z.imag/2.)) + M_LN4;
+		r.imag = atan2(z.imag, z.real);
+	} else {
+		s1.real = z.real - 1.;
+		s1.imag = z.imag;
+		s1 = c_sqrt(s1);
+		s2.real = z.real + 1.;
+		s2.imag = z.imag;
+		s2 = c_sqrt(s2);
+		r.real = asinh(s1.real*s2.real + s1.imag*s2.imag);
+		r.imag = 2.*atan2(s1.imag, s2.real);
+	}
+	return r;
 }
 
 PyDoc_STRVAR(c_acosh_doc,
@@ -53,14 +208,16 @@
 
 
 static Py_complex
-c_asin(Py_complex x)
+c_asin(Py_complex z)
 {
-	/* -i * log[(sqrt(1-x**2) + i*x] */
-	const Py_complex squared = c_prod(x, x);
-	const Py_complex sqrt_1_minus_x_sq = c_sqrt(c_diff(c_one, squared));
-        return c_neg(c_prodi(c_log(
-        		c_sum(sqrt_1_minus_x_sq, c_prodi(x))
-		    )       )     );
+	/* asin(z) = -i asinh(iz) */
+	Py_complex s, r;
+	s.real = -z.imag;
+	s.imag = z.real;
+	s = c_asinh(s);
+	r.real = s.imag;
+	r.imag = -s.real;
+	return r;
 }
 
 PyDoc_STRVAR(c_asin_doc,
@@ -70,13 +227,28 @@
 
 
 static Py_complex
-c_asinh(Py_complex x)
+c_asinh(Py_complex z)
 {
-	Py_complex z;
-	z = c_sqrt(c_half);
-	z = c_log(c_prod(z, c_sum(c_sqrt(c_sum(x, c_i)),
-				  c_sqrt(c_diff(x, c_i)))));
-	return c_sum(z, z);
+	Py_complex s1, s2, r;
+
+        if (fabs(z.real) > LARGE_DOUBLE || fabs(z.imag) > LARGE_DOUBLE) {
+		if (z.imag >= 0.) {
+			r.real = copysign(log(hypot(z.real/2., z.imag/2.)) + M_LN4, z.real);
+		} else {
+			r.real = -copysign(log(hypot(z.real/2., z.imag/2.)) + M_LN4, -z.real);
+		}
+		r.imag = atan2(z.imag, fabs(z.real));
+	} else {
+		s1.real = 1.+z.imag;
+		s1.imag = -z.real;
+		s1 = c_sqrt(s1);
+		s2.real = 1.-z.imag;
+		s2.imag = z.real;
+		s2 = c_sqrt(s2);
+		r.real = asinh(s1.real*s2.imag-s2.real*s1.imag);
+		r.imag = atan2(z.imag, s1.real*s2.real-s1.imag*s2.imag);
+	}
+	return r;
 }
 
 PyDoc_STRVAR(c_asinh_doc,
@@ -86,9 +258,16 @@
 
 
 static Py_complex
-c_atan(Py_complex x)
+c_atan(Py_complex z)
 {
-	return c_prod(c_halfi,c_log(c_quot(c_sum(c_i,x),c_diff(c_i,x))));
+	/* atan(z) = -i atanh(iz) */
+	Py_complex s, r;
+	s.real = -z.imag;
+	s.imag = z.real;
+	s = c_atanh(s);
+	r.real = s.imag;
+	r.imag = -s.real;
+	return r;
 }
 
 PyDoc_STRVAR(c_atan_doc,
@@ -98,9 +277,43 @@
 
 
 static Py_complex
-c_atanh(Py_complex x)
+c_atanh(Py_complex z)
 {
-	return c_prod(c_half,c_log(c_quot(c_sum(c_one,x),c_diff(c_one,x))));
+	Py_complex r;
+	double ay, ratio;
+
+	/* use oddness to reduce to case where real part is >= 0. */
+	if (z.real < 0.) {
+		return c_neg(c_atanh(c_neg(z)));
+	}
+
+	ay = fabs(z.imag);
+	if (z.real > SQRT_LARGE_DOUBLE || ay > SQRT_LARGE_DOUBLE) {
+		if (ay < z.real) {
+			ratio = ay/z.real;
+			r.real = 1./(1.+ratio*ratio)/z.real;
+		} else {
+			ratio = z.real/ay;
+			r.real = ratio/(1.+ratio*ratio)/ay;
+		}
+		/* the two negations in the next line cancel each other out
+		   except when working with unsigned zeros: they're there to
+		   ensure that the branch cut has the correct continuity on
+		   systems that don't support signed zeros */
+		r.imag = -copysign(M_PI_2, -z.imag);
+	} else if (z.real == 1. && ay < R_SQRT_LARGE_DOUBLE) {
+		/* C99 standard says:  atanh(1+/-0.) should be inf +/- 0i */
+		r.real = -log(sqrt(ay)/sqrt(hypot(ay, 2.)));
+		if (ay == 0.) {
+			r.imag = z.imag;
+		} else {
+			r.imag = copysign(atan2(2., -ay)/2, z.imag);
+		}
+	} else {
+		r.real = log1p(4.*z.real/((1-z.real)*(1-z.real) + ay*ay))/4.;
+		r.imag = -atan2(-2.*z.imag, (1-z.real)*(1+z.real) - ay*ay)/2.;
+	}
+	return r;
 }
 
 PyDoc_STRVAR(c_atanh_doc,
@@ -110,11 +323,13 @@
 
 
 static Py_complex
-c_cos(Py_complex x)
+c_cos(Py_complex z)
 {
+	/* cos(z) = cosh(iz) */
 	Py_complex r;
-	r.real = cos(x.real)*cosh(x.imag);
-	r.imag = -sin(x.real)*sinh(x.imag);
+	r.real = -z.imag;
+	r.imag = z.real;
+	r = c_cosh(r);
 	return r;
 }
 
@@ -125,11 +340,20 @@
 
 
 static Py_complex
-c_cosh(Py_complex x)
+c_cosh(Py_complex z)
 {
 	Py_complex r;
-	r.real = cos(x.imag)*cosh(x.real);
-	r.imag = sin(x.imag)*sinh(x.real);
+	double x_minus_one;
+
+	if (fabs(z.real) > EXP_CUTOFF) {
+		/* avoid corner cases involving unnecessary overflow */
+		x_minus_one = z.real - copysign(1., z.real);
+		r.real = cos(z.imag) * cosh(x_minus_one) * M_E;
+		r.imag = sin(z.imag) * sinh(x_minus_one) * M_E;
+	} else {
+		r.real = cos(z.imag) * cosh(z.real);
+		r.imag = sin(z.imag) * sinh(z.real);
+	}
 	return r;
 }
 
@@ -140,12 +364,20 @@
 
 
 static Py_complex
-c_exp(Py_complex x)
+c_exp(Py_complex z)
 {
 	Py_complex r;
-	double l = exp(x.real);
-	r.real = l*cos(x.imag);
-	r.imag = l*sin(x.imag);
+	double l;
+
+	if (z.real > EXP_CUTOFF) {
+		l = exp(z.real-1.);
+		r.real = l*cos(z.imag)*M_E;
+		r.imag = l*sin(z.imag)*M_E;
+	} else {
+		l = exp(z.real);
+		r.real = l*cos(z.imag);
+		r.imag = l*sin(z.imag);
+	}
 	return r;
 }
 
@@ -156,23 +388,75 @@
 
 
 static Py_complex
-c_log(Py_complex x)
+c_log(Py_complex z)
 {
+	/*
+	   The usual formula for the real part is log(hypot(z.real, z.imag)).
+	   There are four situations where this formula is potentially
+	   problematic:
+
+	   (1) the absolute value of z is subnormal.  Then hypot is subnormal,
+	   so has fewer than the usual number of bits of accuracy, hence may
+	   have large relative error.  This then gives a large relative error
+	   in the log.  (But note that this is not a problem if either z.real
+	   or z.imag is zero, in which case hypot(z.real, z.imag) is
+	   represented exactly.)
+
+	   (2) The absolute value of z overflows (e.g. when both z.real and
+	   z.imag are within a factor of 1/sqrt(2) of the largest
+	   representable double).
+
+	   (3) the absolute value of z is close to 1.  In this case it's
+	   unreasonable to expect good accuracy: a change of 1ulp in the real
+	   or imaginary part of z can result in a change of billions of ulps
+	   in the correctly rounded answer.  But there's one special case
+	   where good accuracy *is* achievable with little effort, namely when
+	   fabs(z.real) == 1 or fabs(z.imag) == 1.
+
+	   (4) z = 0.  The simplest thing to do here is to call the
+	   floating-point log with an argument of 0, and let its behaviour
+	   (returning -infinity, signaling a floating-point exception, setting
+	   errno, or whatever) determine that of c_log.  So the usual formula
+	   is fine here.
+
+	   All of situations (1)-(3) above can be resolved by using the
+	   formula
+
+	     log(hypot(x, y)) = log(y) + log1p((x/y)*(x/y))/2
+
+           for 0 < x <= y.
+
+	 */
+
 	Py_complex r;
-	double l = hypot(x.real,x.imag);
-	r.imag = atan2(x.imag, x.real);
-	r.real = log(l);
+	double ax, ay, am, ratio;
+
+	ax = fabs(z.real);
+	ay = fabs(z.imag);
+	if ((0. < ax && ax < R_LARGE_DOUBLE && 0. < ay && ay < R_LARGE_DOUBLE) ||
+	    ax > LARGE_DOUBLE ||
+	    ay > LARGE_DOUBLE ||
+	    ax == 1. ||
+	    ay == 1.) {
+		am = (ax > ay ? ax : ay);        /* am = MAX(ax, ay) */
+		ratio = (ax > ay ? ay : ax)/am;  /* ratio = MIN(ax, ay)/MAX(ax, ay) */
+		r.real = log(am) + log1p(ratio*ratio)/2.;
+	} else {
+		r.real = log(hypot(ax, ay));
+	}
+	r.imag = atan2(z.imag, z.real);
 	return r;
 }
 
 
 static Py_complex
-c_log10(Py_complex x)
+c_log10(Py_complex z)
 {
 	Py_complex r;
-	double l = hypot(x.real,x.imag);
-	r.imag = atan2(x.imag, x.real)/log(10.);
-	r.real = log10(l);
+
+	r = c_log(z);
+	r.real = r.real / M_LN10;
+	r.imag = r.imag / M_LN10;
 	return r;
 }
 
@@ -182,26 +466,19 @@
 "Return the base-10 logarithm of x.");
 
 
-/* internal function not available from Python */
 static Py_complex
-c_prodi(Py_complex x)
+c_sin(Py_complex z)
 {
-	Py_complex r;
-	r.real = -x.imag;
-	r.imag = x.real;
+	/* sin(z) = -i sin(iz) */
+	Py_complex s, r;
+	s.real = -z.imag;
+	s.imag = z.real;
+	s = c_sinh(s);
+	r.real = s.imag;
+	r.imag = -s.real;
 	return r;
 }
 
-
-static Py_complex
-c_sin(Py_complex x)
-{
-	Py_complex r;
-	r.real = sin(x.real) * cosh(x.imag);
-	r.imag = cos(x.real) * sinh(x.imag);
-	return r;
-}
-
 PyDoc_STRVAR(c_sin_doc,
 "sin(x)\n"
 "\n"
@@ -209,12 +486,21 @@
 
 
 static Py_complex
-c_sinh(Py_complex x)
+c_sinh(Py_complex z)
 {
 	Py_complex r;
-	r.real = cos(x.imag) * sinh(x.real);
-	r.imag = sin(x.imag) * cosh(x.real);
+	double x_minus_one;
+
+	if (fabs(z.real) > EXP_CUTOFF) {
+		x_minus_one = z.real - copysign(1., z.real);
+		r.real = cos(z.imag) * sinh(x_minus_one) * M_E;
+		r.imag = sin(z.imag) * cosh(x_minus_one) * M_E;
+	} else {
+		r.real = cos(z.imag) * sinh(z.real);
+		r.imag = sin(z.imag) * cosh(z.real);
+	}
 	return r;
+
 }
 
 PyDoc_STRVAR(c_sinh_doc,
@@ -224,28 +510,66 @@
 
 
 static Py_complex
-c_sqrt(Py_complex x)
+c_sqrt(Py_complex z)
 {
+	/*
+	   Method: use symmetries to reduce to the case when x = z.real and y
+	   = z.imag are nonnegative.  Then the real part of the result is
+	   given by
+
+	     s = sqrt((x + hypot(x, y))/2)
+
+	   and the imaginary part is
+
+	     d = (y/2)/s
+
+	   If either x or y is very large then there's a risk of overflow in
+	   computation of the expression x + hypot(x, y).  We can avoid this
+	   by rewriting the formula for s as:
+
+	     s = 2*sqrt(x/8 + hypot(x/8, y/8))
+
+	   This costs us two extra multiplications, but avoids the overhead of
+	   checking for x and y large.
+
+	   If both x and y are subnormal then hypot(x, y) may also be
+	   subnormal, so will lack full precision.  We solve this by rescaling
+	   x and y by a sufficiently large power of 2 to ensure that x and y
+	   are normal.
+	*/
+
+
 	Py_complex r;
 	double s,d;
-	if (x.real == 0. && x.imag == 0.)
-		r = x;
-	else {
-		s = sqrt(0.5*(fabs(x.real) + hypot(x.real,x.imag)));
-		d = 0.5*x.imag/s;
-		if (x.real > 0.) {
-			r.real = s;
-			r.imag = d;
-		}
-		else if (x.imag >= 0.) {
-			r.real = d;
-			r.imag = s;
-		}
-		else {
-			r.real = -d;
-			r.imag = -s;
-		}
+	double ax, ay;
+
+	if (z.real == 0. && z.imag == 0.) {
+		r.real = 0.;
+		r.imag = z.imag;
+		return r;
 	}
+
+	ax = fabs(z.real);
+	ay = fabs(z.imag);
+
+	if (ax < R_LARGE_DOUBLE && ay < R_LARGE_DOUBLE) {
+		/* here we catch cases where hypot(ax, ay) is subnormal */
+		ax = ldexp(ax, 65);
+		s = ldexp(sqrt(ax + hypot(ax, ldexp(ay, 65))), -33);
+		d = ay/(2.*s);
+	} else {
+		ax /= 8.;
+		s = 2.*sqrt(ax + hypot(ax, ay/8.));
+		d = ay/(2.*s);
+	}
+
+	if (z.real >= 0.) {
+		r.real = s;
+		r.imag = copysign(d, z.imag);
+	} else {
+		r.real = d;
+		r.imag = copysign(s, z.imag);
+	}
 	return r;
 }
 
@@ -256,23 +580,15 @@
 
 
 static Py_complex
-c_tan(Py_complex x)
+c_tan(Py_complex z)
 {
-	Py_complex r;
-	double sr,cr,shi,chi;
-	double rs,is,rc,ic;
-	double d;
-	sr = sin(x.real);
-	cr = cos(x.real);
-	shi = sinh(x.imag);
-	chi = cosh(x.imag);
-	rs = sr * chi;
-	is = cr * shi;
-	rc = cr * chi;
-	ic = -sr * shi;
-	d = rc*rc + ic * ic;
-	r.real = (rs*rc + is*ic) / d;
-	r.imag = (is*rc - rs*ic) / d;
+	/* tan(z) = -i tanh(iz) */
+	Py_complex s, r;
+	s.real = -z.imag;
+	s.imag = z.real;
+	s = c_tanh(s);
+	r.real = s.imag;
+	r.imag = -s.real;
 	return r;
 }
 
@@ -283,23 +599,35 @@
 
 
 static Py_complex
-c_tanh(Py_complex x)
+c_tanh(Py_complex z)
 {
+	/* Formula:
+
+	   tanh(x+iy) = (tanh(x)(1+tan(y)^2) + i tan(y)(1-tanh(x))^2) /
+	   (1+tan(y)^2 tanh(x)^2)
+
+	   To avoid excessive roundoff error, 1-tanh(x)^2 is better computed
+	   as 1/cosh(x)^2.  When abs(x) is large, we approximate 1-tanh(x)^2
+	   by 4 exp(-2*x) instead, to avoid possible overflow in the
+	   computation of cosh(x).
+
+	*/
+
 	Py_complex r;
-	double si,ci,shr,chr;
-	double rs,is,rc,ic;
-	double d;
-	si = sin(x.imag);
-	ci = cos(x.imag);
-	shr = sinh(x.real);
-	chr = cosh(x.real);
-	rs = ci * shr;
-	is = si * chr;
-	rc = ci * chr;
-	ic = si * shr;
-	d = rc*rc + ic*ic;
-	r.real = (rs*rc + is*ic) / d;
-	r.imag = (is*rc - rs*ic) / d;
+	double tx, ty, cx, txty, denom;
+
+	if (fabs(z.real) > EXP_CUTOFF) {
+		r.real = copysign(1., z.real);
+                r.imag = 2.*sin(2.*z.imag)*exp(-2.*fabs(z.real));
+	} else {
+		tx = tanh(z.real);
+		ty = tan(z.imag);
+		cx = 1./cosh(z.real);
+		txty = tx*ty;
+		denom = 1. + txty*txty;
+		r.real = tx*(1.+ty*ty)/denom;
+		r.imag = ((ty/denom)*cx)*cx;
+	}
 	return r;
 }
 
@@ -323,9 +651,9 @@
 	if (PyTuple_GET_SIZE(args) == 2)
 		x = c_quot(x, c_log(y));
 	PyFPE_END_PROTECT(x)
+	Py_ADJUST_ERANGE2(x.real, x.imag);
 	if (errno != 0)
 		return math_error();
-	Py_ADJUST_ERANGE2(x.real, x.imag);
 	return PyComplex_FromCComplex(x);
 }
 
Index: pyconfig.h.in
===================================================================
--- pyconfig.h.in	(revision 59184)
+++ pyconfig.h.in	(working copy)
@@ -37,6 +37,9 @@
 /* Define this if your time.h defines altzone. */
 #undef HAVE_ALTZONE
 
+/* Define to 1 if you have the `asinh' function. */
+#undef HAVE_ASINH
+
 /* Define to 1 if you have the <asm/types.h> header file. */
 #undef HAVE_ASM_TYPES_H
 
@@ -85,6 +88,9 @@
 /* Define to 1 if you have the <conio.h> header file. */
 #undef HAVE_CONIO_H
 
+/* Define to 1 if you have the `copysign' function. */
+#undef HAVE_COPYSIGN
+
 /* Define to 1 if you have the `ctermid' function. */
 #undef HAVE_CTERMID
 
@@ -333,6 +339,9 @@
 /* Define to 1 if you have the <linux/netlink.h> header file. */
 #undef HAVE_LINUX_NETLINK_H
 
+/* Define to 1 if you have the `log1p' function. */
+#undef HAVE_LOG1P
+
 /* Define this if you have the type long long. */
 #undef HAVE_LONG_LONG