diff --git a/Modules/cmathmodule.c b/Modules/cmathmodule.c
index ec48ce8..8e3c31e 100644
--- a/Modules/cmathmodule.c
+++ b/Modules/cmathmodule.c
@@ -3,31 +3,172 @@
 /* much code borrowed from mathmodule.c */
 
 #include "Python.h"
+/* we need DBL_MAX, DBL_MIN, DBL_EPSILON, DBL_MANT_DIG and FLT_RADIX from
+   float.h.  We assume that FLT_RADIX is either 2 or 16. */
+#include <float.h>
 
-#ifndef M_PI
-#define M_PI (3.141592653589793239)
+#if (FLT_RADIX != 2 && FLT_RADIX != 16)
+#error "Modules/cmathmodule.c expects FLT_RADIX to be 2 or 16"
 #endif
 
-/* First, the C functions that do the real work */
+#ifndef M_LN2
+#define M_LN2 (0.6931471805599453094) /* natural log of 2 */
+#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_LN10
+#define M_LN10 (2.302585092994045684) /* natural log of 10 */
+#endif
+
+/*
+   CM_LARGE_DOUBLE is used to avoid spurious overflow in the sqrt, log,
+   inverse trig and inverse hyperbolic trig functions.  Its log is used in the
+   evaluation of exp, cos, cosh, sin, sinh, tan, and tanh to avoid unecessary
+   overflow.
+ */
+
+#define CM_LARGE_DOUBLE (DBL_MAX/4.)
+#define CM_SQRT_LARGE_DOUBLE (sqrt(CM_LARGE_DOUBLE))
+#define CM_LOG_LARGE_DOUBLE (log(CM_LARGE_DOUBLE))
+#define CM_SQRT_DBL_MIN (sqrt(DBL_MIN))
+
+/* 
+   CM_SCALE_UP is an odd integer chosen such that multiplication by
+   2**CM_SCALE_UP is sufficient to turn a subnormal into a normal.
+   CM_SCALE_DOWN is (-(CM_SCALE_UP+1)/2).  These scalings are used to compute
+   square roots accurately when the real and imaginary parts of the argument
+   are subnormal.
+*/
+
+#if FLT_RADIX==2
+#define CM_SCALE_UP (2*(DBL_MANT_DIG/2) + 1)
+#elif FLT_RADIX==16
+#define CM_SCALE_UP (4*DBL_MANT_DIG+1)
+#endif
+#define CM_SCALE_DOWN (-(CM_SCALE_UP+1)/2)
 
 /* 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);
 
+/* Code to deal with special values (infinities, NaNs, etc.). */
+
+/* special_type takes a double and returns an integer code indicating
+   the type of the double as follows:
+*/
+
+enum special_types {
+	ST_NINF,	/* 0, negative infinity */
+	ST_NEG,		/* 1, negative finite number (nonzero) */
+	ST_NZERO,	/* 2, -0. */
+	ST_PZERO,	/* 3, +0. */
+	ST_POS,		/* 4, positive finite number (nonzero) */
+	ST_PINF,	/* 5, positive infinity */
+	ST_NAN,		/* 6, Not a Number */
+};
+
+static enum special_types
+special_type(double d)
+{
+	if (Py_IS_FINITE(d)) {
+		if (d != 0) {
+			if (copysign(1., d) == 1.)
+				return ST_POS;
+			else
+				return ST_NEG;
+		}
+		else {
+			if (copysign(1., d) == 1.)
+				return ST_PZERO;
+			else
+				return ST_NZERO;
+		}
+	}
+	if (Py_IS_NAN(d))
+		return ST_NAN;
+	if (copysign(1., d) == 1.)
+		return ST_PINF;
+	else
+		return ST_NINF;
+}
+
+#define SPECIAL_VALUE(z, table)						\
+	if (!Py_IS_FINITE((z).real) || !Py_IS_FINITE((z).imag)) {	\
+		errno = 0;                                              \
+		return table[special_type((z).real)]	                \
+			    [special_type((z).imag)];			\
+	}
+
+#define P Py_MATH_PI
+#define P14 0.25*Py_MATH_PI
+#define P12 0.5*Py_MATH_PI
+#define P34 0.75*Py_MATH_PI
+#ifdef MS_WINDOWS
+/* On Windows HUGE_VAL is an extern variable and not a constant. Since the
+   special value arrays need a constant we have to roll our own infinity
+   and nan. */
+#  define INF (DBL_MAX*DBL_MAX)
+#  define N (INF*0.)
+#else
+#  define INF Py_HUGE_VAL
+#  define N Py_NAN
+#endif /* MS_WINDOWS */
+#define U -9.5426319407711027e33 /* unlikely value, used as placeholder */
+
+/* First, the C functions that do the real work.  Each of the c_*
+   functions computes and returns the C99 Annex G recommended result
+   and also sets errno as follows: errno = 0 if no floating-point
+   exception is associated with the result; errno = EDOM if C99 Annex
+   G recommends raising divide-by-zero or invalid for this result; and
+   errno = ERANGE where the overflow floating-point signal should be
+   raised.
+*/
+
+static Py_complex acos_special_values[7][7] = {
+  {{P34,INF},{P,INF}, {P,INF}, {P,-INF}, {P,-INF}, {P34,-INF},{N,INF}},
+  {{P12,INF},{U,U},   {U,U},   {U,U},    {U,U},    {P12,-INF},{N,N}},
+  {{P12,INF},{U,U},   {P12,0.},{P12,-0.},{U,U},    {P12,-INF},{P12,N}},
+  {{P12,INF},{U,U},   {P12,0.},{P12,-0.},{U,U},    {P12,-INF},{P12,N}},
+  {{P12,INF},{U,U},   {U,U},   {U,U},    {U,U},    {P12,-INF},{N,N}},
+  {{P14,INF},{0.,INF},{0.,INF},{0.,-INF},{0.,-INF},{P14,-INF},{N,INF}},
+  {{N,INF},  {N,N},   {N,N},   {N,N},    {N,N},    {N,-INF},  {N,N}}
+};
 
 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;
+
+	SPECIAL_VALUE(z, acos_special_values);
+
+	if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_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_LN2*2., z.imag);
+		} else {
+			r.imag = copysign(log(hypot(z.real/2., z.imag/2.)) +
+					  M_LN2*2., -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);
+	}
+	errno = 0;
+	return r;
 }
 
 PyDoc_STRVAR(c_acos_doc,
@@ -36,14 +177,39 @@
 "Return the arc cosine of x.");
 
 
+static Py_complex acosh_special_values[7][7] = {
+  {{INF,-P34},{INF,-P}, {INF,-P}, {INF,P}, {INF,P}, {INF,P34},{INF,N}},
+  {{INF,-P12},{U,U},    {U,U},    {U,U},   {U,U},   {INF,P12},{N,N}},
+  {{INF,-P12},{U,U},    {0.,-P12},{0.,P12},{U,U},   {INF,P12},{N,N}},
+  {{INF,-P12},{U,U},    {0.,-P12},{0.,P12},{U,U},   {INF,P12},{N,N}},
+  {{INF,-P12},{U,U},    {U,U},    {U,U},   {U,U},   {INF,P12},{N,N}},
+  {{INF,-P14},{INF,-0.},{INF,-0.},{INF,0.},{INF,0.},{INF,P14},{INF,N}},
+  {{INF,N},   {N,N},    {N,N},    {N,N},   {N,N},   {INF,N},  {N,N}}
+};
+
 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;
