From 9ab44b509a935011beb8e9108a2271ee728e8ad4 Mon Sep 17 00:00:00 2001 From: Mark Dickinson Date: Wed, 30 Dec 2009 16:22:49 +0000 Subject: Merged revisions 77139-77140 via svnmerge from svn+ssh://pythondev@svn.python.org/python/trunk ........ r77139 | mark.dickinson | 2009-12-30 12:12:23 +0000 (Wed, 30 Dec 2009) | 3 lines Issue #7534: Fix handling of nans, infinities, and negative zero in ** operator, on IEEE 754 platforms. Thanks Marcos Donolo for original patch. ........ r77140 | mark.dickinson | 2009-12-30 12:22:49 +0000 (Wed, 30 Dec 2009) | 1 line Add Marcos Donolo for work on issue 7534 patch. ........ --- Lib/test/ieee754.txt | 6 +- Lib/test/test_float.py | 211 +++++++++++++++++++++++++++++++++++++++++++++++++ Misc/ACKS | 1 + Misc/NEWS | 4 + Objects/floatobject.c | 97 +++++++++++++++++------ 5 files changed, 294 insertions(+), 25 deletions(-) diff --git a/Lib/test/ieee754.txt b/Lib/test/ieee754.txt index 5a41c8f..89bb0c5 100644 --- a/Lib/test/ieee754.txt +++ b/Lib/test/ieee754.txt @@ -72,7 +72,7 @@ False >>> NAN >= 0 False -All operations involving a NaN return a NaN except for the power of *0* and *1*. +All operations involving a NaN return a NaN except for nan**0 and 1**nan. >>> 1 + NAN nan >>> 1 * NAN @@ -81,8 +81,10 @@ nan nan >>> 1 ** NAN 1.0 +>>> NAN ** 0 +1.0 >>> 0 ** NAN -0.0 +nan >>> (1.0 + FI.epsilon) * NAN nan diff --git a/Lib/test/test_float.py b/Lib/test/test_float.py index 0f984a5..e2d1ea0 100644 --- a/Lib/test/test_float.py +++ b/Lib/test/test_float.py @@ -11,6 +11,11 @@ import random, fractions INF = float("inf") NAN = float("nan") +# decorator for skipping tests on non-IEEE 754 platforms +requires_IEEE_754 = unittest.skipUnless( + float.__getformat__("double").startswith("IEEE"), + "test requires IEEE 754 doubles") + #locate file with float format test values test_dir = os.path.dirname(__file__) or os.curdir format_testfile = os.path.join(test_dir, 'formatfloat_testcases.txt') @@ -161,6 +166,212 @@ class GeneralFloatCases(unittest.TestCase): self.assertTrue(s == s, "{%r} not equal to itself" % f) self.assertTrue(d == d, "{%r : None} not equal to itself" % f) + def assertEqualAndEqualSign(self, a, b): + # fail unless a == b and a and b have the same sign bit; + # the only difference from assertEqual is that this test + # distingishes -0.0 and 0.0. + self.assertEqual((a, copysign(1.0, a)), (b, copysign(1.0, b))) + + @requires_IEEE_754 + def test_float_pow(self): + # test builtin pow and ** operator for IEEE 754 special cases. + # Special cases taken from section F.9.4.4 of the C99 specification + + for pow_op in pow, operator.pow: + # x**NAN is NAN for any x except 1 + self.assertTrue(isnan(pow_op(-INF, NAN))) + self.assertTrue(isnan(pow_op(-2.0, NAN))) + self.assertTrue(isnan(pow_op(-1.0, NAN))) + self.assertTrue(isnan(pow_op(-0.5, NAN))) + self.assertTrue(isnan(pow_op(-0.0, NAN))) + self.assertTrue(isnan(pow_op(0.0, NAN))) + self.assertTrue(isnan(pow_op(0.5, NAN))) + self.assertTrue(isnan(pow_op(2.0, NAN))) + self.assertTrue(isnan(pow_op(INF, NAN))) + self.assertTrue(isnan(pow_op(NAN, NAN))) + + # NAN**y is NAN for any y except +-0 + self.assertTrue(isnan(pow_op(NAN, -INF))) + self.assertTrue(isnan(pow_op(NAN, -2.0))) + self.assertTrue(isnan(pow_op(NAN, -1.0))) + self.assertTrue(isnan(pow_op(NAN, -0.5))) + self.assertTrue(isnan(pow_op(NAN, 0.5))) + self.assertTrue(isnan(pow_op(NAN, 1.0))) + self.assertTrue(isnan(pow_op(NAN, 2.0))) + self.assertTrue(isnan(pow_op(NAN, INF))) + + # (+-0)**y raises ZeroDivisionError for y a negative odd integer + self.assertRaises(ZeroDivisionError, pow_op, -0.0, -1.0) + self.assertRaises(ZeroDivisionError, pow_op, 0.0, -1.0) + + # (+-0)**y raises ZeroDivisionError for y finite and negative + # but not an odd integer + self.assertRaises(ZeroDivisionError, pow_op, -0.0, -2.0) + self.assertRaises(ZeroDivisionError, pow_op, -0.0, -0.5) + self.assertRaises(ZeroDivisionError, pow_op, 0.0, -2.0) + self.assertRaises(ZeroDivisionError, pow_op, 0.0, -0.5) + + # (+-0)**y is +-0 for y a positive odd integer + self.assertEqualAndEqualSign(pow_op(-0.0, 1.0), -0.0) + self.assertEqualAndEqualSign(pow_op(0.0, 1.0), 0.0) + + # (+-0)**y is 0 for y finite and positive but not an odd integer + self.assertEqualAndEqualSign(pow_op(-0.0, 0.5), 0.0) + self.assertEqualAndEqualSign(pow_op(-0.0, 2.0), 0.0) + self.assertEqualAndEqualSign(pow_op(0.0, 0.5), 0.0) + self.assertEqualAndEqualSign(pow_op(0.0, 2.0), 0.0) + + # (-1)**+-inf is 1 + self.assertEqualAndEqualSign(pow_op(-1.0, -INF), 1.0) + self.assertEqualAndEqualSign(pow_op(-1.0, INF), 1.0) + + # 1**y is 1 for any y, even if y is an infinity or nan + self.assertEqualAndEqualSign(pow_op(1.0, -INF), 1.0) + self.assertEqualAndEqualSign(pow_op(1.0, -2.0), 1.0) + self.assertEqualAndEqualSign(pow_op(1.0, -1.0), 1.0) + self.assertEqualAndEqualSign(pow_op(1.0, -0.5), 1.0) + self.assertEqualAndEqualSign(pow_op(1.0, -0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(1.0, 0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(1.0, 0.5), 1.0) + self.assertEqualAndEqualSign(pow_op(1.0, 1.0), 1.0) + self.assertEqualAndEqualSign(pow_op(1.0, 2.0), 1.0) + self.assertEqualAndEqualSign(pow_op(1.0, INF), 1.0) + self.assertEqualAndEqualSign(pow_op(1.0, NAN), 1.0) + + # x**+-0 is 1 for any x, even if x is a zero, infinity, or nan + self.assertEqualAndEqualSign(pow_op(-INF, 0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(-2.0, 0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(-1.0, 0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(-0.5, 0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(-0.0, 0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(0.0, 0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(0.5, 0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(1.0, 0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(2.0, 0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(INF, 0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(NAN, 0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(-INF, -0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(-2.0, -0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(-1.0, -0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(-0.5, -0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(-0.0, -0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(0.0, -0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(0.5, -0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(1.0, -0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(2.0, -0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(INF, -0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(NAN, -0.0), 1.0) + + # x**y defers to complex pow for finite negative x and + # non-integral y. + self.assertEqual(type(pow_op(-2.0, -0.5)), complex) + self.assertEqual(type(pow_op(-2.0, 0.5)), complex) + self.assertEqual(type(pow_op(-1.0, -0.5)), complex) + self.assertEqual(type(pow_op(-1.0, 0.5)), complex) + self.assertEqual(type(pow_op(-0.5, -0.5)), complex) + self.assertEqual(type(pow_op(-0.5, 0.5)), complex) + + # x**-INF is INF for abs(x) < 1 + self.assertEqualAndEqualSign(pow_op(-0.5, -INF), INF) + self.assertEqualAndEqualSign(pow_op(-0.0, -INF), INF) + self.assertEqualAndEqualSign(pow_op(0.0, -INF), INF) + self.assertEqualAndEqualSign(pow_op(0.5, -INF), INF) + + # x**-INF is 0 for abs(x) > 1 + self.assertEqualAndEqualSign(pow_op(-INF, -INF), 0.0) + self.assertEqualAndEqualSign(pow_op(-2.0, -INF), 0.0) + self.assertEqualAndEqualSign(pow_op(2.0, -INF), 0.0) + self.assertEqualAndEqualSign(pow_op(INF, -INF), 0.0) + + # x**INF is 0 for abs(x) < 1 + self.assertEqualAndEqualSign(pow_op(-0.5, INF), 0.0) + self.assertEqualAndEqualSign(pow_op(-0.0, INF), 0.0) + self.assertEqualAndEqualSign(pow_op(0.0, INF), 0.0) + self.assertEqualAndEqualSign(pow_op(0.5, INF), 0.0) + + # x**INF is INF for abs(x) > 1 + self.assertEqualAndEqualSign(pow_op(-INF, INF), INF) + self.assertEqualAndEqualSign(pow_op(-2.0, INF), INF) + self.assertEqualAndEqualSign(pow_op(2.0, INF), INF) + self.assertEqualAndEqualSign(pow_op(INF, INF), INF) + + # (-INF)**y is -0.0 for y a negative odd integer + self.assertEqualAndEqualSign(pow_op(-INF, -1.0), -0.0) + + # (-INF)**y is 0.0 for y negative but not an odd integer + self.assertEqualAndEqualSign(pow_op(-INF, -0.5), 0.0) + self.assertEqualAndEqualSign(pow_op(-INF, -2.0), 0.0) + + # (-INF)**y is -INF for y a positive odd integer + self.assertEqualAndEqualSign(pow_op(-INF, 1.0), -INF) + + # (-INF)**y is INF for y positive but not an odd integer + self.assertEqualAndEqualSign(pow_op(-INF, 0.5), INF) + self.assertEqualAndEqualSign(pow_op(-INF, 2.0), INF) + + # INF**y is INF for y positive + self.assertEqualAndEqualSign(pow_op(INF, 0.5), INF) + self.assertEqualAndEqualSign(pow_op(INF, 1.0), INF) + self.assertEqualAndEqualSign(pow_op(INF, 2.0), INF) + + # INF**y is 0.0 for y negative + self.assertEqualAndEqualSign(pow_op(INF, -2.0), 0.0) + self.assertEqualAndEqualSign(pow_op(INF, -1.0), 0.0) + self.assertEqualAndEqualSign(pow_op(INF, -0.5), 0.0) + + # basic checks not covered by the special cases above + self.assertEqualAndEqualSign(pow_op(-2.0, -2.0), 0.25) + self.assertEqualAndEqualSign(pow_op(-2.0, -1.0), -0.5) + self.assertEqualAndEqualSign(pow_op(-2.0, -0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(-2.0, 0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(-2.0, 1.0), -2.0) + self.assertEqualAndEqualSign(pow_op(-2.0, 2.0), 4.0) + self.assertEqualAndEqualSign(pow_op(-1.0, -2.0), 1.0) + self.assertEqualAndEqualSign(pow_op(-1.0, -1.0), -1.0) + self.assertEqualAndEqualSign(pow_op(-1.0, -0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(-1.0, 0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(-1.0, 1.0), -1.0) + self.assertEqualAndEqualSign(pow_op(-1.0, 2.0), 1.0) + self.assertEqualAndEqualSign(pow_op(2.0, -2.0), 0.25) + self.assertEqualAndEqualSign(pow_op(2.0, -1.0), 0.5) + self.assertEqualAndEqualSign(pow_op(2.0, -0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(2.0, 0.0), 1.0) + self.assertEqualAndEqualSign(pow_op(2.0, 1.0), 2.0) + self.assertEqualAndEqualSign(pow_op(2.0, 2.0), 4.0) + + # 1 ** large and -1 ** large; some libms apparently + # have problems with these + self.assertEqualAndEqualSign(pow_op(1.0, -1e100), 1.0) + self.assertEqualAndEqualSign(pow_op(1.0, 1e100), 1.0) + self.assertEqualAndEqualSign(pow_op(-1.0, -1e100), 1.0) + self.assertEqualAndEqualSign(pow_op(-1.0, 1e100), 1.0) + + # check sign for results that underflow to 0 + self.assertEqualAndEqualSign(pow_op(-2.0, -2000.0), 0.0) + self.assertEqual(type(pow_op(-2.0, -2000.5)), complex) + self.assertEqualAndEqualSign(pow_op(-2.0, -2001.0), -0.0) + self.assertEqualAndEqualSign(pow_op(2.0, -2000.0), 0.0) + self.assertEqualAndEqualSign(pow_op(2.0, -2000.5), 0.0) + self.assertEqualAndEqualSign(pow_op(2.0, -2001.0), 0.0) + self.assertEqualAndEqualSign(pow_op(-0.5, 2000.0), 0.0) + self.assertEqual(type(pow_op(-0.5, 2000.5)), complex) + self.assertEqualAndEqualSign(pow_op(-0.5, 2001.0), -0.0) + self.assertEqualAndEqualSign(pow_op(0.5, 2000.0), 0.0) + self.assertEqualAndEqualSign(pow_op(0.5, 2000.5), 0.0) + self.assertEqualAndEqualSign(pow_op(0.5, 2001.0), 0.0) + + # check we don't raise an exception for subnormal results, + # and validate signs. Tests currently disabled, since + # they fail on systems where a subnormal result from pow + # is flushed to zero (e.g. Debian/ia64.) + #self.assertTrue(0.0 < pow_op(0.5, 1048) < 1e-315) + #self.assertTrue(0.0 < pow_op(-0.5, 1048) < 1e-315) + #self.assertTrue(0.0 < pow_op(0.5, 1047) < 1e-315) + #self.assertTrue(0.0 > pow_op(-0.5, 1047) > -1e-315) + #self.assertTrue(0.0 < pow_op(2.0, -1048) < 1e-315) + #self.assertTrue(0.0 < pow_op(-2.0, -1048) < 1e-315) + #self.assertTrue(0.0 < pow_op(2.0, -1047) < 1e-315) + #self.assertTrue(0.0 > pow_op(-2.0, -1047) > -1e-315) class FormatFunctionsTestCase(unittest.TestCase): diff --git a/Misc/ACKS b/Misc/ACKS index 0098bce..d4099ea 100644 --- a/Misc/ACKS +++ b/Misc/ACKS @@ -186,6 +186,7 @@ Yves Dionne Daniel Dittmar Jaromir Dolecek Ismail Donmez +Marcos Donolo Dima Dorfman Cesar Douady Dean Draayer diff --git a/Misc/NEWS b/Misc/NEWS index 6062f0d..923cd30 100644 --- a/Misc/NEWS +++ b/Misc/NEWS @@ -12,6 +12,10 @@ What's New in Python 3.2 Alpha 1? Core and Builtins ----------------- +- Issue #7534: Fix handling of IEEE specials (infinities, nans, + negative zero) in ** operator. The behaviour now conforms to that + described in C99 Annex F. + - Issue #1811: improve accuracy and cross-platform consistency for true division of integers: the result of a/b is now correctly rounded for ints a and b (at least on IEEE 754 platforms), and in diff --git a/Objects/floatobject.c b/Objects/floatobject.c index b281f81..33196ff 100644 --- a/Objects/floatobject.c +++ b/Objects/floatobject.c @@ -671,10 +671,15 @@ float_floor_div(PyObject *v, PyObject *w) return r; } +/* determine whether x is an odd integer or not; assumes that + x is not an infinity or nan. */ +#define DOUBLE_IS_ODD_INTEGER(x) (fmod(fabs(x), 2.0) == 1.0) + static PyObject * float_pow(PyObject *v, PyObject *w, PyObject *z) { double iv, iw, ix; + int negate_result = 0; if ((PyObject *)z != Py_None) { PyErr_SetString(PyExc_TypeError, "pow() 3rd argument not " @@ -686,20 +691,56 @@ float_pow(PyObject *v, PyObject *w, PyObject *z) CONVERT_TO_DOUBLE(w, iw); /* Sort out special cases here instead of relying on pow() */ - if (iw == 0) { /* v**0 is 1, even 0**0 */ + if (iw == 0) { /* v**0 is 1, even 0**0 */ return PyFloat_FromDouble(1.0); } - if (iv == 0.0) { /* 0**w is error if w<0, else 1 */ + if (Py_IS_NAN(iv)) { /* nan**w = nan, unless w == 0 */ + return PyFloat_FromDouble(iv); + } + if (Py_IS_NAN(iw)) { /* v**nan = nan, unless v == 1; 1**nan = 1 */ + return PyFloat_FromDouble(iv == 1.0 ? 1.0 : iw); + } + if (Py_IS_INFINITY(iw)) { + /* v**inf is: 0.0 if abs(v) < 1; 1.0 if abs(v) == 1; inf if + * abs(v) > 1 (including case where v infinite) + * + * v**-inf is: inf if abs(v) < 1; 1.0 if abs(v) == 1; 0.0 if + * abs(v) > 1 (including case where v infinite) + */ + iv = fabs(iv); + if (iv == 1.0) + return PyFloat_FromDouble(1.0); + else if ((iw > 0.0) == (iv > 1.0)) + return PyFloat_FromDouble(fabs(iw)); /* return inf */ + else + return PyFloat_FromDouble(0.0); + } + if (Py_IS_INFINITY(iv)) { + /* (+-inf)**w is: inf for w positive, 0 for w negative; in + * both cases, we need to add the appropriate sign if w is + * an odd integer. + */ + int iw_is_odd = DOUBLE_IS_ODD_INTEGER(iw); + if (iw > 0.0) + return PyFloat_FromDouble(iw_is_odd ? iv : fabs(iv)); + else + return PyFloat_FromDouble(iw_is_odd ? + copysign(0.0, iv) : 0.0); + } + if (iv == 0.0) { /* 0**w is: 0 for w positive, 1 for w zero + (already dealt with above), and an error + if w is negative. */ + int iw_is_odd = DOUBLE_IS_ODD_INTEGER(iw); if (iw < 0.0) { PyErr_SetString(PyExc_ZeroDivisionError, - "0.0 cannot be raised to a negative power"); + "0.0 cannot be raised to a " + "negative power"); return NULL; } - return PyFloat_FromDouble(0.0); - } - if (iv == 1.0) { /* 1**w is 1, even 1**inf and 1**nan */ - return PyFloat_FromDouble(1.0); + /* use correct sign if iw is odd */ + return PyFloat_FromDouble(iw_is_odd ? iv : 0.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. @@ -710,33 +751,41 @@ float_pow(PyObject *v, PyObject *w, PyObject *z) */ return PyComplex_Type.tp_as_number->nb_power(v, w, z); } - /* iw is an exact integer, albeit perhaps a very large one. + /* iw is an exact integer, albeit perhaps a very large + * one. Replace iv by its absolute value and remember + * to negate the pow result if iw is odd. + */ + iv = -iv; + negate_result = DOUBLE_IS_ODD_INTEGER(iw); + } + + if (iv == 1.0) { /* 1**w is 1, even 1**inf and 1**nan */ + /* (-1) ** large_integer also ends up here. Here's an + * extract from the comments for the previous + * implementation explaining why this special case is + * necessary: + * * -1 raised to an exact integer should never be exceptional. * Alas, some libms (chiefly glibc as of early 2003) return * NaN and set EDOM on pow(-1, large_int) if the int doesn't * happen to be representable in a *C* integer. That's a - * bug; we let that slide in math.pow() (which currently - * reflects all platform accidents), but not for Python's **. - */ - if (iv == -1.0 && Py_IS_FINITE(iw)) { - /* Return 1 if iw is even, -1 if iw is odd; there's - * no guarantee that any C integral type is big - * enough to hold iw, so we have to check this - * indirectly. - */ - ix = floor(iw * 0.5) * 2.0; - return PyFloat_FromDouble(ix == iw ? 1.0 : -1.0); - } - /* Else iv != -1.0, and overflow or underflow are possible. - * Unless we're to write pow() ourselves, we have to trust - * the platform to do this correctly. + * bug. */ + return PyFloat_FromDouble(negate_result ? -1.0 : 1.0); } + + /* Now iv and iw are finite, iw is nonzero, and iv is + * positive and not equal to 1.0. We finally allow + * the platform pow to step in and do the rest. + */ errno = 0; PyFPE_START_PROTECT("pow", return NULL) ix = pow(iv, iw); PyFPE_END_PROTECT(ix) Py_ADJUST_ERANGE1(ix); + if (negate_result) + ix = -ix; + if (errno != 0) { /* We don't expect any errno value other than ERANGE, but * the range of libm bugs appears unbounded. @@ -748,6 +797,8 @@ float_pow(PyObject *v, PyObject *w, PyObject *z) return PyFloat_FromDouble(ix); } +#undef DOUBLE_IS_ODD_INTEGER + static PyObject * float_neg(PyFloatObject *v) { -- cgit v0.12