rounding errors in Python floor division

Question

I know rounding errors happen in floating point arithmetic but can somebody explain the reason for this one:

>>> 8.0 / 0.4  # as expected 20.0 >>> floor(8.0 / 0.4)  # int works too 20 >>> 8.0 // 0.4  # expecting 20.0 19.0 

This happens on both Python 2 and 3 on x64.

As far as I see it this is either a bug or a very dumb specification of // since I don't see any reason why the last expression should evaluate to 19.0.

Why isn't a // b simply defined as floor(a / b) ?

EDIT: 8.0 % 0.4 also evaluates to 0.3999999999999996. At least this is consequent since then 8.0 // 0.4 * 0.4 + 8.0 % 0.4 evaluates to 8.0

EDIT: This is not a duplicate of Is floating point math broken? since I am asking why this specific operation is subject to (maybe avoidable) rounding errors, and why a // b isn't defined as / equal to floor(a / b)

REMARK: I guess that the deeper reason why this doesn't work is that floor division is discontinuous and thus has an infinite condition number making it an ill-posed problem. Floor division and floating-point numbers simply are fundamentally incompatible and you should never use // on floats. Just use integers or fractions instead.

Answer 1

As you and khelwood already noticed, 0.4 cannot be exactly represented as a float. Why? It is two fifth (4/10 == 2/5) which does not have a finite binary fraction representation.

Try this:

from fractions import Fraction Fraction('8.0') // Fraction('0.4')     # or equivalently     #     Fraction(8, 1) // Fraction(2, 5)     # or     #     Fraction('8/1') // Fraction('2/5') # 20 

However

Fraction('8') // Fraction(0.4) # 19 

Here, 0.4 is interpreted as a float literal (and thus a floating point binary number) which requires (binary) rounding, and only then converted to the rational number Fraction(3602879701896397, 9007199254740992), which is almost but not exactly 4 / 10. Then the floored division is executed, and because

19 * Fraction(3602879701896397, 9007199254740992) < 8.0 

and

20 * Fraction(3602879701896397, 9007199254740992) > 8.0 

the result is 19, not 20.

The same probably happens for

8.0 // 0.4 

I.e., it seems floored division is determined atomically (but on the only approximate float values of the interpreted float literals).

So why does

floor(8.0 / 0.4) 

give the "right" result? Because there, two rounding errors cancel each other out. First ¹⁾ the division is performed, yielding something slightly smaller than 20.0, but not representable as float. It gets rounded to the closest float, which happens to be 20.0. Only then, the floor operation is performed, but now acting on exactly 20.0, thus not changing the number any more.

---

¹⁾ As Kyle Strand points out, that the exact result is determined then rounded isn't what actually happens low²⁾-level (CPython's C code or even CPU instructions). However, it can be a useful model for determining the expected ³⁾ result.

²⁾ On the lowest ⁴⁾ level, however, this might not be too far off. Some chipsets determine float results by first computing a more precise (but still not exact, simply has some more binary digits) internal floating point result and then rounding to IEEE double precision.

³⁾ "expected" by the Python specification, not necessarily by our intuition.

⁴⁾ Well, lowest level above logic gates. We don't have to consider the quantum mechanics that make semiconductors possible to understand this.

Answer 2

After checking the semi-official sources of the float object in cpython on github (https://github.com/python/cpython/blob/966b24071af1b320a1c7646d33474eeae057c20f/Objects/floatobject.c) one can understand what happens here.

For normal division float_div is called (line 560) which internally converts the python floats to c-doubles, does the division and then converts the resulting double back to a python float. If you simply do that with 8.0/0.4 in c you get:

#include "stdio.h" #include "math.h"  int main(){     double vx = 8.0;     double wx = 0.4;     printf("%lf\n", floor(vx/wx));     printf("%d\n", (int)(floor(vx/wx))); }  // gives: // 20.000000 // 20 

For the floor division, something else happens. Internally, float_floor_div (line 654) gets called, which then calls float_divmod, a function that is supposed to return a tuple of python floats containing the floored division, as well as the mod/remainder, even though the latter is just thrown away by PyTuple_GET_ITEM(t, 0). These values are computed the following way (After conversion to c-doubles):

1. The remainder is computed by using double mod = fmod(numerator, denominator).
2. The numerator is reduced by mod to get a integral value when you then do the division.
3. The result for the floored division is calculated by effectively computing floor((numerator - mod) / denominator)
4. Afterwards, the check already mentioned in @Kasramvd's answer is done. But this only snaps the result of (numerator - mod) / denominator to the nearest integral value.

The reason why this gives a different result is, that fmod(8.0, 0.4) due to floating-point arithmetic gives 0.4 instead of 0.0. Therefore, the result that is computed is actually floor((8.0 - 0.4) / 0.4) = 19 and snapping (8.0 - 0.4) / 0.4) = 19 to the nearest integral value does not fix the error made introduced by the "wrong" result of fmod. You can easily chack that in c as well:

#include "stdio.h" #include "math.h"  int main(){     double vx = 8.0;     double wx = 0.4;     double mod = fmod(vx, wx);     printf("%lf\n", mod);     double div = (vx-mod)/wx;     printf("%lf\n", div); }  // gives: // 0.4 // 19.000000 

I would guess, that they chose this way of computing the floored division to keep the validity of (numerator//divisor)*divisor + fmod(numerator, divisor) = numerator (as mentioned in the link in @0x539's answer), even though this now results in a somewhat unexpected behavior of floor(8.0/0.4) != 8.0//0.4.