Skip to content

Commit a0ce375

Browse files
authored
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.
1 parent a0157b5 commit a0ce375

File tree

5 files changed

+268
-0
lines changed

5 files changed

+268
-0
lines changed

Doc/library/math.rst

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -175,6 +175,27 @@ Number-theoretic and representation functions
175175
of *x* and are floats.
176176

177177

178+
.. function:: remainder(x, y)
179+
180+
Return the IEEE 754-style remainder of *x* with respect to *y*. For
181+
finite *x* and finite nonzero *y*, this is the difference ``x - n*y``,
182+
where ``n`` is the closest integer to the exact value of the quotient ``x /
183+
y``. If ``x / y`` is exactly halfway between two consecutive integers, the
184+
nearest *even* integer is used for ``n``. The remainder ``r = remainder(x,
185+
y)`` thus always satisfies ``abs(r) <= 0.5 * abs(y)``.
186+
187+
Special cases follow IEEE 754: in particular, ``remainder(x, math.inf)`` is
188+
*x* for any finite *x*, and ``remainder(x, 0)`` and
189+
``remainder(math.inf, x)`` raise :exc:`ValueError` for any non-NaN *x*.
190+
If the result of the remainder operation is zero, that zero will have
191+
the same sign as *x*.
192+
193+
On platforms using IEEE 754 binary floating-point, the result of this
194+
operation is always exactly representable: no rounding error is introduced.
195+
196+
.. versionadded:: 3.7
197+
198+
178199
.. function:: trunc(x)
179200

180201
Return the :class:`~numbers.Real` value *x* truncated to an

Doc/whatsnew/3.7.rst

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -110,6 +110,12 @@ Added another argument *monetary* in :meth:`format_string` of :mod:`locale`.
110110
If *monetary* is true, the conversion uses monetary thousands separator and
111111
grouping strings. (Contributed by Garvit in :issue:`10379`.)
112112

113+
math
114+
----
115+
116+
New :func:`~math.remainder` function, implementing the IEEE 754-style remainder
117+
operation. (Contributed by Mark Dickinson in :issue:`29962`.)
118+
113119
os
114120
--
115121

Lib/test/test_math.py

Lines changed: 135 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1000,6 +1000,135 @@ def testRadians(self):
10001000
self.ftest('radians(-45)', math.radians(-45), -math.pi/4)
10011001
self.ftest('radians(0)', math.radians(0), 0)
10021002

