bpo-29962: add math.remainder (#950)
* Implement math.remainder.
* Fix markup for arguments; use double spaces after period.
* Mark up function reference in what's new entry.
* Add comment explaining the calculation in the final branch.
* Fix out-of-order entry in whatsnew.
* Add comment explaining why it's good enough to compare m with c, in spite of possible rounding error.
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 @@
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 @@
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 @@
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 @@
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 876767e..396ada1 100644
--- a/Misc/NEWS
+++ b/Misc/NEWS
@@ -303,6 +303,9 @@
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 @@
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 @@
"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 @@
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},