Merged revisions 62380,62382-62383 via svnmerge from

svn+ssh://pythondev@svn.python.org/python/trunk

........
  r62380 | christian.heimes | 2008-04-19 01:13:07 +0200 (Sat, 19 Apr 2008) | 3 lines

  I finally got the time to update and merge Mark's and my trunk-math branch. 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 :)
........
  r62382 | christian.heimes | 2008-04-19 01:40:40 +0200 (Sat, 19 Apr 2008) | 2 lines

  Added new files to Windows project files
  More Windows related fixes are coming soon
........
  r62383 | christian.heimes | 2008-04-19 01:49:11 +0200 (Sat, 19 Apr 2008) | 1 line

  Stupid me. Py_RETURN_NAN should actually return something ...
........

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
@@ -151,3 +244,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

@@ -128,6 +128,14 @@ Power and logarithmic functions:
    return the natural logarithm of *x* (that is, the logarithm to base *e*).
 
 
+.. 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*.
@@ -135,7 +143,11 @@ Power and logarithmic functions:
 
 .. 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)
@@ -186,6 +198,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.
@@ -210,6 +229,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*.
@@ -219,6 +245,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:
 
 
@@ -231,6 +265,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
@@ -244,9 +279,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

@@ -57,6 +57,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,6 +19,7 @@ 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);
@@ -26,6 +27,7 @@ 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 return 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);

Include/pymath.h

+#ifndef Py_PYMATH_H
+#define Py_PYMATH_H
+
+#include "pyconfig.h" /* include for defines */
+
+#ifdef HAVE_STDINT_H
+#include <stdint.h>
+#endif
+
+/**************************************************************************
+Symbols and macros to supply platform-independent interfaces to mathematical
+functions and constants
+**************************************************************************/
+
+/* Python provides implementations for copysign, acosh, asinh, atanh, 
+ * log1p and hypot in Python/pymath.c just in case your math library doesn't
+ * provide the functions.
+ *
+ *Note: PC/pyconfig.h defines copysign as _copysign
+ */
+#ifndef HAVE_COPYSIGN
+extern double copysign(doube, double);
+#endif
+
+#ifndef HAVE_ACOSH
+extern double acosh(double);
+#endif
+
+#ifndef HAVE_ASINH
+extern double asinh(double);
+#endif
+
+#ifndef HAVE_ATANH
+extern double atanh(double);
+#endif
+
+#ifndef HAVE_LOG1P
+extern double log1p(double);
+#endif
+
+#ifndef HAVE_HYPOT
+extern double hypot(double, double);
+#endif
+
+/* extra declarations */
+#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 *);
+extern double pow(double, double);
+#endif /* __STDC__ */
+#endif /* _MSC_VER */
+
+#ifdef _OSF_SOURCE
+/* OSF1 5.1 doesn't make these available with XOPEN_SOURCE_EXTENDED defined */
+extern int finite(double);
+extern double copysign(double, double);
+#endif
+
+/* High precision defintion of pi and e (Euler)
+ * The values are taken from libc6's math.h.
+ */
+#ifndef Py_MATH_PIl
+#define Py_MATH_PIl 3.1415926535897932384626433832795029L
+#endif
+#ifndef Py_MATH_PI
+#define Py_MATH_PI 3.14159265358979323846
+#endif
+
+#ifndef Py_MATH_El
+#define Py_MATH_El 2.7182818284590452353602874713526625L
+#endif
+
+#ifndef Py_MATH_E
+#define Py_MATH_E 2.7182818284590452354
+#endif
+
+/* Py_IS_NAN(X)
+ * Return 1 if float or double arg is a NaN, else 0.
+ * Caution:
+ *     X is evaluated more than once.
+ *     This may not work on all platforms.  Each platform has *some*
+ *     way to spell this, though -- override in pyconfig.h if you have
+ *     a platform where it doesn't work.
+ * Note: PC/pyconfig.h defines Py_IS_NAN as _isnan
+ */
+#ifndef Py_IS_NAN
+#ifdef HAVE_ISNAN
+#define Py_IS_NAN(X) isnan(X)
+#else
+#define Py_IS_NAN(X) ((X) != (X))
+#endif
+#endif
+
+/* Py_IS_INFINITY(X)
+ * Return 1 if float or double arg is an infinity, else 0.
+ * Caution:
+ *    X is evaluated more than once.
+ *    This implementation may set the underflow flag if |X| is very small;
+ *    it really can't be implemented correctly (& easily) before C99.
+ *    Override in pyconfig.h if you have a better spelling on your platform.
+ * Note: PC/pyconfig.h defines Py_IS_INFINITY as _isinf
+ */
+#ifndef Py_IS_INFINITY
+#ifdef HAVE_ISINF
+#define Py_IS_INFINITY(X) isinf(X)
+#else
+#define Py_IS_INFINITY(X) ((X) && (X)*0.5 == (X))
+#endif
+#endif
+
+/* Py_IS_FINITE(X)
+ * Return 1 if float or double arg is neither infinite nor NAN, else 0.
+ * Some compilers (e.g. VisualStudio) have intrisics for this, so a special
+ * macro for this particular test is useful
+ * Note: PC/pyconfig.h defines Py_IS_FINITE as _finite
+ */
+#ifndef Py_IS_FINITE
+#ifdef HAVE_FINITE
+#define Py_IS_FINITE(X) finite(X)
+#else
+#define Py_IS_FINITE(X) (!Py_IS_INFINITY(X) && !Py_IS_NAN(X))
+#endif
+#endif
+
+/* HUGE_VAL is supposed to expand to a positive double infinity.  Python
+ * uses Py_HUGE_VAL instead because some platforms are broken in this
+ * respect.  We used to embed code in pyport.h to try to worm around that,
+ * but different platforms are broken in conflicting ways.  If you're on
+ * a platform where HUGE_VAL is defined incorrectly, fiddle your Python
+ * config to #define Py_HUGE_VAL to something that works on your platform.
+ */
+#ifndef Py_HUGE_VAL
+#define Py_HUGE_VAL HUGE_VAL
+#endif
+
+/* Py_NAN
+ * A value that evaluates to a NaN. On IEEE 754 platforms INF*0 or
+ * INF/INF works. Define Py_NO_NAN in pyconfig.h if your platform
+ * doesn't support NaNs.
+ */
+#if !defined(Py_NAN) && !defined(Py_NO_NAN)
+#define Py_NAN (Py_HUGE_VAL * 0.)
+#endif
+
+/* Py_OVERFLOWED(X)
+ * Return 1 iff a libm function overflowed.  Set errno to 0 before calling
+ * a libm function, and invoke this macro after, passing the function
+ * result.
+ * Caution:
+ *    This isn't reliable.  C99 no longer requires libm to set errno under
+ *	  any exceptional condition, but does require +- HUGE_VAL return
+ *	  values on overflow.  A 754 box *probably* maps HUGE_VAL to a
+ *	  double infinity, and we're cool if that's so, unless the input
+ *	  was an infinity and an infinity is the expected result.  A C89
+ *	  system sets errno to ERANGE, so we check for that too.  We're
+ *	  out of luck if a C99 754 box doesn't map HUGE_VAL to +Inf, or
+ *	  if the returned result is a NaN, or if a C89 box returns HUGE_VAL
+ *	  in non-overflow cases.
+ *    X is evaluated more than once.
+ * Some platforms have better way to spell this, so expect some #ifdef'ery.
+ *
+ * OpenBSD uses 'isinf()' because a compiler bug on that platform causes
+ * the longer macro version to be mis-compiled. This isn't optimal, and
+ * should be removed once a newer compiler is available on that platform.
+ * The system that had the failure was running OpenBSD 3.2 on Intel, with
+ * gcc 2.95.3.
+ *
+ * According to Tim's checkin, the FreeBSD systems use isinf() to work
+ * around a FPE bug on that platform.
+ */
+#if defined(__FreeBSD__) || defined(__OpenBSD__)
+#define Py_OVERFLOWED(X) isinf(X)
+#else
+#define Py_OVERFLOWED(X) ((X) != 0.0 && (errno == ERANGE ||    \
+					 (X) == Py_HUGE_VAL || \
+					 (X) == -Py_HUGE_VAL))
+#endif
+
+#endif /* Py_PYMATH_H */