+
+	SPECIAL_VALUE(z, acosh_special_values);
+
+	if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {
+		/* avoid unnecessary overflow for large arguments */
+		r.real = log(hypot(z.real/2., z.imag/2.)) + M_LN2*2.;
+		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);
+	}
+	errno = 0;
+	return r;
 }
 
 PyDoc_STRVAR(c_acosh_doc,
@@ -53,14 +219,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,
@@ -69,14 +237,44 @@
 "Return the arc sine of x.");
 
 
+static Py_complex asinh_special_values[7][7] = {
+  {{-INF,-P14},{-INF,-0.},{-INF,-0.},{-INF,0.},{-INF,0.},{-INF,P14},{-INF,N}},
+  {{-INF,-P12},{U,U},     {U,U},     {U,U},    {U,U},    {-INF,P12},{N,N}},
+  {{-INF,-P12},{U,U},     {-0.,-0.}, {-0.,0.}, {U,U},    {-INF,P12},{N,N}},
+  {{INF,-P12}, {U,U},     {0.,-0.},  {0.,0.},  {U,U},    {INF,P12}, {N,N}},
+  {{INF,-P12}, {U,U},     {U,U},     {U,U},    {U,U},    {INF,P12}, {N,N}},
+  {{INF,-P14}, {INF,-0.}, {INF,-0.}, {INF,0.}, {INF,0.}, {INF,P14}, {INF,N}},
+  {{INF,N},    {N,N},     {N,-0.},   {N,0.},   {N,N},    {INF,N},   {N,N}}
+};
+
 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;
+
+	SPECIAL_VALUE(z, asinh_special_values);
+
+	if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {
+		if (z.imag >= 0.) {
+			r.real = copysign(log(hypot(z.real/2., z.imag/2.)) +
+					  M_LN2*2., z.real);
+		} else {
+			r.real = -copysign(log(hypot(z.real/2., z.imag/2.)) +
+					   M_LN2*2., -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);
+	}
+	errno = 0;
+	return r;
 }
 
 PyDoc_STRVAR(c_asinh_doc,
@@ -86,9 +284,37 @@
 
 
 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;
+}
+
+/* Windows screws up atan2 for inf and nan */
+static double
+c_atan2(Py_complex z)
+{
+	if (Py_IS_NAN(z.real) || Py_IS_NAN(z.imag))
+		return Py_NAN;
+	if (Py_IS_INFINITY(z.imag)) {
+		if (Py_IS_INFINITY(z.real)) {
+			if (copysign(1., z.real) == 1.)
+				/* atan2(+-inf, +inf) == +-pi/4 */
+				return copysign(0.25*Py_MATH_PI, z.imag);
+			else
+				/* atan2(+-inf, -inf) == +-pi*3/4 */
+				return copysign(0.75*Py_MATH_PI, z.imag);
+		}
+		/* atan2(+-inf, x) == +-pi/2 for finite x */
+		return copysign(0.5*Py_MATH_PI, z.imag);
+	}
+	return atan2(z.imag, z.real);
 }
 
 PyDoc_STRVAR(c_atan_doc,
@@ -97,10 +323,61 @@
 "Return the arc tangent of x.");
 
 
+static Py_complex atanh_special_values[7][7] = {
+  {{-0.,-P12},{-0.,-P12},{-0.,-P12},{-0.,P12},{-0.,P12},{-0.,P12},{-0.,N}},
+  {{-0.,-P12},{U,U},     {U,U},     {U,U},    {U,U},    {-0.,P12},{N,N}},
+  {{-0.,-P12},{U,U},     {-0.,-0.}, {-0.,0.}, {U,U},    {-0.,P12},{-0.,N}},
+  {{0.,-P12}, {U,U},     {0.,-0.},  {0.,0.},  {U,U},    {0.,P12}, {0.,N}},
+  {{0.,-P12}, {U,U},     {U,U},     {U,U},    {U,U},    {0.,P12}, {N,N}},
+  {{0.,-P12}, {0.,-P12}, {0.,-P12}, {0.,P12}, {0.,P12}, {0.,P12}, {0.,N}},
+  {{0.,-P12}, {N,N},     {N,N},     {N,N},    {N,N},    {0.,P12}, {N,N}}
+};
+
 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, h;
+
+	SPECIAL_VALUE(z, atanh_special_values);
+
+	/* Reduce to case where z.real >= 0., using atanh(z) = -atanh(-z). */
+	if (z.real < 0.) {
+		return c_neg(c_atanh(c_neg(z)));
+	}
+
+	ay = fabs(z.imag);
+	if (z.real > CM_SQRT_LARGE_DOUBLE || ay > CM_SQRT_LARGE_DOUBLE) {
+		/*
+		   if abs(z) is large then we use the approximation
+		   atanh(z) ~ 1/z +/- i*pi/2 (+/- depending on the sign
+		   of z.imag)
+		*/
+		h = hypot(z.real/2., z.imag/2.);  /* safe from overflow */
+		r.real = z.real/4./h/h;
+		/* 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(Py_MATH_PI/2., -z.imag);
+		errno = 0;
+	} else if (z.real == 1. && ay < CM_SQRT_DBL_MIN) {
+		/* C99 standard says:  atanh(1+/-0.) should be inf +/- 0i */
+		if (ay == 0.) {
+			r.real = INF;
+			r.imag = z.imag;
+			errno = EDOM;
+		} else {
+			r.real = -log(sqrt(ay)/sqrt(hypot(ay, 2.)));
+			r.imag = copysign(atan2(2., -ay)/2, z.imag);
+			errno = 0;
+		}
+	} 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.;
+		errno = 0;
+	}
+	return r;
 }
 
 PyDoc_STRVAR(c_atanh_doc,
@@ -110,11 +387,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;
 }
 
@@ -124,12 +403,64 @@
 "Return the cosine of x.");
 
 
+/* cosh(infinity + i*y) needs to be dealt with specially */
+static Py_complex cosh_special_values[7][7] = {
+  {{INF,N},{U,U},{INF,0.}, {INF,-0.},{U,U},{INF,N},{INF,N}},
+  {{N,N},  {U,U},{U,U},    {U,U},    {U,U},{N,N},  {N,N}},
+  {{N,0.}, {U,U},{1.,0.},  {1.,-0.}, {U,U},{N,0.}, {N,0.}},
+  {{N,0.}, {U,U},{1.,-0.}, {1.,0.},  {U,U},{N,0.}, {N,0.}},
+  {{N,N},  {U,U},{U,U},    {U,U},    {U,U},{N,N},  {N,N}},
+  {{INF,N},{U,U},{INF,-0.},{INF,0.}, {U,U},{INF,N},{INF,N}},
+  {{N,N},  {N,N},{N,0.},   {N,0.},   {N,N},{N,N},  {N,N}}
+};
+
 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;
+
+	/* special treatment for cosh(+/-inf + iy) if y is not a NaN */
+	if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
+		if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag) &&
+		    (z.imag != 0.)) {
+			if (z.real > 0) {
+				r.real = copysign(INF, cos(z.imag));
+				r.imag = copysign(INF, sin(z.imag));
+			}
+			else {
+				r.real = copysign(INF, cos(z.imag));
+				r.imag = -copysign(INF, sin(z.imag));
+			}
+		}
+		else {
+			r = cosh_special_values[special_type(z.real)]
+				               [special_type(z.imag)];
+		}
+		/* need to set errno = EDOM if y is +/- infinity and x is not
+		   a NaN */
+		if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real))
+			errno = EDOM;
+		else
+			errno = 0;
+		return r;
+	}
+
+	if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {
+		/* deal correctly with cases where cosh(z.real) overflows but
+		   cosh(z) does not. */
+		x_minus_one = z.real - copysign(1., z.real);
+		r.real = cos(z.imag) * cosh(x_minus_one) * Py_MATH_E;
+		r.imag = sin(z.imag) * sinh(x_minus_one) * Py_MATH_E;
+	} else {
+		r.real = cos(z.imag) * cosh(z.real);
+		r.imag = sin(z.imag) * sinh(z.real);
+	}
+	/* detect overflow, and set errno accordingly */
+	if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
+		errno = ERANGE;
+	else
+		errno = 0;
 	return r;
 }
 
