diff options
| -rw-r--r-- | Doc/library/math.rst | 21 | ||||
| -rw-r--r-- | Doc/whatsnew/3.7.rst | 6 | ||||
| -rw-r--r-- | Lib/test/test_math.py | 135 | ||||
| -rw-r--r-- | Misc/NEWS | 3 | ||||
| -rw-r--r-- | Modules/mathmodule.c | 103 |
5 files changed, 268 insertions, 0 deletions
diff --git a/Doc/library/math.rst b/Doc/library/math.rst index b1f5692..0c53b29 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 12f65ff..b5107ea 100644 --- a/Doc/whatsnew/3.7.rst +++ b/Doc/whatsnew/3.7.rst @@ -110,6 +110,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 -- diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py index eaa41bc..70cb574 100644 --- a/Lib/test/test_math.py +++ b/Lib/test/test_math.py @@ -1000,6 +1000,135 @@ class MathTests(unittest.TestCase): 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 @@ class MathTests(unittest.TestCase): 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 @@ -303,6 +303,9 @@ Extension Modules Library ------- +- bpo-29962: Add math.remainder operation, implementing remainder + as specified in IEEE 754. + - bpo-29649: Improve struct.pack_into() exception messages for problems with the buffer size and offset. Patch by Andrew Nester. diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c index 243560a..c19620d 100644 --- a/Modules/mathmodule.c +++ b/Modules/mathmodule.c @@ -600,6 +600,102 @@ 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); + + /* + 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; + } + else if (m > c) { + 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); + } + 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 +1168,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 +2360,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}, |