I finally got the time to update and merge Mark's and my trunk-math b…

…ranch. The patch is collaborated work of Mark Dickinson and me. It was mostly done a few months ago. The patch fixes a lot of loose ends and edge cases related to operations with NaN, INF, very small values and complex math.

The patch also adds acosh, asinh, atanh, log1p and copysign to all platforms. Finally it fixes differences between platforms like different results or exceptions for edge cases. Have fun :)

Doc/library/cmath.rst

@@ -14,8 +14,81 @@ method: these methods are used to convert the object to a complex or
 floating-point number, respectively, and the function is then applied to the
 result of the conversion.
 
-The functions are:
+.. 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.
+
+
+Complex coordinates
+-------------------
+
+Complex numbers can be expressed by two important coordinate systems.
+Python's :class:`complex` type uses rectangular coordinates where a number
+on the complex plain is defined by two floats, the real part and the imaginary
+part.
+
+Definition::
+
+   z = x + 1j * y
+
+   x := real(z)
+   y := imag(z)
+
+In engineering the polar coordinate system is popular for complex numbers. In
+polar coordinates a complex number is defined by the radius *r* and the phase
+angle *φ*. The radius *r* is the absolute value of the complex, which can be
+viewed as distance from (0, 0). The radius *r* is always 0 or a positive float.
+The phase angle *φ* is the counter clockwise angle from the positive x axis,
+e.g. *1* has the angle *0*, *1j* has the angle *π/2* and *-1* the angle *-π*.
+
+.. note::
+   While :func:`phase` and func:`polar` return *+π* for a negative real they
+   may return *-π* for a complex with a very small negative imaginary
+   part, e.g. *-1-1E-300j*.
+
+
+Definition::
+
+   z = r * exp(1j * φ)
+   z = r * cis(φ)
+
+   r := abs(z) := sqrt(real(z)**2 + imag(z)**2)
+   phi := phase(z) := atan2(imag(z), real(z))
+   cis(φ) := cos(φ) + 1j * sin(φ)
+
+
+.. function:: phase(x)
+
+   Return phase, also known as the argument, of a complex.
+
+   .. versionadded:: 2.6
+
+
+.. function:: polar(x)
+
+   Convert a :class:`complex` from rectangular coordinates to polar 
+   coordinates. The function returns a tuple with the two elements
+   *r* and *phi*. *r* is the distance from 0 and *phi* the phase 
+   angle.
+
+   .. versionadded:: 2.6
+
+
+.. function:: rect(r, phi)
+
+   Convert from polar coordinates to rectangular coordinates and return
+   a :class:`complex`.
+
+   .. versionadded:: 2.6
+
+
+
+cmath functions
+---------------
 
 .. function:: acos(x)
 
@@ -37,30 +110,35 @@ The functions are:
 
 .. 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)
@@ -78,6 +156,21 @@ The functions are:
    Return the exponential value ``e**x``.
 
 
+.. function:: isinf(x)
+
+   Return *True* if the real or the imaginary part of x is positive
+   or negative infinity.
+
+   .. versionadded:: 2.6
+
+
+.. function:: isnan(x)
+
+   Return *True* if the real or imaginary part of x is not a number (NaN).
+
+   .. versionadded:: 2.6
+
+
 .. function:: log(x[, base])
 
    Returns the logarithm of *x* to the given *base*. If the *base* is not
@@ -154,3 +247,4 @@ cuts for numerical purposes, a good reference should be the following:
    nothing's sign bit.  In Iserles, A., and Powell, M. (eds.), The state of the art
    in numerical analysis. Clarendon Press (1987) pp165-211.
 
+

Doc/library/math.rst

@@ -139,14 +139,26 @@ Power and logarithmic functions:
       *base* argument added.
 
 
+.. function:: log1p(x)
+
+   Return the natural logarithm of *1+x* (base *e*). The
+   result is calculated in a way which is accurate for *x* near zero.
+
+   .. versionadded:: 2.6
+
+
 .. function:: log10(x)
 
    Return the base-10 logarithm of *x*.
 
 
 .. function:: pow(x, y)
 
-   Return ``x**y``.
+   Return ``x**y``. ``1.0**y`` returns *1.0*, even for ``1.0**nan``. ``0**y``
+   returns *0.* for all positive *y*, *0* and *NAN*.
+
+   .. versionchanged:: 2.6
+      The outcome of ``1**nan`` and ``0**nan`` was undefined.
 
 
 .. function:: sqrt(x)