@@ -139,13 +470,65 @@
 "Return the hyperbolic cosine of x.");
 
 
+/* exp(infinity + i*y) and exp(-infinity + i*y) need special treatment for
+   finite y */
+static Py_complex exp_special_values[7][7] = {
+  {{0.,0.},{U,U},{0.,-0.}, {0.,0.}, {U,U},{0.,0.},{0.,0.}},
+  {{N,N},  {U,U},{U,U},    {U,U},   {U,U},{N,N},  {N,N}},
+  {{N,N},  {U,U},{1.,-0.}, {1.,0.}, {U,U},{N,N},  {N,N}},
+  {{N,N},  {U,U},{1.,-0.}, {1.,0.}, {U,U},{N,N},  {N,N}},
+  {{N,N},  {U,U},{U,U},    {U,U},   {U,U},{N,N},  {N,N}},
+  {{INF,N},{U,U},{INF,-0.},{INF,0.},{U,U},{INF,N},{INF,N}},
+  {{N,N},  {N,N},{N,-0.},  {N,0.},  {N,N},{N,N},  {N,N}}
+};
+
 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 (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
+		if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)
+		    && (z.imag != 0.)) {
+			if (z.real > 0) {
+				r.real = copysign(INF, cos(z.imag));
+				r.imag = copysign(INF, sin(z.imag));
+			}
+			else {
+				r.real = copysign(0., cos(z.imag));
+				r.imag = copysign(0., sin(z.imag));
+			}
+		}
+		else {
+			r = exp_special_values[special_type(z.real)]
+				              [special_type(z.imag)];
+		}
+		/* need to set errno = EDOM if y is +/- infinity and x is not
+		   a NaN and not -infinity */
+		if (Py_IS_INFINITY(z.imag) &&
+		    (Py_IS_FINITE(z.real) ||
+		     (Py_IS_INFINITY(z.real) && z.real > 0)))
+			errno = EDOM;
+		else
+			errno = 0;
+		return r;
+	}
+
+	if (z.real > CM_LOG_LARGE_DOUBLE) {
+		l = exp(z.real-1.);
+		r.real = l*cos(z.imag)*Py_MATH_E;
+		r.imag = l*sin(z.imag)*Py_MATH_E;
+	} else {
+		l = exp(z.real);
+		r.real = l*cos(z.imag);
+		r.imag = l*sin(z.imag);
+	}
+	/* detect overflow, and set errno accordingly */
+	if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
+		errno = ERANGE;
+	else
+		errno = 0;
 	return r;
 }
 
@@ -155,24 +538,97 @@
 "Return the exponential value e**x.");
 
 
+static Py_complex log_special_values[7][7] = {
+  {{INF,-P34},{INF,-P}, {INF,-P},  {INF,P},  {INF,P}, {INF,P34}, {INF,N}},
+  {{INF,-P12},{U,U},    {U,U},     {U,U},    {U,U},   {INF,P12}, {N,N}},
+  {{INF,-P12},{U,U},    {-INF,-P}, {-INF,P}, {U,U},   {INF,P12}, {N,N}},
+  {{INF,-P12},{U,U},    {-INF,-0.},{-INF,0.},{U,U},   {INF,P12}, {N,N}},
+  {{INF,-P12},{U,U},    {U,U},     {U,U},    {U,U},   {INF,P12}, {N,N}},
+  {{INF,-P14},{INF,-0.},{INF,-0.}, {INF,0.}, {INF,0.},{INF,P14}, {INF,N}},
+  {{INF,N},   {N,N},    {N,N},     {N,N},    {N,N},   {INF,N},   {N,N}}
+};
+
 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 absolute error
+	   in the log.  This can be solved by rescaling z by a suitable power
+	   of 2.
+
+	   (2) the absolute value of z is greater than DBL_MAX (e.g. when both
+	   z.real and z.imag are within a factor of 1/sqrt(2) of DBL_MAX)
+	   Again, rescaling solves this.
+
+	   (3) the absolute value of z is close to 1.  In this case it's
+	   difficult to achieve good accuracy, at least in part because 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.
+
+	   (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.
+
+	 */
+
 	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, an, h;
+
+	SPECIAL_VALUE(z, log_special_values);
+
+	ax = fabs(z.real);
+	ay = fabs(z.imag);
+
+	if (ax > CM_LARGE_DOUBLE || ay > CM_LARGE_DOUBLE) {
+		r.real = log(hypot(ax/2., ay/2.)) + M_LN2;
+	} else if (ax < DBL_MIN && ay < DBL_MIN) {
+		if (ax > 0. || ay > 0.) {
+			/* catch cases where hypot(ax, ay) is subnormal */
+			r.real = log(hypot(ldexp(ax, DBL_MANT_DIG),
+				 ldexp(ay, DBL_MANT_DIG))) - DBL_MANT_DIG*M_LN2;
+		}
+		else {
+			/* log(+/-0. +/- 0i) */
+			r.real = -INF;
+			r.imag = atan2(z.imag, z.real);
+			errno = EDOM;
+			return r;
+		}
+	} else {
+		h = hypot(ax, ay);
+		if (0.71 <= h && h <= 1.73) {
+			am = ax > ay ? ax : ay;  /* max(ax, ay) */
+			an = ax > ay ? ay : ax;  /* min(ax, ay) */
+			r.real = log1p((am-1)*(am+1)+an*an)/2.;
+		} else {
+			r.real = log(h);
+		}
+	}
+	r.imag = atan2(z.imag, z.real);
+	errno = 0;
 	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);
+	int errno_save;
+
+	r = c_log(z);
+	errno_save = errno; /* just in case the divisions affect errno */
+	r.real = r.real / M_LN10;
+	r.imag = r.imag / M_LN10;
+	errno = errno_save;
 	return r;
 }
 
@@ -182,23 +638,16 @@
 "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;
-	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);
+	/* 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;
 }
 
@@ -208,12 +657,63 @@
 "Return the sine of x.");
 
 
+/* sinh(infinity + i*y) needs to be dealt with specially */
+static Py_complex sinh_special_values[7][7] = {
+  {{INF,N},{U,U},{-INF,-0.},{-INF,0.},{U,U},{INF,N},{INF,N}},
+  {{N,N},  {U,U},{U,U},     {U,U},    {U,U},{N,N},  {N,N}},
+  {{0.,N}, {U,U},{-0.,-0.}, {-0.,0.}, {U,U},{0.,N}, {0.,N}},
+  {{0.,N}, {U,U},{0.,-0.},  {0.,0.},  {U,U},{0.,N}, {0.,N}},
+  {{N,N},  {U,U},{U,U},     {U,U},    {U,U},{N,N},  {N,N}},
+  {{INF,N},{U,U},{INF,-0.}, {INF,0.}, {U,U},{INF,N},{INF,N}},
+  {{N,N},  {N,N},{N,-0.},   {N,0.},   {N,N},{N,N},  {N,N}}
+};
+
 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;