Include/pyport.h

@@ -336,123 +336,6 @@ extern "C" {
 #define Py_SAFE_DOWNCAST(VALUE, WIDE, NARROW) (NARROW)(VALUE)
 #endif
 
-/* High precision defintion of pi and e (Euler)
- * The values are taken from libc6's math.h.
- */
-#ifndef Py_MATH_PIl
-#define Py_MATH_PIl 3.1415926535897932384626433832795029L
-#endif
-#ifndef Py_MATH_PI
-#define Py_MATH_PI 3.14159265358979323846
-#endif
-
-#ifndef Py_MATH_El
-#define Py_MATH_El 2.7182818284590452353602874713526625L
-#endif
-
-#ifndef Py_MATH_E
-#define Py_MATH_E 2.7182818284590452354
-#endif
-
-/* Py_IS_NAN(X)
- * Return 1 if float or double arg is a NaN, else 0.
- * Caution:
- *     X is evaluated more than once.
- *     This may not work on all platforms.  Each platform has *some*
- *     way to spell this, though -- override in pyconfig.h if you have
- *     a platform where it doesn't work.
- */
-#ifndef Py_IS_NAN
-#ifdef HAVE_ISNAN
-#define Py_IS_NAN(X) isnan(X)
-#else
-#define Py_IS_NAN(X) ((X) != (X))
-#endif
-#endif
-
-/* Py_IS_INFINITY(X)
- * Return 1 if float or double arg is an infinity, else 0.
- * Caution:
- *    X is evaluated more than once.
- *    This implementation may set the underflow flag if |X| is very small;
- *    it really can't be implemented correctly (& easily) before C99.
- *    Override in pyconfig.h if you have a better spelling on your platform.
- */
-#ifndef Py_IS_INFINITY
-#ifdef HAVE_ISINF
-#define Py_IS_INFINITY(X) isinf(X)
-#else
-#define Py_IS_INFINITY(X) ((X) && (X)*0.5 == (X))
-#endif
-#endif
-
-/* Py_IS_FINITE(X)
- * Return 1 if float or double arg is neither infinite nor NAN, else 0.
- * Some compilers (e.g. VisualStudio) have intrisics for this, so a special
- * macro for this particular test is useful
- */
-#ifndef Py_IS_FINITE
-#ifdef HAVE_FINITE
-#define Py_IS_FINITE(X) finite(X)
-#else
-#define Py_IS_FINITE(X) (!Py_IS_INFINITY(X) && !Py_IS_NAN(X))
-#endif
-#endif
-
-/* HUGE_VAL is supposed to expand to a positive double infinity.  Python
- * uses Py_HUGE_VAL instead because some platforms are broken in this
- * respect.  We used to embed code in pyport.h to try to worm around that,
- * but different platforms are broken in conflicting ways.  If you're on
- * a platform where HUGE_VAL is defined incorrectly, fiddle your Python
- * config to #define Py_HUGE_VAL to something that works on your platform.
- */
-#ifndef Py_HUGE_VAL
-#define Py_HUGE_VAL HUGE_VAL
-#endif
-
-/* Py_NAN
- * A value that evaluates to a NaN. On IEEE 754 platforms INF*0 or
- * INF/INF works. Define Py_NO_NAN in pyconfig.h if your platform
- * doesn't support NaNs.
- */
-#if !defined(Py_NAN) && !defined(Py_NO_NAN)
-#define Py_NAN (Py_HUGE_VAL * 0.)
-#endif
-
-/* Py_OVERFLOWED(X)
- * Return 1 iff a libm function overflowed.  Set errno to 0 before calling
- * a libm function, and invoke this macro after, passing the function
- * result.
- * Caution:
- *    This isn't reliable.  C99 no longer requires libm to set errno under
- *	  any exceptional condition, but does require +- HUGE_VAL return
- *	  values on overflow.  A 754 box *probably* maps HUGE_VAL to a
- *	  double infinity, and we're cool if that's so, unless the input
- *	  was an infinity and an infinity is the expected result.  A C89
- *	  system sets errno to ERANGE, so we check for that too.  We're
- *	  out of luck if a C99 754 box doesn't map HUGE_VAL to +Inf, or
- *	  if the returned result is a NaN, or if a C89 box returns HUGE_VAL
- *	  in non-overflow cases.
- *    X is evaluated more than once.
- * Some platforms have better way to spell this, so expect some #ifdef'ery.
- *
- * OpenBSD uses 'isinf()' because a compiler bug on that platform causes
- * the longer macro version to be mis-compiled. This isn't optimal, and
- * should be removed once a newer compiler is available on that platform.
- * The system that had the failure was running OpenBSD 3.2 on Intel, with
- * gcc 2.95.3.
- *
- * According to Tim's checkin, the FreeBSD systems use isinf() to work
- * around a FPE bug on that platform.
- */
-#if defined(__FreeBSD__) || defined(__OpenBSD__)
-#define Py_OVERFLOWED(X) isinf(X)
-#else
-#define Py_OVERFLOWED(X) ((X) != 0.0 && (errno == ERANGE ||    \
-					 (X) == Py_HUGE_VAL || \
-					 (X) == -Py_HUGE_VAL))
-#endif
-
 /* Py_SET_ERRNO_ON_MATH_ERROR(x)
  * If a libm function did not set errno, but it looks like the result
  * overflowed or not-a-number, set errno to ERANGE or EDOM.  Set errno
@@ -559,15 +442,6 @@ extern pid_t forkpty(int *, char *, struct termios *, struct winsize *);
 #endif /* defined(HAVE_OPENPTY) || defined(HAVE_FORKPTY) */
 
 
