]> granicus.if.org Git - python/commitdiff
Merged revisions 77139-77140 via svnmerge from
authorMark Dickinson <dickinsm@gmail.com>
Wed, 30 Dec 2009 16:22:49 +0000 (16:22 +0000)
committerMark Dickinson <dickinsm@gmail.com>
Wed, 30 Dec 2009 16:22:49 +0000 (16:22 +0000)
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
Lib/test/test_float.py
Misc/ACKS
Misc/NEWS
Objects/floatobject.c

index 5a41c8fe5077a7460752c809662e098dfeada50b..89bb0c50a3d24a1ce4445575c585f93f13f67b0c 100644 (file)
@@ -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
 
index 0f984a5ca20c999009a2da6dda84ed00b923fda7..e2d1ea0a7cd630323b0199fa9bec898c08c2c603 100644 (file)
@@ -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):
index 0098bce42361c89dc132004831fc2da45f4c5014..d4099eabde5b1b6e46ab556170c813bcbdcd1697 100644 (file)
--- 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
index 6062f0d560df5d1779c65332d6acc60728cb68c8..923cd30a8194e93514cbf19321df6a2c31fc21fa 100644 (file)
--- 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
index b281f815dc750a755014b1b8c39ef3e5a13d3217..33196ff476d86ffb7c5886f9ea6f0ae3693c9ca0 100644 (file)
@@ -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)
 {