+
+	/* special treatment for sinh(+/-inf + iy) if y is finite and
+	   nonzero */
+	if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
+		if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)
+		    && (z.imag != 0.)) {
+			if (z.real > 0) {
+				r.real = copysign(INF, cos(z.imag));
+				r.imag = copysign(INF, sin(z.imag));
+			}
+			else {
+				r.real = -copysign(INF, cos(z.imag));
+				r.imag = copysign(INF, sin(z.imag));
+			}
+		}
+		else {
+			r = sinh_special_values[special_type(z.real)]
+				               [special_type(z.imag)];
+		}
+		/* need to set errno = EDOM if y is +/- infinity and x is not
+		   a NaN */
+		if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real))
+			errno = EDOM;
+		else
+			errno = 0;
+		return r;
+	}
+
+	if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {
+		x_minus_one = z.real - copysign(1., z.real);
+		r.real = cos(z.imag) * sinh(x_minus_one) * Py_MATH_E;
+		r.imag = sin(z.imag) * cosh(x_minus_one) * Py_MATH_E;
+	} else {
+		r.real = cos(z.imag) * sinh(z.real);
+		r.imag = sin(z.imag) * cosh(z.real);
+	}
+	/* detect overflow, and set errno accordingly */
+	if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))
+		errno = ERANGE;
+	else
+		errno = 0;
 	return r;
 }
 
@@ -223,29 +723,80 @@
 "Return the hyperbolic sine of x.");
 
 
+static Py_complex sqrt_special_values[7][7] = {
+  {{INF,-INF},{0.,-INF},{0.,-INF},{0.,INF},{0.,INF},{INF,INF},{N,INF}},
+  {{INF,-INF},{U,U},    {U,U},    {U,U},   {U,U},   {INF,INF},{N,N}},
+  {{INF,-INF},{U,U},    {0.,-0.}, {0.,0.}, {U,U},   {INF,INF},{N,N}},
+  {{INF,-INF},{U,U},    {0.,-0.}, {0.,0.}, {U,U},   {INF,INF},{N,N}},
+  {{INF,-INF},{U,U},    {U,U},    {U,U},   {U,U},   {INF,INF},{N,N}},
+  {{INF,-INF},{INF,-0.},{INF,-0.},{INF,0.},{INF,0.},{INF,INF},{INF,N}},
+  {{INF,-INF},{N,N},    {N,N},    {N,N},   {N,N},   {INF,INF},{N,N}}
+};
+
 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/divisions, 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;
+
+	SPECIAL_VALUE(z, sqrt_special_values);
+
+	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 < DBL_MIN && ay < DBL_MIN && (ax > 0. || ay > 0.)) {
+		/* here we catch cases where hypot(ax, ay) is subnormal */
+		ax = ldexp(ax, CM_SCALE_UP);
+		s = ldexp(sqrt(ax + hypot(ax, ldexp(ay, CM_SCALE_UP))),
+			  CM_SCALE_DOWN);
+	} 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);
+	}
+	errno = 0;
 	return r;
 }
 
@@ -256,23 +807,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;
 }
 
@@ -282,24 +825,78 @@
 "Return the tangent of x.");
 
 
+/* tanh(infinity + i*y) needs to be dealt with specially */
+static Py_complex tanh_special_values[7][7] = {
+  {{-1.,0.},{U,U},{-1.,-0.},{-1.,0.},{U,U},{-1.,0.},{-1.,0.}},
+  {{N,N},   {U,U},{U,U},    {U,U},   {U,U},{N,N},   {N,N}},
+  {{N,N},   {U,U},{-0.,-0.},{-0.,0.},{U,U},{N,N},   {N,N}},
+  {{N,N},   {U,U},{0.,-0.}, {0.,0.}, {U,U},{N,N},   {N,N}},
+  {{N,N},   {U,U},{U,U},    {U,U},   {U,U},{N,N},   {N,N}},
+  {{1.,0.}, {U,U},{1.,-0.}, {1.,0.}, {U,U},{1.,0.}, {1.,0.}},
+  {{N,N},   {N,N},{N,-0.},  {N,0.},  {N,N},{N,N},   {N,N}}
+};
+
 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;
+
+	/* special treatment for tanh(+/-inf + iy) if y is finite and
+	   nonzero */
+	if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
+		if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)
+		    && (z.imag != 0.)) {
+			if (z.real > 0) {
+				r.real = 1.0;
+				r.imag = copysign(0.,
+						  2.*sin(z.imag)*cos(z.imag));
+			}
+			else {
+				r.real = -1.0;
+				r.imag = copysign(0.,
+						  2.*sin(z.imag)*cos(z.imag));
+			}
+		}
+		else {
+			r = tanh_special_values[special_type(z.real)]
+				               [special_type(z.imag)];
+		}
+		/* need to set errno = EDOM if z.imag is +/-infinity and
+		   z.real is finite */
+		if (Py_IS_INFINITY(z.imag) && Py_IS_FINITE(z.real))
+			errno = EDOM;
+		else
+			errno = 0;
+		return r;
+	}
+
+	/* danger of overflow in 2.*z.imag !*/
+	if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {
+		r.real = copysign(1., z.real);
+		r.imag = 4.*sin(z.imag)*cos(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;
+	}
+	errno = 0;
 	return r;
 }
 
@@ -308,6 +905,7 @@
 "\n"
 "Return the hyperbolic tangent of x.");
 
+
 static PyObject *
 cmath_log(PyObject *self, PyObject *args)
 {
@@ -325,7 +923,6 @@
 	PyFPE_END_PROTECT(x)
 	if (errno != 0)
 		return math_error();
-	Py_ADJUST_ERANGE2(x.real, x.imag);
 	return PyComplex_FromCComplex(x);
 }
 
@@ -351,18 +948,24 @@
 static PyObject *
 math_1(PyObject *args, Py_complex (*func)(Py_complex))
 {
-	Py_complex x;
+	Py_complex x,r ;
 	if (!PyArg_ParseTuple(args, "D", &x))
 		return NULL;
 	errno = 0;
-	PyFPE_START_PROTECT("complex function", return 0)
-	x = (*func)(x);
-	PyFPE_END_PROTECT(x)
-	Py_ADJUST_ERANGE2(x.real, x.imag);
-	if (errno != 0)
-		return math_error();
-	else
-		return PyComplex_FromCComplex(x);
+	PyFPE_START_PROTECT("complex function", return 0);
+	r = (*func)(x);
+	PyFPE_END_PROTECT(r);
+	if (errno == EDOM) {
+		PyErr_SetString(PyExc_ValueError, "math domain error");
+		return NULL;
+	}
+	else if (errno == ERANGE) {
+		PyErr_SetString(PyExc_OverflowError, "math range error");
+		return NULL;
+	}
+	else {
+		return PyComplex_FromCComplex(r);
+	}
 }
 
 #define FUNC1(stubname, func) \
@@ -386,6 +989,151 @@
 FUNC1(cmath_tan, c_tan)
 FUNC1(cmath_tanh, c_tanh)
 