-/************************
- * WRAPPER FOR <math.h> *
- ************************/
-
-#ifndef HAVE_HYPOT
-extern double hypot(double, double);
-#endif
-
-
 /* On 4.4BSD-descendants, ctype functions serves the whole range of
  * wchar_t character set rather than single byte code points only.
  * This characteristic can break some operations of string object

Lib/test/ieee754.txt

+======================================
+Python IEEE 754 floating point support
+======================================
+
+>>> from sys import float_info as FI
+>>> from math import *
+>>> PI = pi
+>>> E = e
+
+You must never compare two floats with == because you are not going to get
+what you expect. We treat two floats as equal if the difference between them
+is small than epsilon.
+>>> EPS = 1E-15
+>>> def equal(x, y):
+...     """Almost equal helper for floats"""
+...     return abs(x - y) < EPS
+
+
+NaNs and INFs
+=============
+
+In Python 2.6 and newer NaNs (not a number) and infinity can be constructed
+from the strings 'inf' and 'nan'.
+
+>>> INF = float('inf')
+>>> NINF = float('-inf')
+>>> NAN = float('nan')
+
+>>> INF
+inf
+>>> NINF
+-inf
+>>> NAN
+nan
+
+The math module's ``isnan`` and ``isinf`` functions can be used to detect INF
+and NAN:
+>>> isinf(INF), isinf(NINF), isnan(NAN)
+(True, True, True)
+>>> INF == -NINF
+True
+
+Infinity
+--------
+
+Ambiguous operations like ``0 * inf`` or ``inf - inf`` result in NaN.
+>>> INF * 0
+nan
+>>> INF - INF
+nan
+>>> INF / INF
+nan
+
+However unambigous operations with inf return inf:
+>>> INF * INF
+inf
+>>> 1.5 * INF
+inf
+>>> 0.5 * INF
+inf
+>>> INF / 1000
+inf
+
+Not a Number
+------------
+
+NaNs are never equal to another number, even itself
+>>> NAN == NAN
+False
+>>> NAN < 0
+False
+>>> NAN >= 0
+False
+
+All operations involving a NaN return a NaN except for the power of *0* and *1*.
+>>> 1 + NAN
+nan
+>>> 1 * NAN
+nan
+>>> 0 * NAN
+nan
+>>> 1 ** NAN
+1.0
+>>> 0 ** NAN
+0.0
+>>> (1.0 + FI.epsilon) * NAN
+nan
+
+Misc Functions
+==============
+
+The power of 1 raised to x is always 1.0, even for special values like 0,
+infinity and NaN.
+
+>>> pow(1, 0)
+1.0
+>>> pow(1, INF)
+1.0
+>>> pow(1, -INF)
+1.0
+>>> pow(1, NAN)
+1.0
+
+The power of 0 raised to x is defined as 0, if x is positive. Negative
+values are a domain error or zero division error and NaN result in a
+silent NaN.
+
+>>> pow(0, 0)
+1.0
+>>> pow(0, INF)
+0.0
+>>> pow(0, -INF)
+Traceback (most recent call last):
+...
+ValueError: math domain error
+>>> 0 ** -1
+Traceback (most recent call last):
+...
+ZeroDivisionError: 0.0 cannot be raised to a negative power
+>>> pow(0, NAN)
+nan
+
+
+Trigonometric Functions
+=======================
+
+>>> sin(INF)
+Traceback (most recent call last):
+...
+ValueError: math domain error
+>>> sin(NINF)
+Traceback (most recent call last):
+...
+ValueError: math domain error
+>>> sin(NAN)
+nan
+>>> cos(INF)
+Traceback (most recent call last):
+...
+ValueError: math domain error
+>>> cos(NINF)
+Traceback (most recent call last):
+...
+ValueError: math domain error
+>>> cos(NAN)
+nan
+>>> tan(INF)
+Traceback (most recent call last):
+...
+ValueError: math domain error
+>>> tan(NINF)
+Traceback (most recent call last):
+...
+ValueError: math domain error
+>>> tan(NAN)
+nan
+
+Neither pi nor tan are exact, but you can assume that tan(pi/2) is a large value
+and tan(pi) is a very small value:
+>>> tan(PI/2) > 1E10
+True
+>>> -tan(-PI/2) > 1E10
+True
+>>> tan(PI) < 1E-15
+True
+
+>>> asin(NAN), acos(NAN), atan(NAN)
+(nan, nan, nan)
+>>> asin(INF), asin(NINF)
+Traceback (most recent call last):
+...
+ValueError: math domain error
+>>> acos(INF), acos(NINF)
+Traceback (most recent call last):
+...
+ValueError: math domain error
+>>> equal(atan(INF), PI/2), equal(atan(NINF), -PI/2)
+(True, True)
+
+
+Hyberbolic Functions
+====================
+

Lib/test/test_float.py

@@ -2,12 +2,12 @@
 import unittest, struct
 import os
 from test import test_support
+import math
+from math import isinf, isnan
+import operator
 
-def isinf(x):
-    return x * 0.5 == x
-
-def isnan(x):
-    return x != x
+INF = float("inf")
+NAN = float("nan")
 
 class FormatFunctionsTestCase(unittest.TestCase):
 
@@ -239,6 +239,17 @@ class InfNanTest(unittest.TestCase):
         self.assertEqual(str(1e300 * 1e300 * 0), "nan")
         self.assertEqual(str(-1e300 * 1e300 * 0), "nan")
 
+    def notest_float_nan(self):
+        self.assert_(NAN.is_nan())
+        self.failIf(INF.is_nan())
+        self.failIf((0.).is_nan())
+
+    def notest_float_inf(self):
+        self.assert_(INF.is_inf())
+        self.failIf(NAN.is_inf())
+        self.failIf((0.).is_inf())
+
+
 def test_main():
     test_support.run_unittest(
         FormatFunctionsTestCase,

Makefile.pre.in

@@ -276,6 +276,7 @@ PYTHON_OBJS=	\
 		Python/peephole.o \
 		Python/pyarena.o \
 		Python/pyfpe.o \
+		Python/pymath.o \
 		Python/pystate.o \
 		Python/pythonrun.o \
 		Python/structmember.o \
@@ -622,6 +623,7 @@ PYTHON_HEADERS= \
 		Include/pydebug.h \
 		Include/pyerrors.h \
 		Include/pyfpe.h \
+		Include/pymath.h \
 		Include/pygetopt.h \
 		Include/pymem.h \
 		Include/pyport.h \

Objects/complexobject.c

@@ -187,6 +187,38 @@ c_powi(Py_complex x, long n)
 
 }
 
+double
+c_abs(Py_complex z)
+{
+	/* sets errno = ERANGE on overflow;  otherwise errno = 0 */
+	double result;
+
+	if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
+		/* C99 rules: if either the real or the imaginary part is an
+		   infinity, return infinity, even if the other part is a
+		   NaN. */
+		if (Py_IS_INFINITY(z.real)) {
+			result = fabs(z.real);
+			errno = 0;
+			return result;
+		}
+		if (Py_IS_INFINITY(z.imag)) {
+			result = fabs(z.imag);
+			errno = 0;
+			return result;
+		}
+		/* either the real or imaginary part is a NaN,
+		   and neither is infinite. Result should be NaN. */
+		return Py_NAN;
+	}
+	result = hypot(z.real, z.imag);
+	if (!Py_IS_FINITE(result))
+		errno = ERANGE;
+	else
+		errno = 0;
+	return result;
+}
+
 static PyObject *
 complex_subtype_from_c_complex(PyTypeObject *type, Py_complex cval)
 {
@@ -321,8 +353,7 @@ complex_to_buf(char *buf, int bufsz, PyComplexObject *v, int precision)
 		if (!Py_IS_FINITE(v->cval.imag)) {
 			if (Py_IS_NAN(v->cval.imag))
 				strncpy(buf, "nan*j", 6);
-			/* else if (copysign(1, v->cval.imag) == 1) */
-			else if (v->cval.imag > 0)
+			else if (copysign(1, v->cval.imag) == 1)
 				strncpy(buf, "inf*j", 6);
 			else
 				strncpy(buf, "-inf*j", 7);
@@ -578,9 +609,16 @@ static PyObject *
 complex_abs(PyComplexObject *v)
 {
 	double result;
+
 	PyFPE_START_PROTECT("complex_abs", return 0)
-	result = hypot(v->cval.real,v->cval.imag);
+	result = c_abs(v->cval);
 	PyFPE_END_PROTECT(result)
+
+	if (errno == ERANGE) {
+		PyErr_SetString(PyExc_OverflowError,
+				"absolute value too large");
+		return NULL;
+	}
 	return PyFloat_FromDouble(result);
 }
 
@@ -658,9 +696,29 @@ complex_getnewargs(PyComplexObject *v)
 	return Py_BuildValue("(D)", &v->cval);
 }
 
+#if 0
+static PyObject *
+complex_is_finite(PyObject *self)
+{
+	Py_complex c;
+	c = ((PyComplexObject *)self)->cval;
+	return PyBool_FromLong((long)(Py_IS_FINITE(c.real) &&
+				      Py_IS_FINITE(c.imag)));
+}
+
+PyDoc_STRVAR(complex_is_finite_doc,
+"complex.is_finite() -> bool\n"
+"\n"
+"Returns True if the real and the imaginary part is finite.");
+#endif
+
 static PyMethodDef complex_methods[] = {
 	{"conjugate",	(PyCFunction)complex_conjugate,	METH_NOARGS,
 	 complex_conjugate_doc},
+#if 0
+	{"is_finite",	(PyCFunction)complex_is_finite,	METH_NOARGS,
+	 complex_is_finite_doc},
+#endif
 	{"__getnewargs__",	(PyCFunction)complex_getnewargs,	METH_NOARGS},
 	{NULL,		NULL}		/* sentinel */
 };

Objects/floatobject.c

@@ -16,10 +16,6 @@
 #include <ieeefp.h>
 #endif
 
-#if !defined(__STDC__)
-extern double fmod(double, double);
-extern double pow(double, double);
-#endif
 
 #ifdef _OSF_SOURCE
 /* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */
@@ -224,11 +220,11 @@ PyFloat_FromString(PyObject *v)
 			p++;
 		}
 		if (PyOS_strnicmp(p, "inf", 4) == 0) {
-			return PyFloat_FromDouble(sign * Py_HUGE_VAL);
+			Py_RETURN_INF(sign);
 		}
 #ifdef Py_NAN
 		if(PyOS_strnicmp(p, "nan", 4) == 0) {
-			return PyFloat_FromDouble(Py_NAN);
+			Py_RETURN_NAN;
 		}
 #endif
 		PyOS_snprintf(buffer, sizeof(buffer),
@@ -378,110 +374,6 @@ format_float(char *buf, size_t buflen, PyFloatObject *v, int precision)
 	format_double(buf, buflen, PyFloat_AS_DOUBLE(v), precision);
 }
 