1003+
@requires_IEEE_754
1004+
def testRemainder(self):
1005+
from fractions import Fraction
1006+
1007+
def validate_spec(x, y, r):
1008+
"""
1009+
Check that r matches remainder(x, y) according to the IEEE 754
1010+
specification. Assumes that x, y and r are finite and y is nonzero.
1011+
"""
1012+
fx, fy, fr = Fraction(x), Fraction(y), Fraction(r)
1013+
# r should not exceed y/2 in absolute value
1014+
self.assertLessEqual(abs(fr), abs(fy/2))
1015+
# x - r should be an exact integer multiple of y
1016+
n = (fx - fr) / fy
1017+
self.assertEqual(n, int(n))
1018+
if abs(fr) == abs(fy/2):
1019+
# If |r| == |y/2|, n should be even.
1020+
self.assertEqual(n/2, int(n/2))
1021+
1022+
# triples (x, y, remainder(x, y)) in hexadecimal form.
1023+
testcases = [
1024+
# Remainders modulo 1, showing the ties-to-even behaviour.
1025+
'-4.0 1 -0.0',
1026+
'-3.8 1 0.8',
1027+
'-3.0 1 -0.0',
1028+
'-2.8 1 -0.8',
1029+
'-2.0 1 -0.0',
1030+
'-1.8 1 0.8',
1031+
'-1.0 1 -0.0',
1032+
'-0.8 1 -0.8',
1033+
'-0.0 1 -0.0',
1034+
' 0.0 1 0.0',
1035+
' 0.8 1 0.8',
1036+
' 1.0 1 0.0',
1037+
' 1.8 1 -0.8',
1038+
' 2.0 1 0.0',
1039+
' 2.8 1 0.8',
1040+
' 3.0 1 0.0',
1041+
' 3.8 1 -0.8',
1042+
' 4.0 1 0.0',
1043+
1044+
# Reductions modulo 2*pi
1045+
'0x0.0p+0 0x1.921fb54442d18p+2 0x0.0p+0',
1046+
'0x1.921fb54442d18p+0 0x1.921fb54442d18p+2 0x1.921fb54442d18p+0',
1047+
'0x1.921fb54442d17p+1 0x1.921fb54442d18p+2 0x1.921fb54442d17p+1',
1048+
'0x1.921fb54442d18p+1 0x1.921fb54442d18p+2 0x1.921fb54442d18p+1',
1049+
'0x1.921fb54442d19p+1 0x1.921fb54442d18p+2 -0x1.921fb54442d17p+1',
1050+
'0x1.921fb54442d17p+2 0x1.921fb54442d18p+2 -0x0.0000000000001p+2',
1051+
'0x1.921fb54442d18p+2 0x1.921fb54442d18p+2 0x0p0',
1052+
'0x1.921fb54442d19p+2 0x1.921fb54442d18p+2 0x0.0000000000001p+2',
1053+
'0x1.2d97c7f3321d1p+3 0x1.921fb54442d18p+2 0x1.921fb54442d14p+1',
1054+
'0x1.2d97c7f3321d2p+3 0x1.921fb54442d18p+2 -0x1.921fb54442d18p+1',
1055+
'0x1.2d97c7f3321d3p+3 0x1.921fb54442d18p+2 -0x1.921fb54442d14p+1',
1056+
'0x1.921fb54442d17p+3 0x1.921fb54442d18p+2 -0x0.0000000000001p+3',
1057+
'0x1.921fb54442d18p+3 0x1.921fb54442d18p+2 0x0p0',
1058+
'0x1.921fb54442d19p+3 0x1.921fb54442d18p+2 0x0.0000000000001p+3',
1059+
'0x1.f6a7a2955385dp+3 0x1.921fb54442d18p+2 0x1.921fb54442d14p+1',
1060+
'0x1.f6a7a2955385ep+3 0x1.921fb54442d18p+2 0x1.921fb54442d18p+1',
1061+
'0x1.f6a7a2955385fp+3 0x1.921fb54442d18p+2 -0x1.921fb54442d14p+1',
1062+
'0x1.1475cc9eedf00p+5 0x1.921fb54442d18p+2 0x1.921fb54442d10p+1',
1063+
'0x1.1475cc9eedf01p+5 0x1.921fb54442d18p+2 -0x1.921fb54442d10p+1',
1064+
1065+
# Symmetry with respect to signs.
1066+
' 1 0.c 0.4',
1067+
'-1 0.c -0.4',
1068+
' 1 -0.c 0.4',
1069+
'-1 -0.c -0.4',
1070+
' 1.4 0.c -0.4',
1071+
'-1.4 0.c 0.4',
1072+
' 1.4 -0.c -0.4',
1073+
'-1.4 -0.c 0.4',
1074+
1075+
# Huge modulus, to check that the underlying algorithm doesn't
1076+
# rely on 2.0 * modulus being representable.
1077+
'0x1.dp+1023 0x1.4p+1023 0x0.9p+1023',
1078+
'0x1.ep+1023 0x1.4p+1023 -0x0.ap+1023',
1079+
'0x1.fp+1023 0x1.4p+1023 -0x0.9p+1023',
1080+
]
1081+
1082+
for case in testcases:
1083+
with self.subTest(case=case):
1084+
x_hex, y_hex, expected_hex = case.split()
1085+
x = float.fromhex(x_hex)
1086+
y = float.fromhex(y_hex)
1087+
expected = float.fromhex(expected_hex)
1088+
validate_spec(x, y, expected)
1089+
actual = math.remainder(x, y)
1090+
# Cheap way of checking that the floats are
1091+
# as identical as we need them to be.
1092+
self.assertEqual(actual.hex(), expected.hex())
1093+
1094+
# Test tiny subnormal modulus: there's potential for
1095+
# getting the implementation wrong here (for example,
1096+
# by assuming that modulus/2 is exactly representable).
1097+
tiny = float.fromhex('1p-1074') # min +ve subnormal
1098+
for n in range(-25, 25):
1099+
if n == 0:
1100+
continue
1101+
y = n * tiny
1102+
for m in range(100):
1103+
x = m * tiny
1104+
actual = math.remainder(x, y)
1105+
validate_spec(x, y, actual)
1106+
actual = math.remainder(-x, y)
1107+
validate_spec(-x, y, actual)
1108+
1109+
# Special values.
1110+
# NaNs should propagate as usual.
1111+
for value in [NAN, 0.0, -0.0, 2.0, -2.3, NINF, INF]:
1112+
self.assertIsNaN(math.remainder(NAN, value))
1113+
self.assertIsNaN(math.remainder(value, NAN))
1114+
1115+
# remainder(x, inf) is x, for non-nan non-infinite x.
1116+
for value in [-2.3, -0.0, 0.0, 2.3]:
1117+
self.assertEqual(math.remainder(value, INF), value)
1118+
self.assertEqual(math.remainder(value, NINF), value)
1119+
1120+
# remainder(x, 0) and remainder(infinity, x) for non-NaN x are invalid
1121+
# operations according to IEEE 754-2008 7.2(f), and should raise.
1122+
for value in [NINF, -2.3, -0.0, 0.0, 2.3, INF]:
1123+
with self.assertRaises(ValueError):
1124+
math.remainder(INF, value)
1125+
with self.assertRaises(ValueError):
1126+
math.remainder(NINF, value)
1127+
with self.assertRaises(ValueError):
1128+
math.remainder(value, 0.0)
1129+
with self.assertRaises(ValueError):
1130+
math.remainder(value, -0.0)
1131+
10031132
def testSin(self):
10041133
self.assertRaises(TypeError, math.sin)
10051134
self.ftest('sin(0)', math.sin(0), 0)
@@ -1286,6 +1415,12 @@ def test_mtestfile(self):
12861415
self.fail('Failures in test_mtestfile:\n ' +
12871416
'\n '.join(failures))
12881417

