From 55f679fbffd43d2410dc3e261d332cb6c04bfbfb Mon Sep 17 00:00:00 2001 From: Mark Dickinson Date: Sat, 1 Apr 2017 15:36:14 +0100 Subject: [PATCH 1/6] Implement math.remainder. --- Doc/library/math.rst | 21 +++++++ Doc/whatsnew/3.7.rst | 6 ++ Lib/test/test_math.py | 135 ++++++++++++++++++++++++++++++++++++++++++ Misc/NEWS | 3 + Modules/mathmodule.c | 55 +++++++++++++++++ 5 files changed, 220 insertions(+) diff --git a/Doc/library/math.rst b/Doc/library/math.rst index da2b8cc5862791..cfc0f725e9bf1e 100644 --- a/Doc/library/math.rst +++ b/Doc/library/math.rst @@ -175,6 +175,27 @@ Number-theoretic and representation functions of *x* and are floats. +.. function:: remainder(x, y) + + Return the IEEE 754-style remainder of ``x`` with respect to ``y``. For + finite ``x`` and finite nonzero ``y``, this is the difference ``x - n*y``, + where ``n`` is the closest integer to the exact value of the quotient ``x / + y``. If ``x / y`` is exactly halfway between two consecutive integers, the + nearest *even* integer is used for ``n``. The remainder ``r = remainder(x, + y)`` thus always satisfies ``abs(r) <= 0.5 * abs(y)``. + + Special cases follow IEEE 754: in particular, ``remainder(x, math.inf)`` is + ``x`` for any finite ``x``, and ``remainder(x, 0)`` and + ``remainder(math.inf, x)`` raise :exc:`ValueError` for any non-NaN ``x``. + If the result of the remainder operation is zero, that zero will have + the same sign as ``x``. + + On platforms using IEEE 754 binary floating-point, the result of this + operation is always exactly representable: no rounding error is introduced. + + .. versionadded:: 3.7 + + .. function:: trunc(x) Return the :class:`~numbers.Real` value *x* truncated to an diff --git a/Doc/whatsnew/3.7.rst b/Doc/whatsnew/3.7.rst index 19e04bb19efc5e..a55215d17b5817 100644 --- a/Doc/whatsnew/3.7.rst +++ b/Doc/whatsnew/3.7.rst @@ -111,6 +111,12 @@ Serhiy Storchaka in :issue:`28682`.) Added support for :ref:`file descriptors ` in :func:`~os.scandir` on Unix. (Contributed by Serhiy Storchaka in :issue:`25996`.) +math +---- + +New ``remainder`` function, implementing the IEEE 754-style remainder +operation. (Constructed by Mark Dickinson in :issue:`29962`.) + unittest.mock ------------- diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py index eaa41bca3f686e..70cb57465f50f1 100644 --- a/Lib/test/test_math.py +++ b/Lib/test/test_math.py @@ -1000,6 +1000,135 @@ def testRadians(self): self.ftest('radians(-45)', math.radians(-45), -math.pi/4) self.ftest('radians(0)', math.radians(0), 0) + @requires_IEEE_754 + def testRemainder(self): + from fractions import Fraction + + def validate_spec(x, y, r): + """ + Check that r matches remainder(x, y) according to the IEEE 754 + specification. Assumes that x, y and r are finite and y is nonzero. + """ + fx, fy, fr = Fraction(x), Fraction(y), Fraction(r) + # r should not exceed y/2 in absolute value + self.assertLessEqual(abs(fr), abs(fy/2)) + # x - r should be an exact integer multiple of y + n = (fx - fr) / fy + self.assertEqual(n, int(n)) + if abs(fr) == abs(fy/2): + # If |r| == |y/2|, n should be even. + self.assertEqual(n/2, int(n/2)) + + # triples (x, y, remainder(x, y)) in hexadecimal form. + testcases = [ + # Remainders modulo 1, showing the ties-to-even behaviour. + '-4.0 1 -0.0', + '-3.8 1 0.8', + '-3.0 1 -0.0', + '-2.8 1 -0.8', + '-2.0 1 -0.0', + '-1.8 1 0.8', + '-1.0 1 -0.0', + '-0.8 1 -0.8', + '-0.0 1 -0.0', + ' 0.0 1 0.0', + ' 0.8 1 0.8', + ' 1.0 1 0.0', + ' 1.8 1 -0.8', + ' 2.0 1 0.0', + ' 2.8 1 0.8', + ' 3.0 1 0.0', + ' 3.8 1 -0.8', + ' 4.0 1 0.0', + + # Reductions modulo 2*pi + '0x0.0p+0 0x1.921fb54442d18p+2 0x0.0p+0', + '0x1.921fb54442d18p+0 0x1.921fb54442d18p+2 0x1.921fb54442d18p+0', + '0x1.921fb54442d17p+1 0x1.921fb54442d18p+2 0x1.921fb54442d17p+1', + '0x1.921fb54442d18p+1 0x1.921fb54442d18p+2 0x1.921fb54442d18p+1', + '0x1.921fb54442d19p+1 0x1.921fb54442d18p+2 -0x1.921fb54442d17p+1', + '0x1.921fb54442d17p+2 0x1.921fb54442d18p+2 -0x0.0000000000001p+2', + '0x1.921fb54442d18p+2 0x1.921fb54442d18p+2 0x0p0', + '0x1.921fb54442d19p+2 0x1.921fb54442d18p+2 0x0.0000000000001p+2', + '0x1.2d97c7f3321d1p+3 0x1.921fb54442d18p+2 0x1.921fb54442d14p+1', + '0x1.2d97c7f3321d2p+3 0x1.921fb54442d18p+2 -0x1.921fb54442d18p+1', + '0x1.2d97c7f3321d3p+3 0x1.921fb54442d18p+2 -0x1.921fb54442d14p+1', + '0x1.921fb54442d17p+3 0x1.921fb54442d18p+2 -0x0.0000000000001p+3', + '0x1.921fb54442d18p+3 0x1.921fb54442d18p+2 0x0p0', + '0x1.921fb54442d19p+3 0x1.921fb54442d18p+2 0x0.0000000000001p+3', + '0x1.f6a7a2955385dp+3 0x1.921fb54442d18p+2 0x1.921fb54442d14p+1', + '0x1.f6a7a2955385ep+3 0x1.921fb54442d18p+2 0x1.921fb54442d18p+1', + '0x1.f6a7a2955385fp+3 0x1.921fb54442d18p+2 -0x1.921fb54442d14p+1', + '0x1.1475cc9eedf00p+5 0x1.921fb54442d18p+2 0x1.921fb54442d10p+1', + '0x1.1475cc9eedf01p+5 0x1.921fb54442d18p+2 -0x1.921fb54442d10p+1', + + # Symmetry with respect to signs. + ' 1 0.c 0.4', + '-1 0.c -0.4', + ' 1 -0.c 0.4', + '-1 -0.c -0.4', + ' 1.4 0.c -0.4', + '-1.4 0.c 0.4', + ' 1.4 -0.c -0.4', + '-1.4 -0.c 0.4', + + # Huge modulus, to check that the underlying algorithm doesn't + # rely on 2.0 * modulus being representable. + '0x1.dp+1023 0x1.4p+1023 0x0.9p+1023', + '0x1.ep+1023 0x1.4p+1023 -0x0.ap+1023', + '0x1.fp+1023 0x1.4p+1023 -0x0.9p+1023', + ] + + for case in testcases: + with self.subTest(case=case): + x_hex, y_hex, expected_hex = case.split() + x = float.fromhex(x_hex) + y = float.fromhex(y_hex) + expected = float.fromhex(expected_hex) + validate_spec(x, y, expected) + actual = math.remainder(x, y) + # Cheap way of checking that the floats are + # as identical as we need them to be. + self.assertEqual(actual.hex(), expected.hex()) + + # Test tiny subnormal modulus: there's potential for + # getting the implementation wrong here (for example, + # by assuming that modulus/2 is exactly representable). + tiny = float.fromhex('1p-1074') # min +ve subnormal + for n in range(-25, 25): + if n == 0: + continue + y = n * tiny + for m in range(100): + x = m * tiny + actual = math.remainder(x, y) + validate_spec(x, y, actual) + actual = math.remainder(-x, y) + validate_spec(-x, y, actual) + + # Special values. + # NaNs should propagate as usual. + for value in [NAN, 0.0, -0.0, 2.0, -2.3, NINF, INF]: + self.assertIsNaN(math.remainder(NAN, value)) + self.assertIsNaN(math.remainder(value, NAN)) + + # remainder(x, inf) is x, for non-nan non-infinite x. + for value in [-2.3, -0.0, 0.0, 2.3]: + self.assertEqual(math.remainder(value, INF), value) + self.assertEqual(math.remainder(value, NINF), value) + + # remainder(x, 0) and remainder(infinity, x) for non-NaN x are invalid + # operations according to IEEE 754-2008 7.2(f), and should raise. + for value in [NINF, -2.3, -0.0, 0.0, 2.3, INF]: + with self.assertRaises(ValueError): + math.remainder(INF, value) + with self.assertRaises(ValueError): + math.remainder(NINF, value) + with self.assertRaises(ValueError): + math.remainder(value, 0.0) + with self.assertRaises(ValueError): + math.remainder(value, -0.0) + def testSin(self): self.assertRaises(TypeError, math.sin) self.ftest('sin(0)', math.sin(0), 0) @@ -1286,6 +1415,12 @@ def test_mtestfile(self): self.fail('Failures in test_mtestfile:\n ' + '\n '.join(failures)) + # Custom assertions. + + def assertIsNaN(self, value): + if not math.isnan(value): + self.fail("Expected a NaN, got {!r}.".format(value)) + class IsCloseTests(unittest.TestCase): isclose = math.isclose # sublcasses should override this diff --git a/Misc/NEWS b/Misc/NEWS index a9acaf8e62f695..a2a58429231742 100644 --- a/Misc/NEWS +++ b/Misc/NEWS @@ -303,6 +303,9 @@ Extension Modules Library ------- +- bpo-29962: Add math.remainder operation, implementing remainder + as specified in IEEE 754. + - bpo-29931: Fixed comparison check for ipaddress.ip_interface objects. Patch by Sanjay Sundaresan. diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c index 243560a4d8bc43..43bb0d1fa30b63 100644 --- a/Modules/mathmodule.c +++ b/Modules/mathmodule.c @@ -600,6 +600,54 @@ m_atan2(double y, double x) return atan2(y, x); } + +/* IEEE 754-style remainder operation: x - n*y where n*y is the nearest + multiple of y to x, taking n even in the case of a tie. Assuming an IEEE 754 + binary floating-point format, the result is always exact. */ + +static double +m_remainder(double x, double y) +{ + /* Deal with most common case first. */ + if (Py_IS_FINITE(x) && Py_IS_FINITE(y)) { + double absx, absy, c, m, r; + + if (y == 0.0) { + return Py_NAN; + } + + absx = fabs(x); + absy = fabs(y); + m = fmod(absx, absy); + c = absy - m; + if (m < c) { + r = m; + } + else if (m > c) { + r = -c; + } + else { + assert(m == c); + r = m - 2.0 * fmod(0.5 * (absx - m), absy); + } + return copysign(1.0, x) * r; + } + + /* Special values. */ + if (Py_IS_NAN(x)) { + return x; + } + if (Py_IS_NAN(y)) { + return y; + } + if (Py_IS_INFINITY(x)) { + return Py_NAN; + } + assert(Py_IS_INFINITY(y)); + return x; +} + + /* Various platforms (Solaris, OpenBSD) do nonstandard things for log(0), log(-ve), log(NaN). Here are wrappers for log and log10 that deal with @@ -1072,6 +1120,12 @@ FUNC1(log1p, m_log1p, 0, "log1p($module, x, /)\n--\n\n" "Return the natural logarithm of 1+x (base e).\n\n" "The result is computed in a way which is accurate for x near zero.") +FUNC2(remainder, m_remainder, + "remainder($module, x, y, /)\n--\n\n" + "Difference between x and the closest integer multiple of y.\n\n" + "Return x - n*y where n*y is the closest integer multiple of y.\n" + "In the case where x is exactly halfway between two multiples of\n" + "y, the nearest even value of n is used. The result is always exact.") FUNC1(sin, sin, 0, "sin($module, x, /)\n--\n\n" "Return the sine of x (measured in radians).") @@ -2258,6 +2312,7 @@ static PyMethodDef math_methods[] = { MATH_MODF_METHODDEF MATH_POW_METHODDEF MATH_RADIANS_METHODDEF + {"remainder", math_remainder, METH_VARARGS, math_remainder_doc}, {"sin", math_sin, METH_O, math_sin_doc}, {"sinh", math_sinh, METH_O, math_sinh_doc}, {"sqrt", math_sqrt, METH_O, math_sqrt_doc}, From 9d50ace5bae6380fd7b39f2a055b40077b2a3f98 Mon Sep 17 00:00:00 2001 From: Mark Dickinson Date: Sat, 1 Apr 2017 16:12:53 +0100 Subject: [PATCH 2/6] Fix markup for arguments; use double spaces after period. --- Doc/library/math.rst | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/Doc/library/math.rst b/Doc/library/math.rst index cfc0f725e9bf1e..e65cacbf7695af 100644 --- a/Doc/library/math.rst +++ b/Doc/library/math.rst @@ -177,18 +177,18 @@ Number-theoretic and representation functions .. function:: remainder(x, y) - Return the IEEE 754-style remainder of ``x`` with respect to ``y``. For - finite ``x`` and finite nonzero ``y``, this is the difference ``x - n*y``, + Return the IEEE 754-style remainder of *x* with respect to *y*. For + finite *x* and finite nonzero *y*, this is the difference ``x - n*y``, where ``n`` is the closest integer to the exact value of the quotient ``x / - y``. If ``x / y`` is exactly halfway between two consecutive integers, the - nearest *even* integer is used for ``n``. The remainder ``r = remainder(x, + y``. If ``x / y`` is exactly halfway between two consecutive integers, the + nearest *even* integer is used for ``n``. The remainder ``r = remainder(x, y)`` thus always satisfies ``abs(r) <= 0.5 * abs(y)``. Special cases follow IEEE 754: in particular, ``remainder(x, math.inf)`` is - ``x`` for any finite ``x``, and ``remainder(x, 0)`` and - ``remainder(math.inf, x)`` raise :exc:`ValueError` for any non-NaN ``x``. + *x* for any finite *x*, and ``remainder(x, 0)`` and + ``remainder(math.inf, x)`` raise :exc:`ValueError` for any non-NaN *x*. If the result of the remainder operation is zero, that zero will have - the same sign as ``x``. + the same sign as *x*. On platforms using IEEE 754 binary floating-point, the result of this operation is always exactly representable: no rounding error is introduced. From 49aaa2afdbb0fcbdb193e64c9b69c9b68e1536e5 Mon Sep 17 00:00:00 2001 From: Mark Dickinson Date: Sat, 1 Apr 2017 16:37:04 +0100 Subject: [PATCH 3/6] Mark up function reference in what's new entry. --- Doc/whatsnew/3.7.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Doc/whatsnew/3.7.rst b/Doc/whatsnew/3.7.rst index a55215d17b5817..77212008eb0403 100644 --- a/Doc/whatsnew/3.7.rst +++ b/Doc/whatsnew/3.7.rst @@ -114,8 +114,8 @@ on Unix. (Contributed by Serhiy Storchaka in :issue:`25996`.) math ---- -New ``remainder`` function, implementing the IEEE 754-style remainder -operation. (Constructed by Mark Dickinson in :issue:`29962`.) +New :func:`~math.remainder` function, implementing the IEEE 754-style remainder +operation. (Contributed by Mark Dickinson in :issue:`29962`.) unittest.mock ------------- From 2b79f7793f1a001d5a256a99f68abb32e6adaee5 Mon Sep 17 00:00:00 2001 From: Mark Dickinson Date: Sat, 1 Apr 2017 16:46:15 +0100 Subject: [PATCH 4/6] Add comment explaining the calculation in the final branch. --- Modules/mathmodule.c | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c index 43bb0d1fa30b63..3b6db6e6b0839d 100644 --- a/Modules/mathmodule.c +++ b/Modules/mathmodule.c @@ -627,6 +627,33 @@ m_remainder(double x, double y) r = -c; } else { + /* + Here absx is exactly halfway between two multiples of absy, + and we need to choose the even multiple. x now has the form + + absx = n * absy + m + + for some integer n (recalling that m = 0.5*absy at this point). + If n is even we want to return m; if n is odd, we need to + return -m. + + So + + 0.5 * (absx - m) = (n/2) * absy + + and now reducing modulo absy gives us: + + | m, if n is odd + fmod(0.5 * (absx - m), absy) = | + | 0, if n is even + + Now m - 2.0 * fmod(...) gives the desired result: m + if n is even, -m if m is odd. + + Note that all steps in fmod(0.5 * (absx - m), absy) + will be computed exactly, with no rounding error + introduced. + */ assert(m == c); r = m - 2.0 * fmod(0.5 * (absx - m), absy); } From 6ab3fd3315e73db27fc296b5f2290972da34df51 Mon Sep 17 00:00:00 2001 From: Mark Dickinson Date: Sat, 1 Apr 2017 17:22:35 +0100 Subject: [PATCH 5/6] Fix out-of-order entry in whatsnew. --- Doc/whatsnew/3.7.rst | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/Doc/whatsnew/3.7.rst b/Doc/whatsnew/3.7.rst index 77212008eb0403..2bdfa55df17866 100644 --- a/Doc/whatsnew/3.7.rst +++ b/Doc/whatsnew/3.7.rst @@ -102,6 +102,12 @@ Added another argument *monetary* in :meth:`format_string` of :mod:`locale`. If *monetary* is true, the conversion uses monetary thousands separator and grouping strings. (Contributed by Garvit in :issue:`10379`.) +math +---- + +New :func:`~math.remainder` function, implementing the IEEE 754-style remainder +operation. (Contributed by Mark Dickinson in :issue:`29962`.) + os -- @@ -111,12 +117,6 @@ Serhiy Storchaka in :issue:`28682`.) Added support for :ref:`file descriptors ` in :func:`~os.scandir` on Unix. (Contributed by Serhiy Storchaka in :issue:`25996`.) -math ----- - -New :func:`~math.remainder` function, implementing the IEEE 754-style remainder -operation. (Contributed by Mark Dickinson in :issue:`29962`.) - unittest.mock ------------- From 6d80fceddc48a97f88fc783582092705fb051317 Mon Sep 17 00:00:00 2001 From: Mark Dickinson Date: Sun, 2 Apr 2017 11:43:53 +0100 Subject: [PATCH 6/6] Add comment explaining why it's good enough to compare m with c, in spite of possible rounding error. --- Modules/mathmodule.c | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c index 3b6db6e6b0839d..c19620d6bdc3f7 100644 --- a/Modules/mathmodule.c +++ b/Modules/mathmodule.c @@ -619,6 +619,27 @@ m_remainder(double x, double y) absx = fabs(x); absy = fabs(y); m = fmod(absx, absy); + + /* + Warning: some subtlety here. What we *want* to know at this point is + whether the remainder m is less than, equal to, or greater than half + of absy. However, we can't do that comparison directly because we + can't be sure that 0.5*absy is representable (the mutiplication + might incur precision loss due to underflow). So instead we compare + m with the complement c = absy - m: m < 0.5*absy if and only if m < + c, and so on. The catch is that absy - m might also not be + representable, but it turns out that it doesn't matter: + + - if m > 0.5*absy then absy - m is exactly representable, by + Sterbenz's lemma, so m > c + - if m == 0.5*absy then again absy - m is exactly representable + and m == c + - if m < 0.5*absy then either (i) 0.5*absy is exactly representable, + in which case 0.5*absy < absy - m, so 0.5*absy <= c and hence m < + c, or (ii) absy is tiny, either subnormal or in the lowest normal + binade. Then absy - m is exactly representable and again m < c. + */ + c = absy - m; if (m < c) { r = m;