-#ifdef Py_BROKEN_REPR
-/* The following function is based on Tcl_PrintDouble,
- * from tclUtil.c.
- */
-
-#define is_infinite(d)	( (d) > DBL_MAX || (d) < -DBL_MAX )
-#define is_nan(d)		((d) != (d))
-
-static void
-format_double_repr(char *dst, double value)
-{
-    char *p, c;
-    int exp;
-    int signum;
-    char buffer[30];
-
-	/*
-	 * Handle NaN.
-	 */
-
-	if (is_nan(value)) {
-	    strcpy(dst, "nan");
-	    return;
-	}
-
-	/*
-	 * Handle infinities.
-	 */
-
-	if (is_infinite(value)) {
-	    if (value < 0) {
-		strcpy(dst, "-inf");
-	    } else {
-		strcpy(dst, "inf");
-	    }
-	    return;
-	}
-
-	/*
-	 * Ordinary (normal and denormal) values.
-	 */
-
-	exp = _PyFloat_Digits(buffer, value, &signum)+1;
-	if (signum) {
-	    *dst++ = '-';
-	}
-	p = buffer;
-	if (exp < -3 || exp > 17) {
-	    /*
-	     * E format for numbers < 1e-3 or >= 1e17.
-	     */
-
-	    *dst++ = *p++;
-	    c = *p;
-	    if (c != '\0') {
-		*dst++ = '.';
-		while (c != '\0') {
-		    *dst++ = c;
-		    c = *++p;
-		}
-	    }
-	    sprintf(dst, "e%+d", exp-1);
-	} else {
-	    /*
-	     * F format for others.
-	     */
-
-	    if (exp <= 0) {
-		*dst++ = '0';
-	    }
-	    c = *p;
-	    while (exp-- > 0) {
-		if (c != '\0') {
-		    *dst++ = c;
-		    c = *++p;
-		} else {
-		    *dst++ = '0';
-		}
-	    }
-	    *dst++ = '.';
-	    if (c == '\0') {
-		*dst++ = '0';
-	    } else {
-		while (++exp < 0) {
-		    *dst++ = '0';
-		}
-		while (c != '\0') {
-		    *dst++ = c;
-		    c = *++p;
-		}
-	    }
-	    *dst++ = '\0';
-	}
-}
-
-static void
-format_float_repr(char *buf, PyFloatObject *v)
-{
-	assert(PyFloat_Check(v));
-	format_double_repr(buf, PyFloat_AS_DOUBLE(v));
-}
-
-#endif /* Py_BROKEN_REPR */
-
 /* Macro and helper that convert PyObject obj to a C double and store
    the value in dbl.  If conversion to double raises an exception, obj is
    set to NULL, and the function invoking this macro returns NULL.  If
@@ -534,13 +426,8 @@ convert_to_double(PyObject **v, double *dbl)
 static PyObject *
 float_repr(PyFloatObject *v)
 {
-#ifdef Py_BROKEN_REPR
-	char buf[30];
-	format_float_repr(buf, v);
-#else
 	char buf[100];
 	format_float(buf, sizeof(buf), v, PREC_REPR);
-#endif
 
 	return PyUnicode_FromString(buf);
 }
@@ -804,10 +691,13 @@ float_div(PyObject *v, PyObject *w)
 	double a,b;
 	CONVERT_TO_DOUBLE(v, a);
 	CONVERT_TO_DOUBLE(w, b);
+#ifdef Py_NAN
 	if (b == 0.0) {
-		PyErr_SetString(PyExc_ZeroDivisionError, "float division");
+		PyErr_SetString(PyExc_ZeroDivisionError,
+				"float division");
 		return NULL;
 	}
+#endif
 	PyFPE_START_PROTECT("divide", return 0)
 	a = a / b;
 	PyFPE_END_PROTECT(a)
@@ -819,12 +709,15 @@ float_rem(PyObject *v, PyObject *w)
 {
 	double vx, wx;
 	double mod;
- 	CONVERT_TO_DOUBLE(v, vx);
- 	CONVERT_TO_DOUBLE(w, wx);
+	CONVERT_TO_DOUBLE(v, vx);
+	CONVERT_TO_DOUBLE(w, wx);
+#ifdef Py_NAN
 	if (wx == 0.0) {
-		PyErr_SetString(PyExc_ZeroDivisionError, "float modulo");
+		PyErr_SetString(PyExc_ZeroDivisionError,
+				"float modulo");
 		return NULL;
 	}
+#endif
 	PyFPE_START_PROTECT("modulo", return 0)
 	mod = fmod(vx, wx);
 	/* note: checking mod*wx < 0 is incorrect -- underflows to
@@ -928,6 +821,9 @@ float_pow(PyObject *v, PyObject *w, PyObject *z)
 		}
 		return PyFloat_FromDouble(0.0);
 	}
+	if (iv == 1.0) { /* 1**w is 1, even 1**inf and 1**nan */
+		return PyFloat_FromDouble(1.0);
+	}
 	if (iv < 0.0) {
 		/* Whether this is an error is a mess, and bumps into libm
 		 * bugs so we have to figure it out ourselves.
@@ -994,6 +890,57 @@ float_bool(PyFloatObject *v)
 	return v->ob_fval != 0.0;
 }
 
+static PyObject *
+float_is_integer(PyObject *v)
+{
+	double x = PyFloat_AsDouble(v);
+	PyObject *o;
+	
+	if (x == -1.0 && PyErr_Occurred())
+		return NULL;
+	if (!Py_IS_FINITE(x))
+		Py_RETURN_FALSE;
+	PyFPE_START_PROTECT("is_integer", return NULL)
+	o = (floor(x) == x) ? Py_True : Py_False;
+	PyFPE_END_PROTECT(x)
+	if (errno != 0) {
+		PyErr_SetFromErrno(errno == ERANGE ? PyExc_OverflowError :
+						     PyExc_ValueError);
+		return NULL;
+	}
+	Py_INCREF(o);
+	return o;
+}
+
+#if 0
+static PyObject *
+float_is_inf(PyObject *v)
+{
+	double x = PyFloat_AsDouble(v);
+	if (x == -1.0 && PyErr_Occurred())
+		return NULL;
+	return PyBool_FromLong((long)Py_IS_INFINITY(x));
+}
+
+static PyObject *
+float_is_nan(PyObject *v)
+{
+	double x = PyFloat_AsDouble(v);
+	if (x == -1.0 && PyErr_Occurred())
+		return NULL;
+	return PyBool_FromLong((long)Py_IS_NAN(x));
+}
+
+static PyObject *
+float_is_finite(PyObject *v)
+{
+	double x = PyFloat_AsDouble(v);
+	if (x == -1.0 && PyErr_Occurred())
+		return NULL;
+	return PyBool_FromLong((long)Py_IS_FINITE(x));
+}
+#endif
+
 static PyObject *
 float_trunc(PyObject *v)
 {
@@ -1368,7 +1315,7 @@ PyDoc_STRVAR(float__format__doc,
 
 
 static PyMethodDef float_methods[] = {
-  	{"conjugate",	(PyCFunction)float_float,	METH_NOARGS,
+	{"conjugate",	(PyCFunction)float_float,	METH_NOARGS,
 	 "Returns self, the complex conjugate of any float."},
 	{"__trunc__",	(PyCFunction)float_trunc, METH_NOARGS,
          "Returns the Integral closest to x between 0 and x."},
@@ -1377,6 +1324,16 @@ static PyMethodDef float_methods[] = {
          "When an argument is passed, works like built-in round(x, ndigits)."},
 	{"as_integer_ratio", (PyCFunction)float_as_integer_ratio, METH_NOARGS,
 	 float_as_integer_ratio_doc},
+	{"is_integer",	(PyCFunction)float_is_integer,	METH_NOARGS,
+	 "Returns True if the float is an integer."},
+#if 0
+	{"is_inf",	(PyCFunction)float_is_inf,	METH_NOARGS,
+	 "Returns True if the float is positive or negative infinite."},
+	{"is_finite",	(PyCFunction)float_is_finite,	METH_NOARGS,
+	 "Returns True if the float is finite, neither infinite nor NaN."},
+	{"is_nan",	(PyCFunction)float_is_nan,	METH_NOARGS,
+	 "Returns True if the float is not a number (NaN)."},
+#endif
 	{"__getnewargs__",	(PyCFunction)float_getnewargs,	METH_NOARGS},
 	{"__getformat__",	(PyCFunction)float_getformat,	
 	 METH_O|METH_CLASS,		float_getformat_doc},
@@ -1534,10 +1491,6 @@ _PyFloat_Init(void)
 	double_format = detected_double_format;
 	float_format = detected_float_format;
 
-#ifdef Py_BROKEN_REPR	
-	/* Initialize floating point repr */
-	_PyFloat_DigitsInit();
-#endif
 	/* Init float info */
 	if (FloatInfoType.tp_name == 0)
 		PyStructSequence_InitType(&FloatInfoType, &floatinfo_desc);

Objects/longobject.c

@@ -3611,9 +3611,21 @@ long_round(PyObject *self, PyObject *args)
 #undef UNDEF_NDIGITS
 }
 