@@ -197,6 +209,13 @@ Trigonometric functions:
    Return the sine of *x* radians.
 
 
+.. function:: asinh(x)
+
+   Return the inverse hyperbolic sine of *x*, in radians.
+
+   .. versionadded:: 2.6
+
+
 .. function:: tan(x)
 
    Return the tangent of *x* radians.
@@ -221,6 +240,13 @@ Hyperbolic functions:
    Return the hyperbolic cosine of *x*.
 
 
+.. function:: acosh(x)
+
+   Return the inverse hyperbolic cosine of *x*, in radians.
+
+   .. versionadded:: 2.6
+
+
 .. function:: sinh(x)
 
    Return the hyperbolic sine of *x*.
@@ -230,6 +256,14 @@ Hyperbolic functions:
 
    Return the hyperbolic tangent of *x*.
 
+
+.. function:: atanh(x)
+
+   Return the inverse hyperbolic tangent of *x*, in radians.
+
+   .. versionadded:: 2.6
+
+
 The module also defines two mathematical constants:
 
 
@@ -242,6 +276,7 @@ The module also defines two mathematical constants:
 
    The mathematical constant *e*.
 
+
 .. note::
 
    The :mod:`math` module consists mostly of thin wrappers around the platform C
@@ -255,9 +290,17 @@ The module also defines two mathematical constants:
    :exc:`OverflowError` isn't defined, and in cases where ``math.log(0)`` raises
    :exc:`OverflowError`, ``math.log(0L)`` may raise :exc:`ValueError` instead.
 
+   All functions return a quite *NaN* if at least one of the args is *NaN*.
+   Signaling *NaN*s raise an exception. The exception type still depends on the
+   platform and libm implementation. It's usually :exc:`ValueError` for *EDOM*
+   and :exc:`OverflowError` for errno *ERANGE*.
+
+   ..versionchanged:: 2.6
+      In earlier versions of Python the outcome of an operation with NaN as
+      input depended on platform and libm implementation.
+
 
 .. seealso::
 
    Module :mod:`cmath`
       Complex number versions of many of these functions.
-

Include/Python.h

@@ -73,6 +73,7 @@
 #if defined(PYMALLOC_DEBUG) && !defined(WITH_PYMALLOC)
 #error "PYMALLOC_DEBUG requires WITH_PYMALLOC"
 #endif
+#include "pymath.h"
 #include "pymem.h"
 
 #include "object.h"

Include/complexobject.h

@@ -19,13 +19,15 @@ typedef struct {
 #define c_prod _Py_c_prod
 #define c_quot _Py_c_quot
 #define c_pow _Py_c_pow
+#define c_abs _Py_c_abs
 
 PyAPI_FUNC(Py_complex) c_sum(Py_complex, Py_complex);
 PyAPI_FUNC(Py_complex) c_diff(Py_complex, Py_complex);
 PyAPI_FUNC(Py_complex) c_neg(Py_complex);
 PyAPI_FUNC(Py_complex) c_prod(Py_complex, Py_complex);
 PyAPI_FUNC(Py_complex) c_quot(Py_complex, Py_complex);
 PyAPI_FUNC(Py_complex) c_pow(Py_complex, Py_complex);
+PyAPI_FUNC(double) c_abs(Py_complex);
 
 
 /* Complex object interface */

Include/floatobject.h

@@ -21,6 +21,17 @@ PyAPI_DATA(PyTypeObject) PyFloat_Type;
 #define PyFloat_Check(op) PyObject_TypeCheck(op, &PyFloat_Type)
 #define PyFloat_CheckExact(op) (Py_TYPE(op) == &PyFloat_Type)
 
+#ifdef Py_NAN
+#define Py_RETURN_NAN PyFloat_FromDouble(Py_NAN)
+#endif
+
+#define Py_RETURN_INF(sign) do					\
+	if (copysign(1., sign) == 1.) {				\
+		return PyFloat_FromDouble(Py_HUGE_VAL);	\
+	} else {						\
+		return PyFloat_FromDouble(-Py_HUGE_VAL);	\
+	} while(0)
+
 PyAPI_FUNC(double) PyFloat_GetMax(void);
 PyAPI_FUNC(double) PyFloat_GetMin(void);
 PyAPI_FUNC(PyObject *) PyFloat_GetInfo(void);