+static PyObject *
+cmath_phase(PyObject *self, PyObject *args)
+{
+	Py_complex z;
+	double phi;
+	if (!PyArg_ParseTuple(args, "D:phase", &z))
+		return NULL;
+	errno = 0;
+	PyFPE_START_PROTECT("arg function", return 0)
+	phi = c_atan2(z);
+	PyFPE_END_PROTECT(r)
+	if (errno != 0)
+		return math_error();
+	else
+		return PyFloat_FromDouble(phi);
+}
+
+PyDoc_STRVAR(cmath_phase_doc,
+"phase(z) -> float\n\n\
+Return argument, also known as the phase angle, of a complex.");
+
+static PyObject *
+cmath_polar(PyObject *self, PyObject *args)
+{
+	Py_complex z;
+	double r, phi;
+	if (!PyArg_ParseTuple(args, "D:polar", &z))
+		return NULL;
+	PyFPE_START_PROTECT("polar function", return 0)
+	phi = c_atan2(z); /* should not cause any exception */
+	r = c_abs(z); /* sets errno to ERANGE on overflow;  otherwise 0 */
+	PyFPE_END_PROTECT(r)
+	if (errno != 0)
+		return math_error();
+	else
+		return Py_BuildValue("dd", r, phi);
+}
+
+PyDoc_STRVAR(cmath_polar_doc,
+"polar(z) -> r: float, phi: float\n\n\
+Convert a complex from rectangular coordinates to polar coordinates. r is\n\
+the distance from 0 and phi the phase angle.");
+
+/*
+  rect() isn't covered by the C99 standard, but it's not too hard to
+  figure out 'spirit of C99' rules for special value handing:
+
+    rect(x, t) should behave like exp(log(x) + it) for positive-signed x
+    rect(x, t) should behave like -exp(log(-x) + it) for negative-signed x
+    rect(nan, t) should behave like exp(nan + it), except that rect(nan, 0)
+      gives nan +- i0 with the sign of the imaginary part unspecified.
+
+*/
+
+static Py_complex rect_special_values[7][7] = {
+  {{INF,N},{U,U},{-INF,0.},{-INF,-0.},{U,U},{INF,N},{INF,N}},
+  {{N,N},  {U,U},{U,U},    {U,U},     {U,U},{N,N},  {N,N}},
+  {{0.,0.},{U,U},{-0.,0.}, {-0.,-0.}, {U,U},{0.,0.},{0.,0.}},
+  {{0.,0.},{U,U},{0.,-0.}, {0.,0.},   {U,U},{0.,0.},{0.,0.}},
+  {{N,N},  {U,U},{U,U},    {U,U},     {U,U},{N,N},  {N,N}},
+  {{INF,N},{U,U},{INF,-0.},{INF,0.},  {U,U},{INF,N},{INF,N}},
+  {{N,N},  {N,N},{N,0.},   {N,0.},    {N,N},{N,N},  {N,N}}
+};
+
+static PyObject *
+cmath_rect(PyObject *self, PyObject *args)
+{
+	Py_complex z;
+	double r, phi;
+	if (!PyArg_ParseTuple(args, "dd:rect", &r, &phi))
+		return NULL;
+	errno = 0;
+	PyFPE_START_PROTECT("rect function", return 0)
+
+	/* deal with special values */
+	if (!Py_IS_FINITE(r) || !Py_IS_FINITE(phi)) {
+		/* if r is +/-infinity and phi is finite but nonzero then
+		   result is (+-INF +-INF i), but we need to compute cos(phi)
+		   and sin(phi) to figure out the signs. */
+		if (Py_IS_INFINITY(r) && (Py_IS_FINITE(phi)
+					  && (phi != 0.))) {
+			if (r > 0) {
+				z.real = copysign(INF, cos(phi));
+				z.imag = copysign(INF, sin(phi));
+			}
+			else {
+				z.real = -copysign(INF, cos(phi));
+				z.imag = -copysign(INF, sin(phi));
+			}
+		}
+		else {
+			z = rect_special_values[special_type(r)]
+				               [special_type(phi)];
+		}
+		/* need to set errno = EDOM if r is a nonzero number and phi
+		   is infinite */
+		if (r != 0. && !Py_IS_NAN(r) && Py_IS_INFINITY(phi))
+			errno = EDOM;
+		else
+			errno = 0;
+	}
+	else {
+		z.real = r * cos(phi);
+		z.imag = r * sin(phi);
+		errno = 0;
+	}
+
+	PyFPE_END_PROTECT(z)
+	if (errno != 0)
+		return math_error();
+	else
+		return PyComplex_FromCComplex(z);
+}
+
+PyDoc_STRVAR(cmath_rect_doc,
+"rect(r, phi) -> z: complex\n\n\
+Convert from polar coordinates to rectangular coordinates.");
+
+static PyObject *
+cmath_isnan(PyObject *self, PyObject *args)
+{
+	Py_complex z;
+	if (!PyArg_ParseTuple(args, "D:isnan", &z))
+		return NULL;
+	return PyBool_FromLong(Py_IS_NAN(z.real) || Py_IS_NAN(z.imag));
+}
+
+PyDoc_STRVAR(cmath_isnan_doc,
+"isnan(z) -> bool\n\
+Checks if the real or imaginary part of z not a number (NaN)");
+
+static PyObject *
+cmath_isinf(PyObject *self, PyObject *args)
+{
+	Py_complex z;
+	if (!PyArg_ParseTuple(args, "D:isnan", &z))
+		return NULL;
+	return PyBool_FromLong(Py_IS_INFINITY(z.real) ||
+			       Py_IS_INFINITY(z.imag));
+}
+
+PyDoc_STRVAR(cmath_isinf_doc,
+"isinf(z) -> bool\n\
+Checks if the real or imaginary part of z is infinite.");
+
 
 PyDoc_STRVAR(module_doc,
 "This module is always available. It provides access to mathematical\n"
@@ -401,8 +1149,13 @@
 	{"cos",    cmath_cos,   METH_VARARGS, c_cos_doc},
 	{"cosh",   cmath_cosh,  METH_VARARGS, c_cosh_doc},
 	{"exp",    cmath_exp,   METH_VARARGS, c_exp_doc},
+	{"isinf",  cmath_isinf, METH_VARARGS, cmath_isinf_doc},
+	{"isnan",  cmath_isnan, METH_VARARGS, cmath_isnan_doc},
 	{"log",    cmath_log,   METH_VARARGS, cmath_log_doc},
 	{"log10",  cmath_log10, METH_VARARGS, c_log10_doc},
+	{"phase",  cmath_phase, METH_VARARGS, cmath_phase_doc},
+	{"polar",  cmath_polar, METH_VARARGS, cmath_polar_doc},
+	{"rect",   cmath_rect,  METH_VARARGS, cmath_rect_doc},
 	{"sin",    cmath_sin,   METH_VARARGS, c_sin_doc},
 	{"sinh",   cmath_sinh,  METH_VARARGS, c_sinh_doc},
 	{"sqrt",   cmath_sqrt,  METH_VARARGS, c_sqrt_doc},
@@ -421,6 +1174,6 @@
 		return;
 
 	PyModule_AddObject(m, "pi",
-                           PyFloat_FromDouble(atan(1.0) * 4.0));
-	PyModule_AddObject(m, "e", PyFloat_FromDouble(exp(1.0)));
+                           PyFloat_FromDouble(Py_MATH_PI));
+	PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
 }
diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c
index cf2bf64..8c48316 100644
--- a/Modules/mathmodule.c
+++ b/Modules/mathmodule.c
@@ -1,17 +1,60 @@
 /* Math module -- standard C math library functions, pi and e */
 