+#if 0
+static PyObject *
+long_is_finite(PyObject *v)
+{
+	Py_RETURN_TRUE;
+}
+#endif
+
 static PyMethodDef long_methods[] = {
 	{"conjugate",	(PyCFunction)long_long,	METH_NOARGS,
 	 "Returns self, the complex conjugate of any int."},
+#if 0
+	{"is_finite",	(PyCFunction)long_is_finite,	METH_NOARGS,
+	 "Returns always True."},
+#endif
 	{"__trunc__",	(PyCFunction)long_long,	METH_NOARGS,
          "Truncating an Integral returns itself."},
 	{"__floor__",	(PyCFunction)long_long,	METH_NOARGS,

PC/VC6/pythoncore.dsp

@@ -587,6 +587,10 @@ SOURCE=..\..\Python\pyfpe.c
 # End Source File
 # Begin Source File
 
+SOURCE=..\..\Python\pymath.c
+# End Source File
+# Begin Source File
+
 SOURCE=..\..\Python\pystate.c
 # End Source File
 # Begin Source File

PC/VS7.1/pythoncore.vcproj

@@ -697,6 +697,9 @@
 		<File
 			RelativePath="..\..\Python\pyfpe.c">
 		</File>
+		<File
+			RelativePath="..\..\Python\pymath.c">
+		</File>
 		<File
 			RelativePath="..\..\Python\pystate.c">
 		</File>

PC/VS8.0/pythoncore.vcproj

@@ -1706,6 +1706,10 @@
 				RelativePath="..\..\Python\pyfpe.c"
 				>
 			</File>
+			<File
+				RelativePath="..\..\Python\pymath.c"
+				>
+			</File>
 			<File
 				RelativePath="..\..\Python\pystate.c"
 				>

PC/pyconfig.h

@@ -207,12 +207,13 @@ typedef _W64 int ssize_t;
 #endif /* MS_WIN32 && !MS_WIN64 */
 
 typedef int pid_t;
-#define hypot _hypot
 
 #include <float.h>
 #define Py_IS_NAN _isnan
 #define Py_IS_INFINITY(X) (!_finite(X) && !_isnan(X))
 #define Py_IS_FINITE(X) _finite(X)
+#define copysign _copysign
+#define hypot _hypot
 
 #endif /* _MSC_VER */
 
@@ -392,7 +393,7 @@ Py_NO_ENABLE_SHARED to find out.  Also support MS_NO_COREDLL for b/w compat */
 /* Fairly standard from here! */
 
 /* Define to 1 if you have the `copysign' function. */
-/* #define HAVE_COPYSIGN 1*/
+#define HAVE_COPYSIGN 1
 
 /* Define to 1 if you have the `isinf' function. */
 #define HAVE_ISINF 1

PCbuild/pythoncore.vcproj

@@ -870,6 +870,10 @@
 				RelativePath="..\Include\pymactoolbox.h"
 				>
 			</File>
+			<File
+				RelativePath="..\Include\pymath.h"
+				>
+			</File>
 			<File
 				RelativePath="..\Include\pymem.h"
 				>
@@ -1706,6 +1710,10 @@
 				RelativePath="..\Python\pyfpe.c"
 				>
 			</File>
+			<File
+				RelativePath="..\Python\pymath.c"
+				>
+			</File>
 			<File
 				RelativePath="..\Python\pystate.c"
 				>

Python/hypot.c

-/* hypot() replacement */
-
-#include "Python.h"
-
-#ifndef HAVE_HYPOT
-double hypot(double x, double y)
-{
-	double yx;
-
-	x = fabs(x);
-	y = fabs(y);
-	if (x < y) {
-		double temp = x;
-		x = y;
-		y = temp;
-	}
-	if (x == 0.)
-		return 0.;
-	else {
-		yx = y/x;
-		return x*sqrt(1.+yx*yx);
-	}
-}
-#endif /* HAVE_HYPOT */
-

Python/pymath.c

+#include "Python.h"
+
+#ifndef HAVE_HYPOT
+double hypot(double x, double y)
+{
+	double yx;
+
+	x = fabs(x);
+	y = fabs(y);
+	if (x < y) {
+		double temp = x;
+		x = y;
+		y = temp;
+	}
+	if (x == 0.)
+		return 0.;
+	else {
+		yx = y/x;
+		return x*sqrt(1.+yx*yx);
+	}
+}
+#endif /* HAVE_HYPOT */
+
+#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 /* HAVE_COPYSIGN */
+
+#ifndef HAVE_LOG1P
+double
+log1p(double x)
+{
+	/* For x small, we use the following approach.  Let y be the nearest
+	   float to 1+x, then
+
+	     1+x = y * (1 - (y-1-x)/y)
+
+	   so log(1+x) = log(y) + log(1-(y-1-x)/y).  Since (y-1-x)/y is tiny,
+	   the second term is well approximated by (y-1-x)/y.  If abs(x) >=
+	   DBL_EPSILON/2 or the rounding-mode is some form of round-to-nearest
+	   then y-1-x will be exactly representable, and is computed exactly
+	   by (y-1)-x.
+
+	   If abs(x) < DBL_EPSILON/2 and the rounding mode is not known to be
+	   round-to-nearest then this method is slightly dangerous: 1+x could
+	   be rounded up to 1+DBL_EPSILON instead of down to 1, and in that
+	   case y-1-x will not be exactly representable any more and the
+	   result can be off by many ulps.  But this is easily fixed: for a
+	   floating-point number |x| < DBL_EPSILON/2., the closest
+	   floating-point number to log(1+x) is exactly x.
+	*/
+
+	double y;
+	if (fabs(x) < DBL_EPSILON/2.) {
+		return x;
+	} else if (-0.5 <= x && x <= 1.) {
+		/* WARNING: it's possible than an overeager 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 /* HAVE_LOG1P */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+static const double ln2 = 6.93147180559945286227E-01;
+static const double two_pow_m28 = 3.7252902984619141E-09; /* 2**-28 */
+static const double two_pow_p28 = 268435456.0; /* 2**28 */
+static const double zero = 0.0;
+
+/* asinh(x)
+ * Method :
+ *	Based on 
+ *		asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
+ *	we have
+ *	asinh(x) := x  if  1+x*x=1,
+ *		 := sign(x)*(log(x)+ln2)) for large |x|, else
+ *		 := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
+ *		 := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))  
+ */
+
+#ifndef HAVE_ASINH
+double
+asinh(double x)
+{	
+	double w;
+	double absx = fabs(x);
+
+	if (Py_IS_NAN(x) || Py_IS_INFINITY(x)) {
+		return x+x;
+	}
+	if (absx < two_pow_m28) {	/* |x| < 2**-28 */
+		return x;	/* return x inexact except 0 */
+	} 
+	if (absx > two_pow_p28) {	/* |x| > 2**28 */
+		w = log(absx)+ln2;
+	}
+	else if (absx > 2.0) {		/* 2 < |x| < 2**28 */
+		w = log(2.0*absx + 1.0 / (sqrt(x*x + 1.0) + absx));
+	}
+	else {				/* 2**-28 <= |x| < 2= */
+		double t = x*x;
+		w = log1p(absx + t / (1.0 + sqrt(1.0 + t)));
+	}
+	return copysign(w, x);
+	
+}
+#endif /* HAVE_ASINH */
+
+/* acosh(x)
+ * Method :
+ *      Based on
+ *	      acosh(x) = log [ x + sqrt(x*x-1) ]
+ *      we have
+ *	      acosh(x) := log(x)+ln2, if x is large; else
+ *	      acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
+ *	      acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
+ *
+ * Special cases:
+ *      acosh(x) is NaN with signal if x<1.
+ *      acosh(NaN) is NaN without signal.
+ */
+
+#ifndef HAVE_ACOSH
+double
+acosh(double x)
+{
+	if (Py_IS_NAN(x)) {
+		return x+x;
+	}
+	if (x < 1.) {			/* x < 1;  return a signaling NaN */
+		errno = EDOM;
+#ifdef Py_NAN
+		return Py_NAN;
+#else
+		return (x-x)/(x-x);
+#endif
+	}
+	else if (x >= two_pow_p28) {	/* x > 2**28 */
+		if (Py_IS_INFINITY(x)) {
+			return x+x;
+		} else {
+			return log(x)+ln2;	/* acosh(huge)=log(2x) */
+		}
+	}
+	else if (x == 1.) {
+		return 0.0;			/* acosh(1) = 0 */
+	}
+	else if (x > 2.) {			/* 2 < x < 2**28 */
+		double t = x*x;
+		return log(2.0*x - 1.0 / (x + sqrt(t - 1.0)));
+	}
+	else {				/* 1 < x <= 2 */
+		double t = x - 1.0;
+		return log1p(t + sqrt(2.0*t + t*t));
+	}
+}
+#endif /* HAVE_ACOSH */
+
+/* atanh(x)
+ * Method :
+ *    1.Reduced x to positive by atanh(-x) = -atanh(x)
+ *    2.For x>=0.5
+ *		  1	      2x			  x
+ *      atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
+ *		  2	     1 - x		      1 - x
+ *
+ *      For x<0.5
+ *      atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
+ *
+ * Special cases:
+ *      atanh(x) is NaN if |x| >= 1 with signal;
+ *      atanh(NaN) is that NaN with no signal;
+ *
+ */
+
+#ifndef HAVE_ATANH
+double
+atanh(double x)
+{
+	double absx;
+	double t;
+
+	if (Py_IS_NAN(x)) {
+		return x+x;
+	}
+	absx = fabs(x);
+	if (absx >= 1.) {		/* |x| >= 1 */
+		errno = EDOM;
+#ifdef Py_NAN
+		return Py_NAN;
+#else
+		return x/zero;
+#endif
+	}
+	if (absx < two_pow_m28) {	/* |x| < 2**-28 */
+		return x;
+	}
+	if (absx < 0.5) {		/* |x| < 0.5 */
+		t = absx+absx;
+		t = 0.5 * log1p(t + t*absx / (1.0 - absx));
+	} 
+	else {				/* 0.5 <= |x| <= 1.0 */
+		t = 0.5 * log1p((absx + absx) / (1.0 - absx));
+	}
+	return copysign(t, x);
+}
+#endif /* HAVE_ATANH */

configure

 #! /bin/sh
-# From configure.in Revision: 62003 .
+# From configure.in Revision: 62146 .
 # Guess values for system-dependent variables and create Makefiles.
 # Generated by GNU Autoconf 2.61 for python 3.0.
 #