float.__mod__ is unstable: if n % k == k then n % k % k == 0

[cls-oi]

Bug report

Bug description:

As an example:

>>> import math
>>> n = math.atan2(-1.0, 4503599627370496.0)
>>> k = 2 * math.pi
>>> n
-2.220446049250313e-16
>>> k
6.283185307179586
>>> n % k
6.283185307179586
>>> n % k % k
0.0
>>>

This looks to be due to:

cpython/Objects/floatobject.c

Lines 619 to 625 in 569fc68

	mod = fmod(vx, wx);
	if (mod) {
	/* ensure the remainder has the same sign as the denominator */
	if ((wx < 0) != (mod < 0)) {
	mod += wx;
	}
	}

On the first modulo we have mod < 0 but mod + wx == wx, so n % k == k, and then on the second modulo we have fmod(wx, wx) == 0, so n % k % k == 0.

CPython versions tested on:

3.12

Operating systems tested on:

Linux

[picnixz]

This would be correct because of rounding. Even if n != 0, we have n + k == k, which is why you get n % k == k (here the mod would have the same sign as k) (and then k % k == copysign(0, k))

[skirpichev]

Documentation warns about this: "Return the floating-point remainder of x / y, as defined by the platform C library function fmod(x, y). Note that the Python expression x % y may not return the same result. The intent of the C standard is that fmod(x, y) be exactly (mathematically; to infinite precision) equal to x - n*y for some integer n such that the result has the same sign as x and magnitude less than abs(y). Python’s x % y returns a result with the sign of y instead, and may not be exactly computable for float arguments. For example, fmod(-1e-100, 1e100) is -1e-100, but the result of Python’s -1e-100 % 1e100 is 1e100-1e-100, which cannot be represented exactly as a float, and rounds to the surprising 1e100. For this reason, function fmod() is generally preferred when working with floats, while Python’s x % y is preferred when working with integers."

I don't think it's fixable if we intent to keep property of being same sign with denominator.

[cls-oi]

Thanks for the citation - it's good to know that this is a known issue, but I'm still confused by the specific choice of behaviour given that there's a simple alternative that would avoid this pitfall -

I don't think it's fixable if we intent to keep property of being same sign with denominator.

I would understand if n % k == k weren't fixable because n and n % k must have the same sign, if it weren't for the fact that n % k % k must also then have the same sign, and n % k % k == 0.0.

The alternative would be if, in the example above, n % k instead returned 0.0. In that case we'd have n % k == n % k % k and also the desirable property that n % k < k.

This would be achievable, for example, by adding an additional check in the code above:

         /* ensure the remainder has the same sign as the denominator */
         if ((wx < 0) != (mod < 0)) {
             mod += wx;
+            if (mod == wx) {
+                mod = copysign(0.0, wx);
+            }
         }
     }
     else {

An alternative might be to change the existing initial mod check:

         return NULL;
     }
     mod = fmod(vx, wx);
-    if (mod) {
+    if (mod + wx != wx) {
         /* ensure the remainder has the same sign as the denominator */
         if ((wx < 0) != (mod < 0)) {
             mod += wx;

Is there any chance we could consider this alternative behaviour, given that the existing behaviour calls for an explicit warning in the documentation discouraging its use? I'd be happy to raise a PR for further discussion if there's any chance of its consideration.

[tim-one]

An invariant the implementation is trying hard to maintain is that after q, r = divmod(x, y), x == q*y + r as closely as possible.

>>> n = -2.220446049250313e-16
>>> k = 6.283185307179586
>>> q, r = divmod(n, k)
>>> n, k
-2.220446049250313e-16, 6.283185307179586)
>>> q, r
-1.0, 6.283185307179586)
>>> q * k + r
0.0

So the invariant is approximately true. % and floor division are intimately related in Python; you can't fiddle one independently of the other.

If we changed r to be 0, then to maintain the invariant we'd have to make n // k 0 too. But that's even worse: the floor of n/k is -1, to infinite precision. That's inarguable.

So something had to "give". It's not possible to simultaneously satisfy all desirable properties.

As the docs suggest, don't use % for floats if you can't live with this.

[cls-oi]

Thanks for the explanation. I suppose what I would expect is n // k == math.trunc(n / k) == 0, as with integer division in C - but I'm aware that isn't a change that's going to happen (and I'm sure has its own justifications), so I'll concede. I only wish %'s shortcomings were more apparent without having read the documentation for math.fmod, but I've no suggestions on that front.

Edit: I just found Why Python's Integer Division Floors, which also mentions:

Tim Peters, who knows where all Python's floating point skeletons are buried, has expressed some worry about my desire to extend these rules to floating point modulo. He's probably right; the truncate-towards-negative-infinity rule can cause precision loss for x%1.0 when x is a very small negative number. But that's not enough for me to break integer modulo, and // is tightly coupled to that.

Well, so it is, I suppose.

[tim-one]

Ya, as Guido wrote, using floor division for integers in this context is universally accepted by mathematicians. It saves a lot of pain. For example, if an analog clock says "1", what did it say 2 hours ago? -1 % 12 returning 11 is far more helpful than returning -1.

But if a language is determined to pretend that floats are a superset of ints (which they aren't, in computers), then -1.0 % 12.0 has to do the same.

BTW, the docs should be updated to mention the newer math.remainder(). That's what most people using floats really want. Like fmod(), it's also exact, but it doesn't care about the signs of its inputs. Instead it returns the difference from the closest integer multiple of the modulus, regardless of whether that's in the $-\infty$ or $+\infty$ direction.

>>> from math import remainder
>>> remainder(1.1, 1)
0.10000000000000009
>>> remainder(1.9, 1)
-0.10000000000000009
>>> remainder(1e-100, 1)
1e-100
>>> remainder(-1e-100, 1)
-1e-100

[cls-oi]

Ah, math.remainder does actually seem to be the best tool for the job. Thanks for the tip!