+/* Here are some comments from Tim Peters, extracted from the
+   discussion attached to http://bugs.python.org/issue1640.  They
+   describe the general aims of the math module with respect to
+   special values, IEEE-754 floating-point exceptions, and Python
+   exceptions.
+
+These are the "spirit of 754" rules:
+
+1. If the mathematical result is a real number, but of magnitude too
+large to approximate by a machine float, overflow is signaled and the
+result is an infinity (with the appropriate sign).
+
+2. If the mathematical result is a real number, but of magnitude too
+small to approximate by a machine float, underflow is signaled and the
+result is a zero (with the appropriate sign).
+
+3. At a singularity (a value x such that the limit of f(y) as y
+approaches x exists and is an infinity), "divide by zero" is signaled
+and the result is an infinity (with the appropriate sign).  This is
+complicated a little by that the left-side and right-side limits may
+not be the same; e.g., 1/x approaches +inf or -inf as x approaches 0
+from the positive or negative directions.  In that specific case, the
+sign of the zero determines the result of 1/0.
+
+4. At a point where a function has no defined result in the extended
+reals (i.e., the reals plus an infinity or two), invalid operation is
+signaled and a NaN is returned.
+
+And these are what Python has historically /tried/ to do (but not
+always successfully, as platform libm behavior varies a lot):
+
+For #1, raise OverflowError.
+
+For #2, return a zero (with the appropriate sign if that happens by
+accident ;-)).
+
+For #3 and #4, raise ValueError.  It may have made sense to raise
+Python's ZeroDivisionError in #3, but historically that's only been
+raised for division by zero and mod by zero.
+
+*/
+
+/*
+   In general, on an IEEE-754 platform the aim is to follow the C99
+   standard, including Annex 'F', whenever possible.  Where the
+   standard recommends raising the 'divide-by-zero' or 'invalid'
+   floating-point exceptions, Python should raise a ValueError.  Where
+   the standard recommends raising 'overflow', Python should raise an
+   OverflowError.  In all other circumstances a value should be
+   returned.
+ */
+
 #include "Python.h"
 #include "longintrepr.h" /* just for SHIFT */
 
-#ifndef _MSC_VER
-#ifndef __STDC__
-extern double fmod (double, double);
-extern double frexp (double, int *);
-extern double ldexp (double, int);
-extern double modf (double, double *);
-#endif /* __STDC__ */
-#endif /* _MSC_VER */
-
 #ifdef _OSF_SOURCE
 /* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */
 extern double copysign(double, double);
@@ -52,41 +95,111 @@
 	return result;
 }
 
+/*
+   math_1 is used to wrap a libm function f that takes a double
+   arguments and returns a double.
+
+   The error reporting follows these rules, which are designed to do
+   the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
+   platforms.
+
+   - a NaN result from non-NaN inputs causes ValueError to be raised
+   - an infinite result from finite inputs causes OverflowError to be
+     raised if can_overflow is 1, or raises ValueError if can_overflow
+     is 0.
+   - if the result is finite and errno == EDOM then ValueError is
+     raised
+   - if the result is finite and nonzero and errno == ERANGE then
+     OverflowError is raised
+
+   The last rule is used to catch overflow on platforms which follow
+   C89 but for which HUGE_VAL is not an infinity.
+
+   For the majority of one-argument functions these rules are enough
+   to ensure that Python's functions behave as specified in 'Annex F'
+   of the C99 standard, with the 'invalid' and 'divide-by-zero'
+   floating-point exceptions mapping to Python's ValueError and the
+   'overflow' floating-point exception mapping to OverflowError.
+   math_1 only works for functions that don't have singularities *and*
+   the possibility of overflow; fortunately, that covers everything we
+   care about right now.
+*/
+
 static PyObject *
 math_1_to_whatever(PyObject *arg, double (*func) (double),
-                   PyObject *(*from_double_func) (double))
+                   PyObject *(*from_double_func) (double),
+                   int can_overflow)
 {
-	double x = PyFloat_AsDouble(arg);
+	double x, r;
+	x = PyFloat_AsDouble(arg);
 	if (x == -1.0 && PyErr_Occurred())
 		return NULL;
 	errno = 0;
-	PyFPE_START_PROTECT("in math_1", return 0)
-	x = (*func)(x);
-	PyFPE_END_PROTECT(x)
-	Py_SET_ERRNO_ON_MATH_ERROR(x);
-	if (errno && is_error(x))
+	PyFPE_START_PROTECT("in math_1", return 0);
+	r = (*func)(x);
+	PyFPE_END_PROTECT(r);
+	if (Py_IS_NAN(r)) {
+		if (!Py_IS_NAN(x))
+			errno = EDOM;
+		else
+			errno = 0;
+	}
+	else if (Py_IS_INFINITY(r)) {
+		if (Py_IS_FINITE(x))
+			errno = can_overflow ? ERANGE : EDOM;
+		else
+			errno = 0;
+	}
+	if (errno && is_error(r))
 		return NULL;
 	else
-        	return (*from_double_func)(x);
+        	return (*from_double_func)(r);
+}
+
+/*
+   math_2 is used to wrap a libm function f that takes two double
+   arguments and returns a double.
+
+   The error reporting follows these rules, which are designed to do
+   the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
+   platforms.
+
+   - a NaN result from non-NaN inputs causes ValueError to be raised
+   - an infinite result from finite inputs causes OverflowError to be
+     raised.
+   - if the result is finite and errno == EDOM then ValueError is
+     raised
+   - if the result is finite and nonzero and errno == ERANGE then
+     OverflowError is raised
+
+   The last rule is used to catch overflow on platforms which follow
+   C89 but for which HUGE_VAL is not an infinity.
+
+   For most two-argument functions (copysign, fmod, hypot, atan2)
+   these rules are enough to ensure that Python's functions behave as
+   specified in 'Annex F' of the C99 standard, with the 'invalid' and
+   'divide-by-zero' floating-point exceptions mapping to Python's
+   ValueError and the 'overflow' floating-point exception mapping to
+   OverflowError.
+*/
+
+static PyObject *
+math_1(PyObject *arg, double (*func) (double), int can_overflow)
+{
+	return math_1_to_whatever(arg, func, PyFloat_FromDouble, can_overflow);
 }
 
 static PyObject *
-math_1(PyObject *arg, double (*func) (double))
+math_1_to_int(PyObject *arg, double (*func) (double), int can_overflow)
 {
-	return math_1_to_whatever(arg, func, PyFloat_FromDouble);
-}
-
-static PyObject *
-math_1_to_int(PyObject *arg, double (*func) (double))
-{
-	return math_1_to_whatever(arg, func, PyLong_FromDouble);
+	return math_1_to_whatever(arg, func, PyLong_FromDouble, can_overflow);
 }
 
 static PyObject *
 math_2(PyObject *args, double (*func) (double, double), char *funcname)
 {
 	PyObject *ox, *oy;
-	double x, y;
+	double x, y, r;
 	if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy))
 		return NULL;
 	x = PyFloat_AsDouble(ox);
@@ -94,19 +207,30 @@
 	if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
 		return NULL;
 	errno = 0;
-	PyFPE_START_PROTECT("in math_2", return 0)
-	x = (*func)(x, y);
-	PyFPE_END_PROTECT(x)
-	Py_SET_ERRNO_ON_MATH_ERROR(x);
-	if (errno && is_error(x))
+	PyFPE_START_PROTECT("in math_2", return 0);
+	r = (*func)(x, y);
+	PyFPE_END_PROTECT(r);
+	if (Py_IS_NAN(r)) {
+		if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
+			errno = EDOM;
+		else
+			errno = 0;
+	}
+	else if (Py_IS_INFINITY(r)) {
+		if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
+			errno = ERANGE;
+		else
+			errno = 0;
+	}
+	if (errno && is_error(r))
 		return NULL;
 	else
-		return PyFloat_FromDouble(x);
+		return PyFloat_FromDouble(r);
 }
 