1418+
# Custom assertions.
1419+
1420+
def assertIsNaN(self, value):
1421+
if not math.isnan(value):
1422+
self.fail("Expected a NaN, got {!r}.".format(value))
1423+
12891424

12901425
class IsCloseTests(unittest.TestCase):
12911426
isclose = math.isclose # sublcasses should override this

Misc/NEWS

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -303,6 +303,9 @@ Extension Modules
303303
Library
304304
-------
305305

306+
- bpo-29962: Add math.remainder operation, implementing remainder
307+
as specified in IEEE 754.
308+
306309
- bpo-29649: Improve struct.pack_into() exception messages for problems
307310
with the buffer size and offset. Patch by Andrew Nester.
308311

Modules/mathmodule.c

Lines changed: 103 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -600,6 +600,102 @@ m_atan2(double y, double x)
600600
return atan2(y, x);
601601
}
602602

603+
604+
/* IEEE 754-style remainder operation: x - n*y where n*y is the nearest
605+
multiple of y to x, taking n even in the case of a tie. Assuming an IEEE 754
606+
binary floating-point format, the result is always exact. */
607+
608+
static double
609+
m_remainder(double x, double y)
610+
{
611+
/* Deal with most common case first. */
612+
if (Py_IS_FINITE(x) && Py_IS_FINITE(y)) {
613+
double absx, absy, c, m, r;
614+
615+
if (y == 0.0) {
616+
return Py_NAN;
617+
}
618+
619+
absx = fabs(x);
620+
absy = fabs(y);
621+
m = fmod(absx, absy);
622+
623+
/*
624+
Warning: some subtlety here. What we *want* to know at this point is
625+
whether the remainder m is less than, equal to, or greater than half
626+
of absy. However, we can't do that comparison directly because we
627+
can't be sure that 0.5*absy is representable (the mutiplication
628+
might incur precision loss due to underflow). So instead we compare
629+
m with the complement c = absy - m: m < 0.5*absy if and only if m <
630+
c, and so on. The catch is that absy - m might also not be
631+
representable, but it turns out that it doesn't matter:
632+
633+
- if m > 0.5*absy then absy - m is exactly representable, by
634+
Sterbenz's lemma, so m > c
635+
- if m == 0.5*absy then again absy - m is exactly representable
636+
and m == c
637+
- if m < 0.5*absy then either (i) 0.5*absy is exactly representable,
638+
in which case 0.5*absy < absy - m, so 0.5*absy <= c and hence m <
639+
c, or (ii) absy is tiny, either subnormal or in the lowest normal
640+
binade. Then absy - m is exactly representable and again m < c.
641+
*/
642+
643+
c = absy - m;
644+
if (m < c) {
645+
r = m;
646+
}
647+
else if (m > c) {
648+
r = -c;
649+
}
650+
else {
651+
/*
652+
Here absx is exactly halfway between two multiples of absy,
653+
and we need to choose the even multiple. x now has the form
654+
655+
absx = n * absy + m
656+
657+
for some integer n (recalling that m = 0.5*absy at this point).
658+
If n is even we want to return m; if n is odd, we need to
659+
return -m.
660+
661+
So
662+
663+
0.5 * (absx - m) = (n/2) * absy
664+
665+
and now reducing modulo absy gives us:
666+
667+
| m, if n is odd
668+
fmod(0.5 * (absx - m), absy) = |
669+
| 0, if n is even
670+
671+
Now m - 2.0 * fmod(...) gives the desired result: m
672+
if n is even, -m if m is odd.
673+
674+
Note that all steps in fmod(0.5 * (absx - m), absy)
675+
will be computed exactly, with no rounding error
676+
introduced.
677+
*/
678+
assert(m == c);
679+
r = m - 2.0 * fmod(0.5 * (absx - m), absy);
680+
}
681+
return copysign(1.0, x) * r;
682+
}
683+
684+
/* Special values. */
685+
if (Py_IS_NAN(x)) {
686+
return x;
687+
}
688+
if (Py_IS_NAN(y)) {
689+
return y;
690+
}
691+
if (Py_IS_INFINITY(x)) {
692+
return Py_NAN;
693+
}
694+
assert(Py_IS_INFINITY(y));
695+
return x;
696+
}
697+
698+
603699
/*
604700
Various platforms (Solaris, OpenBSD) do nonstandard things for log(0),
605701
log(-ve), log(NaN). Here are wrappers for log and log10 that deal with
@@ -1072,6 +1168,12 @@ FUNC1(log1p, m_log1p, 0,
10721168
"log1p($module, x, /)\n--\n\n"
10731169
"Return the natural logarithm of 1+x (base e).\n\n"
10741170
"The result is computed in a way which is accurate for x near zero.")
1171+
FUNC2(remainder, m_remainder,
1172+
"remainder($module, x, y, /)\n--\n\n"
1173+
"Difference between x and the closest integer multiple of y.\n\n"
1174+
"Return x - n*y where n*y is the closest integer multiple of y.\n"
1175+
"In the case where x is exactly halfway between two multiples of\n"
1176+
"y, the nearest even value of n is used. The result is always exact.")
10751177
FUNC1(sin, sin, 0,
10761178
"sin($module, x, /)\n--\n\n"
10771179
"Return the sine of x (measured in radians).")
@@ -2258,6 +2360,7 @@ static PyMethodDef math_methods[] = {
22582360
MATH_MODF_METHODDEF
22592361
MATH_POW_METHODDEF
22602362
MATH_RADIANS_METHODDEF
2363+
{"remainder", math_remainder, METH_VARARGS, math_remainder_doc},
22612364
{"sin", math_sin, METH_O, math_sin_doc},
22622365
{"sinh", math_sinh, METH_O, math_sinh_doc},
22632366
{"sqrt", math_sqrt, METH_O, math_sqrt_doc},

0 commit comments

Comments
 (0)