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 :)

5 files changed, 196 insertions, 126 deletions

diff --git a/Include/Python.h b/Include/Python.h
index 2ff4c31..60258ef 100644
--- a/Include/Python.h
+++ b/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"
diff --git a/Include/complexobject.h b/Include/complexobject.h
index b036444..84b6d8b 100644
--- a/Include/complexobject.h
+++ b/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 */
diff --git a/Include/floatobject.h b/Include/floatobject.h
index 0e5f6a9..91814f2 100644
--- a/Include/floatobject.h
+++ b/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);
diff --git a/Include/pymath.h b/Include/pymath.h
new file mode 100644
index 0000000..a3735c2
--- /dev/null
+++ b/Include/pymath.h
@@ -0,0 +1,182 @@
+#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 */
diff --git a/Include/pyport.h b/Include/pyport.h
index ce760d2..4656f8c 100644
--- a/Include/pyport.h
+++ b/Include/pyport.h
@@ -353,123 +353,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
@@ -601,15 +484,6 @@ extern int fdatasync(int);
 #endif /* 0 */
 
 
-/************************
- * 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