-#define FUNC1(funcname, func, docstring) \
+#define FUNC1(funcname, func, can_overflow, docstring)			\
 	static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
-		return math_1(args, func); \
+		return math_1(args, func, can_overflow);		    \
 	}\
         PyDoc_STRVAR(math_##funcname##_doc, docstring);
 
@@ -116,15 +240,21 @@
 	}\
         PyDoc_STRVAR(math_##funcname##_doc, docstring);
 
-FUNC1(acos, acos,
+FUNC1(acos, acos, 0,
       "acos(x)\n\nReturn the arc cosine (measured in radians) of x.")
-FUNC1(asin, asin,
+FUNC1(acosh, acosh, 0,
+      "acosh(x)\n\nReturn the hyperbolic arc cosine (measured in radians) of x.")
+FUNC1(asin, asin, 0,
       "asin(x)\n\nReturn the arc sine (measured in radians) of x.")
-FUNC1(atan, atan,
+FUNC1(asinh, asinh, 0,
+      "asinh(x)\n\nReturn the hyperbolic arc sine (measured in radians) of x.")
+FUNC1(atan, atan, 0,
       "atan(x)\n\nReturn the arc tangent (measured in radians) of x.")
 FUNC2(atan2, atan2,
       "atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n"
       "Unlike atan(y/x), the signs of both x and y are considered.")
+FUNC1(atanh, atanh, 0,
+      "atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.")
 
 static PyObject * math_ceil(PyObject *self, PyObject *number) {
 	static PyObject *ceil_str = NULL;
@@ -138,7 +268,7 @@
 
 	method = _PyType_Lookup(Py_TYPE(number), ceil_str);
 	if (method == NULL)
-		return math_1_to_int(number, ceil);
+		return math_1_to_int(number, ceil, 0);
 	else
 		return PyObject_CallFunction(method, "O", number);
 }
@@ -147,23 +277,15 @@
 	     "ceil(x)\n\nReturn the ceiling of x as an int.\n"
 	     "This is the smallest integral value >= x.");
 
-FUNC1(cos, cos,
-      "cos(x)\n\nReturn the cosine of x (measured in radians).")
-FUNC1(cosh, cosh,
-      "cosh(x)\n\nReturn the hyperbolic cosine of x.")
-
-#ifdef MS_WINDOWS
-#  define copysign _copysign
-#  define HAVE_COPYSIGN 1
-#endif
-#ifdef HAVE_COPYSIGN
 FUNC2(copysign, copysign,
-      "copysign(x,y)\n\nReturn x with the sign of y.");
-#endif
-
-FUNC1(exp, exp,
+      "copysign(x,y)\n\nReturn x with the sign of y.")
+FUNC1(cos, cos, 0,
+      "cos(x)\n\nReturn the cosine of x (measured in radians).")
+FUNC1(cosh, cosh, 1,
+      "cosh(x)\n\nReturn the hyperbolic cosine of x.")
+FUNC1(exp, exp, 1,
       "exp(x)\n\nReturn e raised to the power of x.")
-FUNC1(fabs, fabs,
+FUNC1(fabs, fabs, 0,
       "fabs(x)\n\nReturn the absolute value of the float x.")
 
 static PyObject * math_floor(PyObject *self, PyObject *number) {
@@ -178,7 +300,7 @@
 
 	method = _PyType_Lookup(Py_TYPE(number), floor_str);
 	if (method == NULL)
-        	return math_1_to_int(number, floor);
+        	return math_1_to_int(number, floor, 0);
 	else
 		return PyObject_CallFunction(method, "O", number);
 }
@@ -187,22 +309,18 @@
 	     "floor(x)\n\nReturn the floor of x as an int.\n"
 	     "This is the largest integral value <= x.");
 
-FUNC2(fmod, fmod,
-      "fmod(x,y)\n\nReturn fmod(x, y), according to platform C."
-      "  x % y may differ.")
-FUNC2(hypot, hypot,
-      "hypot(x,y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).")
-FUNC2(pow, pow,
-      "pow(x,y)\n\nReturn x**y (x to the power of y).")
-FUNC1(sin, sin,
+FUNC1(log1p, log1p, 1,
+      "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n\
+      The result is computed in a way which is accurate for x near zero.")
+FUNC1(sin, sin, 0,
       "sin(x)\n\nReturn the sine of x (measured in radians).")
-FUNC1(sinh, sinh,
+FUNC1(sinh, sinh, 1,
       "sinh(x)\n\nReturn the hyperbolic sine of x.")
-FUNC1(sqrt, sqrt,
+FUNC1(sqrt, sqrt, 0,
       "sqrt(x)\n\nReturn the square root of x.")
-FUNC1(tan, tan,
+FUNC1(tan, tan, 0,
       "tan(x)\n\nReturn the tangent of x (measured in radians).")
-FUNC1(tanh, tanh,
+FUNC1(tanh, tanh, 0,
       "tanh(x)\n\nReturn the hyperbolic tangent of x.")
 
 static PyObject *
@@ -244,13 +362,17 @@
 	double x = PyFloat_AsDouble(arg);
 	if (x == -1.0 && PyErr_Occurred())
 		return NULL;
-	errno = 0;
-	x = frexp(x, &i);
-	Py_SET_ERRNO_ON_MATH_ERROR(x);
-	if (errno && is_error(x))
-		return NULL;
-	else
-		return Py_BuildValue("(di)", x, i);
+	/* deal with special cases directly, to sidestep platform
+	   differences */
+	if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) {
+		i = 0;
+	}
+	else {
+		PyFPE_START_PROTECT("in math_frexp", return 0);
+		x = frexp(x, &i);
+		PyFPE_END_PROTECT(x);
+	}
+	return Py_BuildValue("(di)", x, i);
 }
 
 PyDoc_STRVAR(math_frexp_doc,
@@ -263,19 +385,24 @@
 static PyObject *
 math_ldexp(PyObject *self, PyObject *args)
 {
-	double x;
+	double x, r;
 	int exp;
 	if (! PyArg_ParseTuple(args, "di:ldexp", &x, &exp))
 		return NULL;
 	errno = 0;
-	PyFPE_START_PROTECT("ldexp", return 0)
-	x = ldexp(x, exp);
-	PyFPE_END_PROTECT(x)
-	Py_SET_ERRNO_ON_MATH_ERROR(x);
-	if (errno && is_error(x))
+	PyFPE_START_PROTECT("in math_ldexp", return 0)
+	r = ldexp(x, exp);
+	PyFPE_END_PROTECT(r)
+	if (Py_IS_FINITE(x) && Py_IS_INFINITY(r))
+		errno = ERANGE;
+	/* Windows MSVC8 sets errno = EDOM on ldexp(NaN, i);
+	   we unset it to avoid raising a ValueError here. */
+	if (errno == EDOM)
+		errno = 0;
+	if (errno && is_error(r))
 		return NULL;
 	else
-		return PyFloat_FromDouble(x);
+		return PyFloat_FromDouble(r);
 }
 
 PyDoc_STRVAR(math_ldexp_doc,
@@ -288,12 +415,10 @@
 	if (x == -1.0 && PyErr_Occurred())
 		return NULL;
 	errno = 0;
+	PyFPE_START_PROTECT("in math_modf", return 0);
 	x = modf(x, &y);
-	Py_SET_ERRNO_ON_MATH_ERROR(x);
-	if (errno && is_error(x))
-		return NULL;
-	else
-		return Py_BuildValue("(dd)", x, y);
+	PyFPE_END_PROTECT(x);
+	return Py_BuildValue("(dd)", x, y);
 }
 
 PyDoc_STRVAR(math_modf_doc,
@@ -332,7 +457,7 @@
 	}
 
 	/* Else let libm handle it by itself. */
-	return math_1(arg, func);
+	return math_1(arg, func, 0);
 }
 
 static PyObject *
@@ -375,6 +500,141 @@
 PyDoc_STRVAR(math_log10_doc,
 "log10(x) -> the base 10 logarithm of x.");
 
+static PyObject *
+math_fmod(PyObject *self, PyObject *args)
+{
+	PyObject *ox, *oy;
+	double r, x, y;
+	if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy))
+		return NULL;
+	x = PyFloat_AsDouble(ox);
+	y = PyFloat_AsDouble(oy);
+	if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
+		return NULL;
+	/* fmod(x, +/-Inf) returns x for finite x. */
+	if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))
+		return PyFloat_FromDouble(x);
+	errno = 0;
+	PyFPE_START_PROTECT("in math_fmod", return 0);
+	r = fmod(x, y);
+	PyFPE_END_PROTECT(r);
+	if (Py_IS_NAN(r)) {
+		if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
+			errno = EDOM;
+		else
+			errno = 0;
+	}
+	if (errno && is_error(r))
+		return NULL;
+	else
+		return PyFloat_FromDouble(r);
+}
+
+PyDoc_STRVAR(math_fmod_doc,
+"fmod(x,y)\n\nReturn fmod(x, y), according to platform C."
+"  x % y may differ.");
+
+static PyObject *
+math_hypot(PyObject *self, PyObject *args)
+{
+	PyObject *ox, *oy;
+	double r, x, y;
+	if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy))
+		return NULL;
+	x = PyFloat_AsDouble(ox);
+	y = PyFloat_AsDouble(oy);
+	if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
+		return NULL;
+	/* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */
+	if (Py_IS_INFINITY(x))
+		return PyFloat_FromDouble(fabs(x));
+	if (Py_IS_INFINITY(y))
+		return PyFloat_FromDouble(fabs(y));
+	errno = 0;
+	PyFPE_START_PROTECT("in math_hypot", return 0);
+	r = hypot(x, y);
+	PyFPE_END_PROTECT(r);
+	if (Py_IS_NAN(r)) {
+		if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
+			errno = EDOM;
+		else
+			errno = 0;
+	}
+	else if (Py_IS_INFINITY(r)) {
+		if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
+			errno = ERANGE;
+		else
+			errno = 0;
+	}
+	if (errno && is_error(r))
+		return NULL;
+	else
+		return PyFloat_FromDouble(r);
+}
+
+PyDoc_STRVAR(math_hypot_doc,
+"hypot(x,y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).");
+
+/* pow can't use math_2, but needs its own wrapper: the problem is
+   that an infinite result can arise either as a result of overflow
+   (in which case OverflowError should be raised) or as a result of
+   e.g. 0.**-5. (for which ValueError needs to be raised.)
+*/
+
+static PyObject *
+math_pow(PyObject *self, PyObject *args)
+{
+	PyObject *ox, *oy;
+	double r, x, y;
+
+	if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy))
+		return NULL;
+	x = PyFloat_AsDouble(ox);
+	y = PyFloat_AsDouble(oy);
+	if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
+		return NULL;
+	/* 1**x and x**0 return 1., even if x is a NaN or infinity. */
+	if (x == 1.0 || y == 0.0)
+	        return PyFloat_FromDouble(1.);
+	errno = 0;
+	PyFPE_START_PROTECT("in math_pow", return 0);
+	r = pow(x, y);
+	PyFPE_END_PROTECT(r);
+	if (Py_IS_NAN(r)) {
+		if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
+			errno = EDOM;
+		else
+			errno = 0;
+	}
+	/* an infinite result arises either from:
+
+	   (A) (+/-0.)**negative,
+	   (B) overflow of x**y with both x and y finite (and x nonzero)
+	   (C) (+/-inf)**positive, or
+	   (D) x**inf with |x| > 1, or x**-inf with |x| < 1.
+
+	   In case (A) we want ValueError to be raised.  In case (B)
+	   OverflowError should be raised.  In cases (C) and (D) the infinite
+	   result should be returned.
+	*/
+	else if (Py_IS_INFINITY(r)) {
+		if (x == 0.)
+			errno = EDOM;
+		else if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
+			errno = ERANGE;
+		else
+			errno = 0;
+	}
+
+	if (errno && is_error(r))
+		return NULL;
+	else
+		return PyFloat_FromDouble(r);
+}
+
+PyDoc_STRVAR(math_pow_doc,
+"pow(x,y)\n\nReturn x**y (x to the power of y).");
+
 static const double degToRad = Py_MATH_PI / 180.0;
 static const double radToDeg = 180.0 / Py_MATH_PI;
 
