[Python-checkins] r62384 - in python/branches/py3k: Doc/library/cmath.rst Doc/library/math.rst Include/Python.h Include/complexobject.h Include/floatobject.h Include/pymath.h Include/pyport.h Lib/test/cmath_testcases.txt Lib/test/ieee754.txt Lib/test/test_cmath.py Lib/test/test_float.py Lib/test/test_math.py Makefile.pre.in Modules/cmathmodule.c Modules/mathmodule.c Objects/complexobject.c Objects/doubledigits.c Objects/floatobject.c Objects/longobject.c PC/VC6/pythoncore.dsp PC/VS7.1/pythoncore.vcproj PC/VS8.0/pythoncore.vcproj PC/pyconfig.h PCbuild/pythoncore.vcproj Python/hypot.c Python/pymath.c configure
christian.heimes
python-checkins at python.org
Sat Apr 19 02:31:42 CEST 2008
Author: christian.heimes
Date: Sat Apr 19 02:31:39 2008
New Revision: 62384
Log:
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 ...
........
Added:
python/branches/py3k/Include/pymath.h
- copied unchanged from r62380, /python/trunk/Include/pymath.h
python/branches/py3k/Lib/test/cmath_testcases.txt
- copied unchanged from r62380, /python/trunk/Lib/test/cmath_testcases.txt
python/branches/py3k/Lib/test/ieee754.txt
- copied unchanged from r62380, /python/trunk/Lib/test/ieee754.txt
python/branches/py3k/Python/pymath.c
- copied unchanged from r62380, /python/trunk/Python/pymath.c
Removed:
python/branches/py3k/Objects/doubledigits.c
python/branches/py3k/Python/hypot.c
Modified:
python/branches/py3k/ (props changed)
python/branches/py3k/Doc/library/cmath.rst
python/branches/py3k/Doc/library/math.rst
python/branches/py3k/Include/Python.h
python/branches/py3k/Include/complexobject.h
python/branches/py3k/Include/floatobject.h
python/branches/py3k/Include/pyport.h
python/branches/py3k/Lib/test/test_cmath.py
python/branches/py3k/Lib/test/test_float.py
python/branches/py3k/Lib/test/test_math.py
python/branches/py3k/Makefile.pre.in
python/branches/py3k/Modules/cmathmodule.c
python/branches/py3k/Modules/mathmodule.c
python/branches/py3k/Objects/complexobject.c
python/branches/py3k/Objects/floatobject.c
python/branches/py3k/Objects/longobject.c
python/branches/py3k/PC/VC6/pythoncore.dsp
python/branches/py3k/PC/VS7.1/pythoncore.vcproj
python/branches/py3k/PC/VS8.0/pythoncore.vcproj
python/branches/py3k/PC/pyconfig.h
python/branches/py3k/PCbuild/pythoncore.vcproj
python/branches/py3k/configure
Modified: python/branches/py3k/Doc/library/cmath.rst
==============================================================================
--- python/branches/py3k/Doc/library/cmath.rst (original)
+++ python/branches/py3k/Doc/library/cmath.rst Sat Apr 19 02:31:39 2008
@@ -14,8 +14,81 @@
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 @@
.. 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 @@
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 @@
nothing's sign bit. In Iserles, A., and Powell, M. (eds.), The state of the art
in numerical analysis. Clarendon Press (1987) pp165-211.
+
Modified: python/branches/py3k/Doc/library/math.rst
==============================================================================
--- python/branches/py3k/Doc/library/math.rst (original)
+++ python/branches/py3k/Doc/library/math.rst Sat Apr 19 02:31:39 2008
@@ -128,6 +128,14 @@
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 @@
.. 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 @@
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 @@
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 @@
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 mathematical constant *e*.
+
.. note::
The :mod:`math` module consists mostly of thin wrappers around the platform C
@@ -244,9 +279,17 @@
: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.
-
Modified: python/branches/py3k/Include/Python.h
==============================================================================
--- python/branches/py3k/Include/Python.h (original)
+++ python/branches/py3k/Include/Python.h Sat Apr 19 02:31:39 2008
@@ -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"
Modified: python/branches/py3k/Include/complexobject.h
==============================================================================
--- python/branches/py3k/Include/complexobject.h (original)
+++ python/branches/py3k/Include/complexobject.h Sat Apr 19 02:31:39 2008
@@ -19,6 +19,7 @@
#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_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 */
Modified: python/branches/py3k/Include/floatobject.h
==============================================================================
--- python/branches/py3k/Include/floatobject.h (original)
+++ python/branches/py3k/Include/floatobject.h Sat Apr 19 02:31:39 2008
@@ -21,6 +21,17 @@
#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);
Modified: python/branches/py3k/Include/pyport.h
==============================================================================
--- python/branches/py3k/Include/pyport.h (original)
+++ python/branches/py3k/Include/pyport.h Sat Apr 19 02:31:39 2008
@@ -336,123 +336,6 @@
#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 @@
#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
Modified: python/branches/py3k/Lib/test/test_cmath.py
==============================================================================
--- python/branches/py3k/Lib/test/test_cmath.py (original)
+++ python/branches/py3k/Lib/test/test_cmath.py Sat Apr 19 02:31:39 2008
@@ -1,6 +1,81 @@
from test.test_support import run_unittest
+from test.test_math import parse_testfile, test_file
import unittest
+import os, sys
import cmath, math
+from cmath import phase, polar, rect, pi
+
+INF = float('inf')
+NAN = float('nan')
+
+complex_zeros = [complex(x, y) for x in [0.0, -0.0] for y in [0.0, -0.0]]
+complex_infinities = [complex(x, y) for x, y in [
+ (INF, 0.0), # 1st quadrant
+ (INF, 2.3),
+ (INF, INF),
+ (2.3, INF),
+ (0.0, INF),
+ (-0.0, INF), # 2nd quadrant
+ (-2.3, INF),
+ (-INF, INF),
+ (-INF, 2.3),
+ (-INF, 0.0),
+ (-INF, -0.0), # 3rd quadrant
+ (-INF, -2.3),
+ (-INF, -INF),
+ (-2.3, -INF),
+ (-0.0, -INF),
+ (0.0, -INF), # 4th quadrant
+ (2.3, -INF),
+ (INF, -INF),
+ (INF, -2.3),
+ (INF, -0.0)
+ ]]
+complex_nans = [complex(x, y) for x, y in [
+ (NAN, -INF),
+ (NAN, -2.3),
+ (NAN, -0.0),
+ (NAN, 0.0),
+ (NAN, 2.3),
+ (NAN, INF),
+ (-INF, NAN),
+ (-2.3, NAN),
+ (-0.0, NAN),
+ (0.0, NAN),
+ (2.3, NAN),
+ (INF, NAN)
+ ]]
+
+def almostEqualF(a, b, rel_err=2e-15, abs_err = 5e-323):
+ """Determine whether floating-point values a and b are equal to within
+ a (small) rounding error. The default values for rel_err and
+ abs_err are chosen to be suitable for platforms where a float is
+ represented by an IEEE 754 double. They allow an error of between
+ 9 and 19 ulps."""
+
+ # special values testing
+ if math.isnan(a):
+ return math.isnan(b)
+ if math.isinf(a):
+ return a == b
+
+ # if both a and b are zero, check whether they have the same sign
+ # (in theory there are examples where it would be legitimate for a
+ # and b to have opposite signs; in practice these hardly ever
+ # occur).
+ if not a and not b:
+ return math.copysign(1., a) == math.copysign(1., b)
+
+ # if a-b overflows, or b is infinite, return False. Again, in
+ # theory there are examples where a is within a few ulps of the
+ # max representable float, and then b could legitimately be
+ # infinite. In practice these examples are rare.
+ try:
+ absolute_error = abs(b-a)
+ except OverflowError:
+ return False
+ else:
+ return absolute_error <= max(abs_err, rel_err * abs(a))
class CMathTests(unittest.TestCase):
# list of all functions in cmath
@@ -12,25 +87,51 @@
test_functions.append(lambda x : cmath.log(x, 1729. + 0j))
test_functions.append(lambda x : cmath.log(14.-27j, x))
- def cAssertAlmostEqual(self, a, b, rel_eps = 1e-10, abs_eps = 1e-100):
- """Check that two complex numbers are almost equal."""
- # the two complex numbers are considered almost equal if
- # either the relative error is <= rel_eps or the absolute error
- # is tiny, <= abs_eps.
- if a == b == 0:
- return
- absolute_error = abs(a-b)
- relative_error = absolute_error/max(abs(a), abs(b))
- if relative_error > rel_eps and absolute_error > abs_eps:
- self.fail("%s and %s are not almost equal" % (a, b))
+ def setUp(self):
+ self.test_values = open(test_file)
+
+ def tearDown(self):
+ self.test_values.close()
+
+ def rAssertAlmostEqual(self, a, b, rel_err = 2e-15, abs_err = 5e-323):
+ """Check that two floating-point numbers are almost equal."""
+
+ # special values testing
+ if math.isnan(a):
+ if math.isnan(b):
+ return
+ self.fail("%s should be nan" % repr(b))
+
+ if math.isinf(a):
+ if a == b:
+ return
+ self.fail("finite result where infinity excpected: "
+ "expected %s, got %s" % (repr(a), repr(b)))
+
+ if not a and not b:
+ if math.atan2(a, -1.) != math.atan2(b, -1.):
+ self.fail("zero has wrong sign: expected %s, got %s" %
+ (repr(a), repr(b)))
+
+ # test passes if either the absolute error or the relative
+ # error is sufficiently small. The defaults amount to an
+ # error of between 9 ulps and 19 ulps on an IEEE-754 compliant
+ # machine.
+
+ try:
+ absolute_error = abs(b-a)
+ except OverflowError:
+ pass
+ else:
+ if absolute_error <= max(abs_err, rel_err * abs(a)):
+ return
+ self.fail("%s and %s are not sufficiently close" % (repr(a), repr(b)))
def test_constants(self):
e_expected = 2.71828182845904523536
pi_expected = 3.14159265358979323846
- self.assertAlmostEqual(cmath.pi, pi_expected, places=9,
- msg="cmath.pi is %s; should be %s" % (cmath.pi, pi_expected))
- self.assertAlmostEqual(cmath.e, e_expected, places=9,
- msg="cmath.e is %s; should be %s" % (cmath.e, e_expected))
+ self.assertAlmostEqual(cmath.pi, pi_expected)
+ self.assertAlmostEqual(cmath.e, e_expected)
def test_user_object(self):
# Test automatic calling of __complex__ and __float__ by cmath
@@ -109,13 +210,13 @@
for f in self.test_functions:
# usual usage
- self.cAssertAlmostEqual(f(MyComplex(cx_arg)), f(cx_arg))
- self.cAssertAlmostEqual(f(MyComplexOS(cx_arg)), f(cx_arg))
+ self.assertEqual(f(MyComplex(cx_arg)), f(cx_arg))
+ self.assertEqual(f(MyComplexOS(cx_arg)), f(cx_arg))
# other combinations of __float__ and __complex__
- self.cAssertAlmostEqual(f(FloatAndComplex()), f(cx_arg))
- self.cAssertAlmostEqual(f(FloatAndComplexOS()), f(cx_arg))
- self.cAssertAlmostEqual(f(JustFloat()), f(flt_arg))
- self.cAssertAlmostEqual(f(JustFloatOS()), f(flt_arg))
+ self.assertEqual(f(FloatAndComplex()), f(cx_arg))
+ self.assertEqual(f(FloatAndComplexOS()), f(cx_arg))
+ self.assertEqual(f(JustFloat()), f(flt_arg))
+ self.assertEqual(f(JustFloatOS()), f(flt_arg))
# TypeError should be raised for classes not providing
# either __complex__ or __float__, even if they provide
# __int__, __long__ or __index__. An old-style class
@@ -138,7 +239,7 @@
# functions, by virtue of providing a __float__ method
for f in self.test_functions:
for arg in [2, 2.]:
- self.cAssertAlmostEqual(f(arg), f(arg.__float__()))
+ self.assertEqual(f(arg), f(arg.__float__()))
# but strings should give a TypeError
for f in self.test_functions:
@@ -182,12 +283,201 @@
float_fn = getattr(math, fn)
complex_fn = getattr(cmath, fn)
for v in values:
- self.cAssertAlmostEqual(float_fn(v), complex_fn(v))
+ z = complex_fn(v)
+ self.rAssertAlmostEqual(float_fn(v), z.real)
+ self.assertEqual(0., z.imag)
# test two-argument version of log with various bases
for base in [0.5, 2., 10.]:
for v in positive:
- self.cAssertAlmostEqual(cmath.log(v, base), math.log(v, base))
+ z = cmath.log(v, base)
+ self.rAssertAlmostEqual(math.log(v, base), z.real)
+ self.assertEqual(0., z.imag)
+
+ def test_specific_values(self):
+ if not float.__getformat__("double").startswith("IEEE"):
+ return
+
+ def rect_complex(z):
+ """Wrapped version of rect that accepts a complex number instead of
+ two float arguments."""
+ return cmath.rect(z.real, z.imag)
+
+ def polar_complex(z):
+ """Wrapped version of polar that returns a complex number instead of
+ two floats."""
+ return complex(*polar(z))
+
+ for id, fn, ar, ai, er, ei, flags in parse_testfile(test_file):
+ arg = complex(ar, ai)
+ expected = complex(er, ei)
+ if fn == 'rect':
+ function = rect_complex
+ elif fn == 'polar':
+ function = polar_complex
+ else:
+ function = getattr(cmath, fn)
+ if 'divide-by-zero' in flags or 'invalid' in flags:
+ try:
+ actual = function(arg)
+ except ValueError:
+ continue
+ else:
+ test_str = "%s: %s(complex(%r, %r))" % (id, fn, ar, ai)
+ self.fail('ValueError not raised in test %s' % test_str)
+
+ if 'overflow' in flags:
+ try:
+ actual = function(arg)
+ except OverflowError:
+ continue
+ else:
+ test_str = "%s: %s(complex(%r, %r))" % (id, fn, ar, ai)
+ self.fail('OverflowError not raised in test %s' % test_str)
+
+ actual = function(arg)
+
+ if 'ignore-real-sign' in flags:
+ actual = complex(abs(actual.real), actual.imag)
+ expected = complex(abs(expected.real), expected.imag)
+ if 'ignore-imag-sign' in flags:
+ actual = complex(actual.real, abs(actual.imag))
+ expected = complex(expected.real, abs(expected.imag))
+
+ # for the real part of the log function, we allow an
+ # absolute error of up to 2e-15.
+ if fn in ('log', 'log10'):
+ real_abs_err = 2e-15
+ else:
+ real_abs_err = 5e-323
+
+ if not (almostEqualF(expected.real, actual.real,
+ abs_err = real_abs_err) and
+ almostEqualF(expected.imag, actual.imag)):
+ error_message = (
+ "%s: %s(complex(%r, %r))\n" % (id, fn, ar, ai) +
+ "Expected: complex(%r, %r)\n" %
+ (expected.real, expected.imag) +
+ "Received: complex(%r, %r)\n" %
+ (actual.real, actual.imag) +
+ "Received value insufficiently close to expected value.")
+ self.fail(error_message)
+
+ def assertCISEqual(self, a, b):
+ eps = 1E-7
+ if abs(a[0] - b[0]) > eps or abs(a[1] - b[1]) > eps:
+ self.fail((a ,b))
+
+ def test_polar(self):
+ self.assertCISEqual(polar(0), (0., 0.))
+ self.assertCISEqual(polar(1.), (1., 0.))
+ self.assertCISEqual(polar(-1.), (1., pi))
+ self.assertCISEqual(polar(1j), (1., pi/2))
+ self.assertCISEqual(polar(-1j), (1., -pi/2))
+
+ def test_phase(self):
+ self.assertAlmostEqual(phase(0), 0.)
+ self.assertAlmostEqual(phase(1.), 0.)
+ self.assertAlmostEqual(phase(-1.), pi)
+ self.assertAlmostEqual(phase(-1.+1E-300j), pi)
+ self.assertAlmostEqual(phase(-1.-1E-300j), -pi)
+ self.assertAlmostEqual(phase(1j), pi/2)
+ self.assertAlmostEqual(phase(-1j), -pi/2)
+
+ # zeros
+ self.assertEqual(phase(complex(0.0, 0.0)), 0.0)
+ self.assertEqual(phase(complex(0.0, -0.0)), -0.0)
+ self.assertEqual(phase(complex(-0.0, 0.0)), pi)
+ self.assertEqual(phase(complex(-0.0, -0.0)), -pi)
+
+ # infinities
+ self.assertAlmostEqual(phase(complex(-INF, -0.0)), -pi)
+ self.assertAlmostEqual(phase(complex(-INF, -2.3)), -pi)
+ self.assertAlmostEqual(phase(complex(-INF, -INF)), -0.75*pi)
+ self.assertAlmostEqual(phase(complex(-2.3, -INF)), -pi/2)
+ self.assertAlmostEqual(phase(complex(-0.0, -INF)), -pi/2)
+ self.assertAlmostEqual(phase(complex(0.0, -INF)), -pi/2)
+ self.assertAlmostEqual(phase(complex(2.3, -INF)), -pi/2)
+ self.assertAlmostEqual(phase(complex(INF, -INF)), -pi/4)
+ self.assertEqual(phase(complex(INF, -2.3)), -0.0)
+ self.assertEqual(phase(complex(INF, -0.0)), -0.0)
+ self.assertEqual(phase(complex(INF, 0.0)), 0.0)
+ self.assertEqual(phase(complex(INF, 2.3)), 0.0)
+ self.assertAlmostEqual(phase(complex(INF, INF)), pi/4)
+ self.assertAlmostEqual(phase(complex(2.3, INF)), pi/2)
+ self.assertAlmostEqual(phase(complex(0.0, INF)), pi/2)
+ self.assertAlmostEqual(phase(complex(-0.0, INF)), pi/2)
+ self.assertAlmostEqual(phase(complex(-2.3, INF)), pi/2)
+ self.assertAlmostEqual(phase(complex(-INF, INF)), 0.75*pi)
+ self.assertAlmostEqual(phase(complex(-INF, 2.3)), pi)
+ self.assertAlmostEqual(phase(complex(-INF, 0.0)), pi)
+
+ # real or imaginary part NaN
+ for z in complex_nans:
+ self.assert_(math.isnan(phase(z)))
+
+ def test_abs(self):
+ # zeros
+ for z in complex_zeros:
+ self.assertEqual(abs(z), 0.0)
+
+ # infinities
+ for z in complex_infinities:
+ self.assertEqual(abs(z), INF)
+
+ # real or imaginary part NaN
+ self.assertEqual(abs(complex(NAN, -INF)), INF)
+ self.assert_(math.isnan(abs(complex(NAN, -2.3))))
+ self.assert_(math.isnan(abs(complex(NAN, -0.0))))
+ self.assert_(math.isnan(abs(complex(NAN, 0.0))))
+ self.assert_(math.isnan(abs(complex(NAN, 2.3))))
+ self.assertEqual(abs(complex(NAN, INF)), INF)
+ self.assertEqual(abs(complex(-INF, NAN)), INF)
+ self.assert_(math.isnan(abs(complex(-2.3, NAN))))
+ self.assert_(math.isnan(abs(complex(-0.0, NAN))))
+ self.assert_(math.isnan(abs(complex(0.0, NAN))))
+ self.assert_(math.isnan(abs(complex(2.3, NAN))))
+ self.assertEqual(abs(complex(INF, NAN)), INF)
+ self.assert_(math.isnan(abs(complex(NAN, NAN))))
+
+ # result overflows
+ if float.__getformat__("double").startswith("IEEE"):
+ self.assertRaises(OverflowError, abs, complex(1.4e308, 1.4e308))
+
+ def assertCEqual(self, a, b):
+ eps = 1E-7
+ if abs(a.real - b[0]) > eps or abs(a.imag - b[1]) > eps:
+ self.fail((a ,b))
+
+ def test_rect(self):
+ self.assertCEqual(rect(0, 0), (0, 0))
+ self.assertCEqual(rect(1, 0), (1., 0))
+ self.assertCEqual(rect(1, -pi), (-1., 0))
+ self.assertCEqual(rect(1, pi/2), (0, 1.))
+ self.assertCEqual(rect(1, -pi/2), (0, -1.))
+
+ def test_isnan(self):
+ self.failIf(cmath.isnan(1))
+ self.failIf(cmath.isnan(1j))
+ self.failIf(cmath.isnan(INF))
+ self.assert_(cmath.isnan(NAN))
+ self.assert_(cmath.isnan(complex(NAN, 0)))
+ self.assert_(cmath.isnan(complex(0, NAN)))
+ self.assert_(cmath.isnan(complex(NAN, NAN)))
+ self.assert_(cmath.isnan(complex(NAN, INF)))
+ self.assert_(cmath.isnan(complex(INF, NAN)))
+
+ def test_isinf(self):
+ self.failIf(cmath.isinf(1))
+ self.failIf(cmath.isinf(1j))
+ self.failIf(cmath.isinf(NAN))
+ self.assert_(cmath.isinf(INF))
+ self.assert_(cmath.isinf(complex(INF, 0)))
+ self.assert_(cmath.isinf(complex(0, INF)))
+ self.assert_(cmath.isinf(complex(INF, INF)))
+ self.assert_(cmath.isinf(complex(NAN, INF)))
+ self.assert_(cmath.isinf(complex(INF, NAN)))
+
def test_main():
run_unittest(CMathTests)
Modified: python/branches/py3k/Lib/test/test_float.py
==============================================================================
--- python/branches/py3k/Lib/test/test_float.py (original)
+++ python/branches/py3k/Lib/test/test_float.py Sat Apr 19 02:31:39 2008
@@ -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 @@
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,
Modified: python/branches/py3k/Lib/test/test_math.py
==============================================================================
--- python/branches/py3k/Lib/test/test_math.py (original)
+++ python/branches/py3k/Lib/test/test_math.py Sat Apr 19 02:31:39 2008
@@ -4,9 +4,45 @@
from test.test_support import run_unittest, verbose
import unittest
import math
+import os
+import sys
-seps='1e-05'
-eps = eval(seps)
+eps = 1E-05
+NAN = float('nan')
+INF = float('inf')
+NINF = float('-inf')
+
+# locate file with test values
+if __name__ == '__main__':
+ file = sys.argv[0]
+else:
+ file = __file__
+test_dir = os.path.dirname(file) or os.curdir
+test_file = os.path.join(test_dir, 'cmath_testcases.txt')
+
+def parse_testfile(fname):
+ """Parse a file with test values
+
+ Empty lines or lines starting with -- are ignored
+ yields id, fn, arg_real, arg_imag, exp_real, exp_imag
+ """
+ with open(fname) as fp:
+ for line in fp:
+ # skip comment lines and blank lines
+ if line.startswith('--') or not line.strip():
+ continue
+
+ lhs, rhs = line.split('->')
+ id, fn, arg_real, arg_imag = lhs.split()
+ rhs_pieces = rhs.split()
+ exp_real, exp_imag = rhs_pieces[0], rhs_pieces[1]
+ flags = rhs_pieces[2:]
+
+ yield (id, fn,
+ float(arg_real), float(arg_imag),
+ float(exp_real), float(exp_imag),
+ flags
+ )
class MathTests(unittest.TestCase):
@@ -28,18 +64,57 @@
self.ftest('acos(-1)', math.acos(-1), math.pi)
self.ftest('acos(0)', math.acos(0), math.pi/2)
self.ftest('acos(1)', math.acos(1), 0)
+ self.assertRaises(ValueError, math.acos, INF)
+ self.assertRaises(ValueError, math.acos, NINF)
+ self.assert_(math.isnan(math.acos(NAN)))
+
+ def testAcosh(self):
+ self.assertRaises(TypeError, math.acosh)
+ self.ftest('acosh(1)', math.acosh(1), 0)
+ self.ftest('acosh(2)', math.acosh(2), 1.3169578969248168)
+ self.assertRaises(ValueError, math.acosh, 0)
+ self.assertRaises(ValueError, math.acosh, -1)
+ self.assertEquals(math.acosh(INF), INF)
+ self.assertRaises(ValueError, math.acosh, NINF)
+ self.assert_(math.isnan(math.acosh(NAN)))
def testAsin(self):
self.assertRaises(TypeError, math.asin)
self.ftest('asin(-1)', math.asin(-1), -math.pi/2)
self.ftest('asin(0)', math.asin(0), 0)
self.ftest('asin(1)', math.asin(1), math.pi/2)
+ self.assertRaises(ValueError, math.asin, INF)
+ self.assertRaises(ValueError, math.asin, NINF)
+ self.assert_(math.isnan(math.asin(NAN)))
+
+ def testAsinh(self):
+ self.assertRaises(TypeError, math.asinh)
+ self.ftest('asinh(0)', math.asinh(0), 0)
+ self.ftest('asinh(1)', math.asinh(1), 0.88137358701954305)
+ self.ftest('asinh(-1)', math.asinh(-1), -0.88137358701954305)
+ self.assertEquals(math.asinh(INF), INF)
+ self.assertEquals(math.asinh(NINF), NINF)
+ self.assert_(math.isnan(math.asinh(NAN)))
def testAtan(self):
self.assertRaises(TypeError, math.atan)
self.ftest('atan(-1)', math.atan(-1), -math.pi/4)
self.ftest('atan(0)', math.atan(0), 0)
self.ftest('atan(1)', math.atan(1), math.pi/4)
+ self.ftest('atan(inf)', math.atan(INF), math.pi/2)
+ self.ftest('atan(-inf)', math.atan(-INF), -math.pi/2)
+ self.assert_(math.isnan(math.atan(NAN)))
+
+ def testAtanh(self):
+ self.assertRaises(TypeError, math.atan)
+ self.ftest('atanh(0)', math.atanh(0), 0)
+ self.ftest('atanh(0.5)', math.atanh(0.5), 0.54930614433405489)
+ self.ftest('atanh(-0.5)', math.atanh(-0.5), -0.54930614433405489)
+ self.assertRaises(ValueError, math.atanh, 1)
+ self.assertRaises(ValueError, math.atanh, -1)
+ self.assertRaises(ValueError, math.atanh, INF)
+ self.assertRaises(ValueError, math.atanh, NINF)
+ self.assert_(math.isnan(math.atanh(NAN)))
def testAtan2(self):
self.assertRaises(TypeError, math.atan2)
@@ -58,6 +133,9 @@
self.ftest('ceil(-0.5)', math.ceil(-0.5), 0)
self.ftest('ceil(-1.0)', math.ceil(-1.0), -1)
self.ftest('ceil(-1.5)', math.ceil(-1.5), -1)
+ #self.assertEquals(math.ceil(INF), INF)
+ #self.assertEquals(math.ceil(NINF), NINF)
+ #self.assert_(math.isnan(math.ceil(NAN)))
class TestCeil:
def __ceil__(self):
@@ -72,17 +150,55 @@
self.assertRaises(TypeError, math.ceil, t)
self.assertRaises(TypeError, math.ceil, t, 0)
+ if float.__getformat__("double").startswith("IEEE"):
+ def testCopysign(self):
+ self.assertRaises(TypeError, math.copysign)
+ # copysign should let us distinguish signs of zeros
+ self.assertEquals(copysign(1., 0.), 1.)
+ self.assertEquals(copysign(1., -0.), -1.)
+ self.assertEquals(copysign(INF, 0.), INF)
+ self.assertEquals(copysign(INF, -0.), NINF)
+ self.assertEquals(copysign(NINF, 0.), INF)
+ self.assertEquals(copysign(NINF, -0.), NINF)
+ # and of infinities
+ self.assertEquals(copysign(1., INF), 1.)
+ self.assertEquals(copysign(1., NINF), -1.)
+ self.assertEquals(copysign(INF, INF), INF)
+ self.assertEquals(copysign(INF, NINF), NINF)
+ self.assertEquals(copysign(NINF, INF), INF)
+ self.assertEquals(copysign(NINF, NINF), NINF)
+ self.assert_(math.isnan(copysign(NAN, 1.)))
+ self.assert_(math.isnan(copysign(NAN, INF)))
+ self.assert_(math.isnan(copysign(NAN, NINF)))
+ self.assert_(math.isnan(copysign(NAN, NAN)))
+ # copysign(INF, NAN) may be INF or it may be NINF, since
+ # we don't know whether the sign bit of NAN is set on any
+ # given platform.
+ self.assert_(math.isinf(copysign(INF, NAN)))
+ # similarly, copysign(2., NAN) could be 2. or -2.
+ self.assertEquals(abs(copysign(2., NAN)), 2.)
+
def testCos(self):
self.assertRaises(TypeError, math.cos)
self.ftest('cos(-pi/2)', math.cos(-math.pi/2), 0)
self.ftest('cos(0)', math.cos(0), 1)
self.ftest('cos(pi/2)', math.cos(math.pi/2), 0)
self.ftest('cos(pi)', math.cos(math.pi), -1)
+ try:
+ self.assert_(math.isnan(math.cos(INF)))
+ self.assert_(math.isnan(math.cos(NINF)))
+ except ValueError:
+ self.assertRaises(ValueError, math.cos, INF)
+ self.assertRaises(ValueError, math.cos, NINF)
+ self.assert_(math.isnan(math.cos(NAN)))
def testCosh(self):
self.assertRaises(TypeError, math.cosh)
self.ftest('cosh(0)', math.cosh(0), 1)
self.ftest('cosh(2)-2*cosh(1)**2', math.cosh(2)-2*math.cosh(1)**2, -1) # Thanks to Lambert
+ self.assertEquals(math.cosh(INF), INF)
+ self.assertEquals(math.cosh(NINF), INF)
+ self.assert_(math.isnan(math.cosh(NAN)))
def testDegrees(self):
self.assertRaises(TypeError, math.degrees)
@@ -95,6 +211,9 @@
self.ftest('exp(-1)', math.exp(-1), 1/math.e)
self.ftest('exp(0)', math.exp(0), 1)
self.ftest('exp(1)', math.exp(1), math.e)
+ self.assertEquals(math.exp(INF), INF)
+ self.assertEquals(math.exp(NINF), 0.)
+ self.assert_(math.isnan(math.exp(NAN)))
def testFabs(self):
self.assertRaises(TypeError, math.fabs)
@@ -115,6 +234,9 @@
# This fails on some platforms - so check it here
self.ftest('floor(1.23e167)', math.floor(1.23e167), 1.23e167)
self.ftest('floor(-1.23e167)', math.floor(-1.23e167), -1.23e167)
+ #self.assertEquals(math.ceil(INF), INF)
+ #self.assertEquals(math.ceil(NINF), NINF)
+ #self.assert_(math.isnan(math.floor(NAN)))
class TestFloor:
def __floor__(self):
@@ -137,6 +259,19 @@
self.ftest('fmod(-10,1)', math.fmod(-10,1), 0)
self.ftest('fmod(-10,0.5)', math.fmod(-10,0.5), 0)
self.ftest('fmod(-10,1.5)', math.fmod(-10,1.5), -1)
+ self.assert_(math.isnan(math.fmod(NAN, 1.)))
+ self.assert_(math.isnan(math.fmod(1., NAN)))
+ self.assert_(math.isnan(math.fmod(NAN, NAN)))
+ self.assertRaises(ValueError, math.fmod, 1., 0.)
+ self.assertRaises(ValueError, math.fmod, INF, 1.)
+ self.assertRaises(ValueError, math.fmod, NINF, 1.)
+ self.assertRaises(ValueError, math.fmod, INF, 0.)
+ self.assertEquals(math.fmod(3.0, INF), 3.0)
+ self.assertEquals(math.fmod(-3.0, INF), -3.0)
+ self.assertEquals(math.fmod(3.0, NINF), 3.0)
+ self.assertEquals(math.fmod(-3.0, NINF), -3.0)
+ self.assertEquals(math.fmod(0.0, 3.0), 0.0)
+ self.assertEquals(math.fmod(0.0, NINF), 0.0)
def testFrexp(self):
self.assertRaises(TypeError, math.frexp)
@@ -152,10 +287,20 @@
testfrexp('frexp(1)', math.frexp(1), (0.5, 1))
testfrexp('frexp(2)', math.frexp(2), (0.5, 2))
+ self.assertEquals(math.frexp(INF)[0], INF)
+ self.assertEquals(math.frexp(NINF)[0], NINF)
+ self.assert_(math.isnan(math.frexp(NAN)[0]))
+
def testHypot(self):
self.assertRaises(TypeError, math.hypot)
self.ftest('hypot(0,0)', math.hypot(0,0), 0)
self.ftest('hypot(3,4)', math.hypot(3,4), 5)
+ self.assertEqual(math.hypot(NAN, INF), INF)
+ self.assertEqual(math.hypot(INF, NAN), INF)
+ self.assertEqual(math.hypot(NAN, NINF), INF)
+ self.assertEqual(math.hypot(NINF, NAN), INF)
+ self.assert_(math.isnan(math.hypot(1.0, NAN)))
+ self.assert_(math.isnan(math.hypot(NAN, -2.0)))
def testLdexp(self):
self.assertRaises(TypeError, math.ldexp)
@@ -163,6 +308,13 @@
self.ftest('ldexp(1,1)', math.ldexp(1,1), 2)
self.ftest('ldexp(1,-1)', math.ldexp(1,-1), 0.5)
self.ftest('ldexp(-1,1)', math.ldexp(-1,1), -2)
+ self.assertRaises(OverflowError, math.ldexp, 1., 1000000)
+ self.assertRaises(OverflowError, math.ldexp, -1., 1000000)
+ self.assertEquals(math.ldexp(1., -1000000), 0.)
+ self.assertEquals(math.ldexp(-1., -1000000), -0.)
+ self.assertEquals(math.ldexp(INF, 30), INF)
+ self.assertEquals(math.ldexp(NINF, -213), NINF)
+ self.assert_(math.isnan(math.ldexp(NAN, 0)))
def testLog(self):
self.assertRaises(TypeError, math.log)
@@ -172,12 +324,31 @@
self.ftest('log(32,2)', math.log(32,2), 5)
self.ftest('log(10**40, 10)', math.log(10**40, 10), 40)
self.ftest('log(10**40, 10**20)', math.log(10**40, 10**20), 2)
+ self.assertEquals(math.log(INF), INF)
+ self.assertRaises(ValueError, math.log, NINF)
+ self.assert_(math.isnan(math.log(NAN)))
+
+ def testLog1p(self):
+ self.assertRaises(TypeError, math.log1p)
+ self.ftest('log1p(1/e -1)', math.log1p(1/math.e-1), -1)
+ self.ftest('log1p(0)', math.log1p(0), 0)
+ self.ftest('log1p(e-1)', math.log1p(math.e-1), 1)
+ self.ftest('log1p(1)', math.log1p(1), math.log(2))
+ self.assertEquals(math.log1p(INF), INF)
+ self.assertRaises(ValueError, math.log1p, NINF)
+ self.assert_(math.isnan(math.log1p(NAN)))
+ n= 2**90
+ self.assertAlmostEquals(math.log1p(n), 62.383246250395075)
+ self.assertAlmostEquals(math.log1p(n), math.log1p(float(n)))
def testLog10(self):
self.assertRaises(TypeError, math.log10)
self.ftest('log10(0.1)', math.log10(0.1), -1)
self.ftest('log10(1)', math.log10(1), 0)
self.ftest('log10(10)', math.log10(10), 1)
+ self.assertEquals(math.log(INF), INF)
+ self.assertRaises(ValueError, math.log10, NINF)
+ self.assert_(math.isnan(math.log10(NAN)))
def testModf(self):
self.assertRaises(TypeError, math.modf)
@@ -191,12 +362,35 @@
testmodf('modf(1.5)', math.modf(1.5), (0.5, 1.0))
testmodf('modf(-1.5)', math.modf(-1.5), (-0.5, -1.0))
+ self.assertEquals(math.modf(INF), (0.0, INF))
+ self.assertEquals(math.modf(NINF), (-0.0, NINF))
+
+ modf_nan = math.modf(NAN)
+ self.assert_(math.isnan(modf_nan[0]))
+ self.assert_(math.isnan(modf_nan[1]))
+
def testPow(self):
self.assertRaises(TypeError, math.pow)
self.ftest('pow(0,1)', math.pow(0,1), 0)
self.ftest('pow(1,0)', math.pow(1,0), 1)
self.ftest('pow(2,1)', math.pow(2,1), 2)
self.ftest('pow(2,-1)', math.pow(2,-1), 0.5)
+ self.assertEqual(math.pow(INF, 1), INF)
+ self.assertEqual(math.pow(NINF, 1), NINF)
+ self.assertEqual((math.pow(1, INF)), 1.)
+ self.assertEqual((math.pow(1, NINF)), 1.)
+ self.assert_(math.isnan(math.pow(NAN, 1)))
+ self.assert_(math.isnan(math.pow(2, NAN)))
+ self.assert_(math.isnan(math.pow(0, NAN)))
+ self.assertEqual(math.pow(1, NAN), 1)
+ self.assertEqual(1**NAN, 1)
+ self.assertEqual(1**INF, 1)
+ self.assertEqual(1**NINF, 1)
+ self.assertEqual(1**0, 1)
+ self.assertEqual(1.**NAN, 1)
+ self.assertEqual(1.**INF, 1)
+ self.assertEqual(1.**NINF, 1)
+ self.assertEqual(1.**0, 1)
def testRadians(self):
self.assertRaises(TypeError, math.radians)
@@ -209,29 +403,52 @@
self.ftest('sin(0)', math.sin(0), 0)
self.ftest('sin(pi/2)', math.sin(math.pi/2), 1)
self.ftest('sin(-pi/2)', math.sin(-math.pi/2), -1)
+ try:
+ self.assert_(math.isnan(math.sin(INF)))
+ self.assert_(math.isnan(math.sin(NINF)))
+ except ValueError:
+ self.assertRaises(ValueError, math.sin, INF)
+ self.assertRaises(ValueError, math.sin, NINF)
+ self.assert_(math.isnan(math.sin(NAN)))
def testSinh(self):
self.assertRaises(TypeError, math.sinh)
self.ftest('sinh(0)', math.sinh(0), 0)
self.ftest('sinh(1)**2-cosh(1)**2', math.sinh(1)**2-math.cosh(1)**2, -1)
self.ftest('sinh(1)+sinh(-1)', math.sinh(1)+math.sinh(-1), 0)
+ self.assertEquals(math.sinh(INF), INF)
+ self.assertEquals(math.sinh(-INF), -INF)
+ self.assert_(math.isnan(math.sinh(NAN)))
def testSqrt(self):
self.assertRaises(TypeError, math.sqrt)
self.ftest('sqrt(0)', math.sqrt(0), 0)
self.ftest('sqrt(1)', math.sqrt(1), 1)
self.ftest('sqrt(4)', math.sqrt(4), 2)
+ self.assertEquals(math.sqrt(INF), INF)
+ self.assertRaises(ValueError, math.sqrt, NINF)
+ self.assert_(math.isnan(math.sqrt(NAN)))
def testTan(self):
self.assertRaises(TypeError, math.tan)
self.ftest('tan(0)', math.tan(0), 0)
self.ftest('tan(pi/4)', math.tan(math.pi/4), 1)
self.ftest('tan(-pi/4)', math.tan(-math.pi/4), -1)
+ try:
+ self.assert_(math.isnan(math.tan(INF)))
+ self.assert_(math.isnan(math.tan(NINF)))
+ except:
+ self.assertRaises(ValueError, math.tan, INF)
+ self.assertRaises(ValueError, math.tan, NINF)
+ self.assert_(math.isnan(math.tan(NAN)))
def testTanh(self):
self.assertRaises(TypeError, math.tanh)
self.ftest('tanh(0)', math.tanh(0), 0)
self.ftest('tanh(1)+tanh(-1)', math.tanh(1)+math.tanh(-1), 0)
+ self.ftest('tanh(inf)', math.tanh(INF), 1)
+ self.ftest('tanh(-inf)', math.tanh(NINF), -1)
+ self.assert_(math.isnan(math.tanh(NAN)))
def test_trunc(self):
self.assertEqual(math.trunc(1), 1)
@@ -326,9 +543,27 @@
else:
self.fail("sqrt(-1) didn't raise ValueError")
+ def test_testfile(self):
+ if not float.__getformat__("double").startswith("IEEE"):
+ return
+ for id, fn, ar, ai, er, ei, flags in parse_testfile(test_file):
+ # Skip if either the input or result is complex, or if
+ # flags is nonempty
+ if ai != 0. or ei != 0. or flags:
+ continue
+ if fn in ['rect', 'polar']:
+ # no real versions of rect, polar
+ continue
+ func = getattr(math, fn)
+ result = func(ar)
+ self.ftest("%s:%s(%r)" % (id, fn, ar), result, er)
def test_main():
- run_unittest(MathTests)
+ from doctest import DocFileSuite
+ suite = unittest.TestSuite()
+ suite.addTest(unittest.makeSuite(MathTests))
+ suite.addTest(DocFileSuite("ieee754.txt"))
+ run_unittest(suite)
if __name__ == '__main__':
test_main()
Modified: python/branches/py3k/Makefile.pre.in
==============================================================================
--- python/branches/py3k/Makefile.pre.in (original)
+++ python/branches/py3k/Makefile.pre.in Sat Apr 19 02:31:39 2008
@@ -276,6 +276,7 @@
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 @@
Include/pydebug.h \
Include/pyerrors.h \
Include/pyfpe.h \
+ Include/pymath.h \
Include/pygetopt.h \
Include/pymem.h \
Include/pyport.h \
Modified: python/branches/py3k/Modules/cmathmodule.c
==============================================================================
--- python/branches/py3k/Modules/cmathmodule.c (original)
+++ python/branches/py3k/Modules/cmathmodule.c Sat Apr 19 02:31:39 2008
@@ -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
+
+#ifndef M_LN10
+#define M_LN10 (2.302585092994045684) /* natural log of 10 */
+#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};
+/*
+ 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));
}
Modified: python/branches/py3k/Modules/mathmodule.c
==============================================================================
--- python/branches/py3k/Modules/mathmodule.c (original)
+++ python/branches/py3k/Modules/mathmodule.c Sat Apr 19 02:31:39 2008
@@ -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))
+math_1(PyObject *arg, double (*func) (double), int can_overflow)
{
- return math_1_to_whatever(arg, func, PyFloat_FromDouble);
+ return math_1_to_whatever(arg, func, PyFloat_FromDouble, can_overflow);
}
static PyObject *
-math_1_to_int(PyObject *arg, double (*func) (double))
+math_1_to_int(PyObject *arg, double (*func) (double), int can_overflow)
{
- 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,
+FUNC2(copysign, copysign,
+ "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,
+FUNC1(cosh, cosh, 1,
"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,
+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);
-
- if (!(v = PyFloat_FromDouble(Py_MATH_E)))
- goto finally;
- if (PyDict_SetItemString(d, "e", 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));
- finally:
+ finally:
return;
}
Modified: python/branches/py3k/Objects/complexobject.c
==============================================================================
--- python/branches/py3k/Objects/complexobject.c (original)
+++ python/branches/py3k/Objects/complexobject.c Sat Apr 19 02:31:39 2008
@@ -187,6 +187,38 @@
}
+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 @@
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 @@
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 @@
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 */
};
Deleted: python/branches/py3k/Objects/doubledigits.c
==============================================================================
--- python/branches/py3k/Objects/doubledigits.c Sat Apr 19 02:31:39 2008
+++ (empty file)
@@ -1,601 +0,0 @@
-/* Free-format floating point printer
- *
- * Based on "Floating-Point Printer Sample Code", by Robert G. Burger,
- * http://www.cs.indiana.edu/~burger/fp/index.html
- */
-
-#include "Python.h"
-
-#if defined(__alpha) || defined(__i386) || defined(_M_IX86) || defined(_M_X64) || defined(_M_IA64)
-#define LITTLE_ENDIAN_IEEE_DOUBLE
-#elif !(defined(__ppc__) || defined(sparc) || defined(__sgi) || defined(_IBMR2) || defined(hpux))
-#error unknown machine type
-#endif
-
-#if defined(_M_IX86)
-#define UNSIGNED64 unsigned __int64
-#elif defined(__alpha)
-#define UNSIGNED64 unsigned long
-#else
-#define UNSIGNED64 unsigned long long
-#endif
-
-#ifndef U32
-#define U32 unsigned int
-#endif
-
-/* exponent stored + 1024, hidden bit to left of decimal point */
-#define bias 1023
-#define bitstoright 52
-#define m1mask 0xf
-#define hidden_bit 0x100000
-#ifdef LITTLE_ENDIAN_IEEE_DOUBLE
-struct dblflt {
- unsigned int m4: 16;
- unsigned int m3: 16;
- unsigned int m2: 16;
- unsigned int m1: 4;
- unsigned int e: 11;
- unsigned int s: 1;
-};
-#else
-/* Big Endian IEEE Double Floats */
-struct dblflt {
- unsigned int s: 1;
- unsigned int e: 11;
- unsigned int m1: 4;
- unsigned int m2: 16;
- unsigned int m3: 16;
- unsigned int m4: 16;
-};
-#endif
-#define float_radix 2.147483648e9
-
-
-typedef UNSIGNED64 Bigit;
-#define BIGSIZE 24
-#define MIN_E -1074
-#define MAX_FIVE 325
-#define B_P1 ((Bigit)1 << 52)
-
-typedef struct {
- int l;
- Bigit d[BIGSIZE];
-} Bignum;
-
-static Bignum R, S, MP, MM, five[MAX_FIVE];
-static Bignum S2, S3, S4, S5, S6, S7, S8, S9;
-static int ruf, k, s_n, use_mp, qr_shift, sl, slr;
-
-static void mul10(Bignum *x);
-static void big_short_mul(Bignum *x, Bigit y, Bignum *z);
-/*
-static void print_big(Bignum *x);
-*/
-static int estimate(int n);
-static void one_shift_left(int y, Bignum *z);
-static void short_shift_left(Bigit x, int y, Bignum *z);
-static void big_shift_left(Bignum *x, int y, Bignum *z);
-static int big_comp(Bignum *x, Bignum *y);
-static int sub_big(Bignum *x, Bignum *y, Bignum *z);
-static void add_big(Bignum *x, Bignum *y, Bignum *z);
-static int add_cmp(void);
-static int qr(void);
-
-/*static int _PyFloat_Digits(char *buf, double v, int *signum);*/
-/*static void _PyFloat_DigitsInit(void);*/
-
-#define ADD(x, y, z, k) {\
- Bigit x_add, z_add;\
- x_add = (x);\
- if ((k))\
- z_add = x_add + (y) + 1, (k) = (z_add <= x_add);\
- else\
- z_add = x_add + (y), (k) = (z_add < x_add);\
- (z) = z_add;\
-}
-
-#define SUB(x, y, z, b) {\
- Bigit x_sub, y_sub;\
- x_sub = (x); y_sub = (y);\
- if ((b))\
- (z) = x_sub - y_sub - 1, b = (y_sub >= x_sub);\
- else\
- (z) = x_sub - y_sub, b = (y_sub > x_sub);\
-}
-
-#define MUL(x, y, z, k) {\
- Bigit x_mul, low, high;\
- x_mul = (x);\
- low = (x_mul & 0xffffffff) * (y) + (k);\
- high = (x_mul >> 32) * (y) + (low >> 32);\
- (k) = high >> 32;\
- (z) = (low & 0xffffffff) | (high << 32);\
-}
-
-#define SLL(x, y, z, k) {\
- Bigit x_sll = (x);\
- (z) = (x_sll << (y)) | (k);\
- (k) = x_sll >> (64 - (y));\
-}
-
-static void
-mul10(Bignum *x)
-{
- int i, l;
- Bigit *p, k;
-
- l = x->l;
- for (i = l, p = &x->d[0], k = 0; i >= 0; i--)
- MUL(*p, 10, *p++, k);
- if (k != 0)
- *p = k, x->l = l+1;
-}
-
-static void
-big_short_mul(Bignum *x, Bigit y, Bignum *z)
-{
- int i, xl, zl;
- Bigit *xp, *zp, k;
- U32 high, low;
-
- xl = x->l;
- xp = &x->d[0];
- zl = xl;
- zp = &z->d[0];
- high = y >> 32;
- low = y & 0xffffffff;
- for (i = xl, k = 0; i >= 0; i--, xp++, zp++) {
- Bigit xlow, xhigh, z0, t, c, z1;
- xlow = *xp & 0xffffffff;
- xhigh = *xp >> 32;
- z0 = (xlow * low) + k; /* Cout is (z0 < k) */
- t = xhigh * low;
- z1 = (xlow * high) + t;
- c = (z1 < t);
- t = z0 >> 32;
- z1 += t;
- c += (z1 < t);
- *zp = (z1 << 32) | (z0 & 0xffffffff);
- k = (xhigh * high) + (c << 32) + (z1 >> 32) + (z0 < k);
- }
- if (k != 0)
- *zp = k, zl++;
- z->l = zl;
-}
-
-/*
-static void
-print_big(Bignum *x)
-{
- int i;
- Bigit *p;
-
- printf("#x");
- i = x->l;
- p = &x->d[i];
- for (p = &x->d[i]; i >= 0; i--) {
- Bigit b = *p--;
- printf("%08x%08x", (int)(b >> 32), (int)(b & 0xffffffff));
- }
-}
-*/
-
-static int
-estimate(int n)
-{
- if (n < 0)
- return (int)(n*0.3010299956639812);
- else
- return 1+(int)(n*0.3010299956639811);
-}
-
-static void
-one_shift_left(int y, Bignum *z)
-{
- int n, m, i;
- Bigit *zp;
-
- n = y / 64;
- m = y % 64;
- zp = &z->d[0];
- for (i = n; i > 0; i--) *zp++ = 0;
- *zp = (Bigit)1 << m;
- z->l = n;
-}
-
-static void
-short_shift_left(Bigit x, int y, Bignum *z)
-{
- int n, m, i, zl;
- Bigit *zp;
-
- n = y / 64;
- m = y % 64;
- zl = n;
- zp = &(z->d[0]);
- for (i = n; i > 0; i--) *zp++ = 0;
- if (m == 0)
- *zp = x;
- else {
- Bigit high = x >> (64 - m);
- *zp = x << m;
- if (high != 0)
- *++zp = high, zl++;
- }
- z->l = zl;
-}
-
-static void
-big_shift_left(Bignum *x, int y, Bignum *z)
-{
- int n, m, i, xl, zl;
- Bigit *xp, *zp, k;
-
- n = y / 64;
- m = y % 64;
- xl = x->l;
- xp = &(x->d[0]);
- zl = xl + n;
- zp = &(z->d[0]);
- for (i = n; i > 0; i--) *zp++ = 0;
- if (m == 0)
- for (i = xl; i >= 0; i--) *zp++ = *xp++;
- else {
- for (i = xl, k = 0; i >= 0; i--)
- SLL(*xp++, m, *zp++, k);
- if (k != 0)
- *zp = k, zl++;
- }
- z->l = zl;
-}
-
-
-static int
-big_comp(Bignum *x, Bignum *y)
-{
- int i, xl, yl;
- Bigit *xp, *yp;
-
- xl = x->l;
- yl = y->l;
- if (xl > yl) return 1;
- if (xl < yl) return -1;
- xp = &x->d[xl];
- yp = &y->d[xl];
- for (i = xl; i >= 0; i--, xp--, yp--) {
- Bigit a = *xp;
- Bigit b = *yp;
-
- if (a > b) return 1;
- else if (a < b) return -1;
- }
- return 0;
-}
-
-static int
-sub_big(Bignum *x, Bignum *y, Bignum *z)
-{
- int xl, yl, zl, b, i;
- Bigit *xp, *yp, *zp;
-
- xl = x->l;
- yl = y->l;
- if (yl > xl) return 1;
- xp = &x->d[0];
- yp = &y->d[0];
- zp = &z->d[0];
-
- for (i = yl, b = 0; i >= 0; i--)
- SUB(*xp++, *yp++, *zp++, b);
- for (i = xl-yl; b && i > 0; i--) {
- Bigit x_sub;
- x_sub = *xp++;
- *zp++ = x_sub - 1;
- b = (x_sub == 0);
- }
- for (; i > 0; i--) *zp++ = *xp++;
- if (b) return 1;
- zl = xl;
- while (*--zp == 0) zl--;
- z->l = zl;
- return 0;
-}
-
-static void
-add_big(Bignum *x, Bignum *y, Bignum *z)
-{
- int xl, yl, k, i;
- Bigit *xp, *yp, *zp;
-
- xl = x->l;
- yl = y->l;
- if (yl > xl) {
- int tl;
- Bignum *tn;
- tl = xl; xl = yl; yl = tl;
- tn = x; x = y; y = tn;
- }
-
- xp = &x->d[0];
- yp = &y->d[0];
- zp = &z->d[0];
-
- for (i = yl, k = 0; i >= 0; i--)
- ADD(*xp++, *yp++, *zp++, k);
- for (i = xl-yl; k && i > 0; i--) {
- Bigit z_add;
- z_add = *xp++ + 1;
- k = (z_add == 0);
- *zp++ = z_add;
- }
- for (; i > 0; i--) *zp++ = *xp++;
- if (k)
- *zp = 1, z->l = xl+1;
- else
- z->l = xl;
-}
-
-static int
-add_cmp()
-{
- int rl, ml, sl, suml;
- static Bignum sum;
-
- rl = R.l;
- ml = (use_mp ? MP.l : MM.l);
- sl = S.l;
-
- suml = rl >= ml ? rl : ml;
- if ((sl > suml+1) || ((sl == suml+1) && (S.d[sl] > 1))) return -1;
- if (sl < suml) return 1;
-
- add_big(&R, (use_mp ? &MP : &MM), &sum);
- return big_comp(&sum, &S);
-}
-
-static int
-qr()
-{
- if (big_comp(&R, &S5) < 0)
- if (big_comp(&R, &S2) < 0)
- if (big_comp(&R, &S) < 0)
- return 0;
- else {
- sub_big(&R, &S, &R);
- return 1;
- }
- else if (big_comp(&R, &S3) < 0) {
- sub_big(&R, &S2, &R);
- return 2;
- }
- else if (big_comp(&R, &S4) < 0) {
- sub_big(&R, &S3, &R);
- return 3;
- }
- else {
- sub_big(&R, &S4, &R);
- return 4;
- }
- else if (big_comp(&R, &S7) < 0)
- if (big_comp(&R, &S6) < 0) {
- sub_big(&R, &S5, &R);
- return 5;
- }
- else {
- sub_big(&R, &S6, &R);
- return 6;
- }
- else if (big_comp(&R, &S9) < 0)
- if (big_comp(&R, &S8) < 0) {
- sub_big(&R, &S7, &R);
- return 7;
- }
- else {
- sub_big(&R, &S8, &R);
- return 8;
- }
- else {
- sub_big(&R, &S9, &R);
- return 9;
- }
-}
-
-#define OUTDIG(d) { *buf++ = (d) + '0'; *buf = 0; return k; }
-
-int
-_PyFloat_Digits(char *buf, double v, int *signum)
-{
- struct dblflt *x;
- int sign, e, f_n, m_n, i, d, tc1, tc2;
- Bigit f;
-
- /* decompose float into sign, mantissa & exponent */
- x = (struct dblflt *)&v;
- sign = x->s;
- e = x->e;
- f = (Bigit)(x->m1 << 16 | x->m2) << 32 | (U32)(x->m3 << 16 | x->m4);
- if (e != 0) {
- e = e - bias - bitstoright;
- f |= (Bigit)hidden_bit << 32;
- }
- else if (f != 0)
- /* denormalized */
- e = 1 - bias - bitstoright;
-
- *signum = sign;
- if (f == 0) {
- *buf++ = '0';
- *buf = 0;
- return 0;
- }
-
- ruf = !(f & 1); /* ruf = (even? f) */
-
- /* Compute the scaling factor estimate, k */
- if (e > MIN_E)
- k = estimate(e+52);
- else {
- int n;
- Bigit y;
-
- for (n = e+52, y = (Bigit)1 << 52; f < y; n--) y >>= 1;
- k = estimate(n);
- }
-
- if (e >= 0)
- if (f != B_P1)
- use_mp = 0, f_n = e+1, s_n = 1, m_n = e;
- else
- use_mp = 1, f_n = e+2, s_n = 2, m_n = e;
- else
- if ((e == MIN_E) || (f != B_P1))
- use_mp = 0, f_n = 1, s_n = 1-e, m_n = 0;
- else
- use_mp = 1, f_n = 2, s_n = 2-e, m_n = 0;
-
- /* Scale it! */
- if (k == 0) {
- short_shift_left(f, f_n, &R);
- one_shift_left(s_n, &S);
- one_shift_left(m_n, &MM);
- if (use_mp) one_shift_left(m_n+1, &MP);
- qr_shift = 1;
- }
- else if (k > 0) {
- s_n += k;
- if (m_n >= s_n)
- f_n -= s_n, m_n -= s_n, s_n = 0;
- else
- f_n -= m_n, s_n -= m_n, m_n = 0;
- short_shift_left(f, f_n, &R);
- big_shift_left(&five[k-1], s_n, &S);
- one_shift_left(m_n, &MM);
- if (use_mp) one_shift_left(m_n+1, &MP);
- qr_shift = 0;
- }
- else {
- Bignum *power = &five[-k-1];
-
- s_n += k;
- big_short_mul(power, f, &S);
- big_shift_left(&S, f_n, &R);
- one_shift_left(s_n, &S);
- big_shift_left(power, m_n, &MM);
- if (use_mp) big_shift_left(power, m_n+1, &MP);
- qr_shift = 1;
- }
-
- /* fixup */
- if (add_cmp() <= -ruf) {
- k--;
- mul10(&R);
- mul10(&MM);
- if (use_mp) mul10(&MP);
- }
-
- /*
- printf("\nk = %d\n", k);
- printf("R = "); print_big(&R);
- printf("\nS = "); print_big(&S);
- printf("\nM- = "); print_big(&MM);
- if (use_mp) printf("\nM+ = "), print_big(&MP);
- putchar('\n');
- fflush(0);
- */
-
- if (qr_shift) {
- sl = s_n / 64;
- slr = s_n % 64;
- }
- else {
- big_shift_left(&S, 1, &S2);
- add_big(&S2, &S, &S3);
- big_shift_left(&S2, 1, &S4);
- add_big(&S4, &S, &S5);
- add_big(&S4, &S2, &S6);
- add_big(&S4, &S3, &S7);
- big_shift_left(&S4, 1, &S8);
- add_big(&S8, &S, &S9);
- }
-
-again:
- if (qr_shift) { /* Take advantage of the fact that S = (ash 1 s_n) */
- if (R.l < sl)
- d = 0;
- else if (R.l == sl) {
- Bigit *p;
-
- p = &R.d[sl];
- d = *p >> slr;
- *p &= ((Bigit)1 << slr) - 1;
- for (i = sl; (i > 0) && (*p == 0); i--) p--;
- R.l = i;
- }
- else {
- Bigit *p;
-
- p = &R.d[sl+1];
- d = *p << (64 - slr) | *(p-1) >> slr;
- p--;
- *p &= ((Bigit)1 << slr) - 1;
- for (i = sl; (i > 0) && (*p == 0); i--) p--;
- R.l = i;
- }
- }
- else /* We need to do quotient-remainder */
- d = qr();
-
- tc1 = big_comp(&R, &MM) < ruf;
- tc2 = add_cmp() > -ruf;
- if (!tc1)
- if (!tc2) {
- mul10(&R);
- mul10(&MM);
- if (use_mp) mul10(&MP);
- *buf++ = d + '0';
- goto again;
- }
- else
- OUTDIG(d+1)
- else
- if (!tc2)
- OUTDIG(d)
- else {
- big_shift_left(&R, 1, &MM);
- if (big_comp(&MM, &S) == -1)
- OUTDIG(d)
- else
- OUTDIG(d+1)
- }
-}
-
-void
-_PyFloat_DigitsInit()
-{
- int n, i, l;
- Bignum *b;
- Bigit *xp, *zp, k;
-
- five[0].l = l = 0;
- five[0].d[0] = 5;
- for (n = MAX_FIVE-1, b = &five[0]; n > 0; n--) {
- xp = &b->d[0];
- b++;
- zp = &b->d[0];
- for (i = l, k = 0; i >= 0; i--)
- MUL(*xp++, 5, *zp++, k);
- if (k != 0)
- *zp = k, l++;
- b->l = l;
- }
-
- /*
- for (n = 1, b = &five[0]; n <= MAX_FIVE; n++) {
- big_shift_left(b++, n, &R);
- print_big(&R);
- putchar('\n');
- }
- fflush(0);
- */
-}
Modified: python/branches/py3k/Objects/floatobject.c
==============================================================================
--- python/branches/py3k/Objects/floatobject.c (original)
+++ python/branches/py3k/Objects/floatobject.c Sat Apr 19 02:31:39 2008
@@ -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 @@
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_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 @@
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 @@
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 @@
{
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 @@
}
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.
@@ -995,6 +891,57 @@
}
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)
{
double x = PyFloat_AsDouble(v);
@@ -1368,7 +1315,7 @@
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 @@
"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 @@
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);
Modified: python/branches/py3k/Objects/longobject.c
==============================================================================
--- python/branches/py3k/Objects/longobject.c (original)
+++ python/branches/py3k/Objects/longobject.c Sat Apr 19 02:31:39 2008
@@ -3611,9 +3611,21 @@
#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,
Modified: python/branches/py3k/PC/VC6/pythoncore.dsp
==============================================================================
--- python/branches/py3k/PC/VC6/pythoncore.dsp (original)
+++ python/branches/py3k/PC/VC6/pythoncore.dsp Sat Apr 19 02:31:39 2008
@@ -587,6 +587,10 @@
# 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
Modified: python/branches/py3k/PC/VS7.1/pythoncore.vcproj
==============================================================================
--- python/branches/py3k/PC/VS7.1/pythoncore.vcproj (original)
+++ python/branches/py3k/PC/VS7.1/pythoncore.vcproj Sat Apr 19 02:31:39 2008
@@ -698,6 +698,9 @@
RelativePath="..\..\Python\pyfpe.c">
</File>
<File
+ RelativePath="..\..\Python\pymath.c">
+ </File>
+ <File
RelativePath="..\..\Python\pystate.c">
</File>
<File
Modified: python/branches/py3k/PC/VS8.0/pythoncore.vcproj
==============================================================================
--- python/branches/py3k/PC/VS8.0/pythoncore.vcproj (original)
+++ python/branches/py3k/PC/VS8.0/pythoncore.vcproj Sat Apr 19 02:31:39 2008
@@ -1707,6 +1707,10 @@
>
</File>
<File
+ RelativePath="..\..\Python\pymath.c"
+ >
+ </File>
+ <File
RelativePath="..\..\Python\pystate.c"
>
</File>
Modified: python/branches/py3k/PC/pyconfig.h
==============================================================================
--- python/branches/py3k/PC/pyconfig.h (original)
+++ python/branches/py3k/PC/pyconfig.h Sat Apr 19 02:31:39 2008
@@ -207,12 +207,13 @@
#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 @@
/* 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
Modified: python/branches/py3k/PCbuild/pythoncore.vcproj
==============================================================================
--- python/branches/py3k/PCbuild/pythoncore.vcproj (original)
+++ python/branches/py3k/PCbuild/pythoncore.vcproj Sat Apr 19 02:31:39 2008
@@ -871,6 +871,10 @@
>
</File>
<File
+ RelativePath="..\Include\pymath.h"
+ >
+ </File>
+ <File
RelativePath="..\Include\pymem.h"
>
</File>
@@ -1707,6 +1711,10 @@
>
</File>
<File
+ RelativePath="..\Python\pymath.c"
+ >
+ </File>
+ <File
RelativePath="..\Python\pystate.c"
>
</File>
Deleted: python/branches/py3k/Python/hypot.c
==============================================================================
--- python/branches/py3k/Python/hypot.c Sat Apr 19 02:31:39 2008
+++ (empty file)
@@ -1,25 +0,0 @@
-/* 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 */
-
Modified: python/branches/py3k/configure
==============================================================================
--- python/branches/py3k/configure (original)
+++ python/branches/py3k/configure Sat Apr 19 02:31:39 2008
@@ -1,5 +1,5 @@
#! /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.
#
More information about the Python-checkins
mailing list