@@ -428,16 +688,16 @@
 "isinf(x) -> bool\n\
 Checks if float x is infinite (positive or negative)");
 
-
 static PyMethodDef math_methods[] = {
 	{"acos",	math_acos,	METH_O,		math_acos_doc},
+	{"acosh",	math_acosh,	METH_O,		math_acosh_doc},
 	{"asin",	math_asin,	METH_O,		math_asin_doc},
+	{"asinh",	math_asinh,	METH_O,		math_asinh_doc},
 	{"atan",	math_atan,	METH_O,		math_atan_doc},
 	{"atan2",	math_atan2,	METH_VARARGS,	math_atan2_doc},
+	{"atanh",	math_atanh,	METH_O,		math_atanh_doc},
 	{"ceil",	math_ceil,	METH_O,		math_ceil_doc},
-#ifdef HAVE_COPYSIGN
 	{"copysign",	math_copysign,	METH_VARARGS,	math_copysign_doc},
-#endif
 	{"cos",		math_cos,	METH_O,		math_cos_doc},
 	{"cosh",	math_cosh,	METH_O,		math_cosh_doc},
 	{"degrees",	math_degrees,	METH_O,		math_degrees_doc},
@@ -451,6 +711,7 @@
 	{"isnan",	math_isnan,	METH_O,		math_isnan_doc},
 	{"ldexp",	math_ldexp,	METH_VARARGS,	math_ldexp_doc},
 	{"log",		math_log,	METH_VARARGS,	math_log_doc},
+	{"log1p",	math_log1p,	METH_O,		math_log1p_doc},
 	{"log10",	math_log10,	METH_O,		math_log10_doc},
 	{"modf",	math_modf,	METH_O,		math_modf_doc},
 	{"pow",		math_pow,	METH_VARARGS,	math_pow_doc},
@@ -472,27 +733,15 @@
 PyMODINIT_FUNC
 initmath(void)
 {
-	PyObject *m, *d, *v;
+	PyObject *m;
 
 	m = Py_InitModule3("math", math_methods, module_doc);
 	if (m == NULL)
 		goto finally;
-	d = PyModule_GetDict(m);
-	if (d == NULL)
-		goto finally;
 
-        if (!(v = PyFloat_FromDouble(Py_MATH_PI)))
-                goto finally;
-	if (PyDict_SetItemString(d, "pi", v) < 0)
-                goto finally;
-	Py_DECREF(v);
+	PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));
+	PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
 
-        if (!(v = PyFloat_FromDouble(Py_MATH_E)))
-                goto finally;
-	if (PyDict_SetItemString(d, "e", v) < 0)
-                goto finally;
-	Py_DECREF(v);
-
-  finally:
+    finally:
 	return;
 }