Python地板除法中的舍入错误

Question

我知道浮点运算中会发生舍入错误，但有人可以解释一下原因：

>>> 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

这发生在 x64 上的 Python 2 和 3 上。

在我看来，这要么是一个错误，要么是 // 的一个非常愚蠢的规范，因为我看不出最后一个表达式应该评估为 19.0 的任何理由。

为什么a // b 不简单定义为floor(a / b)？

编辑：8.0 % 0.4 也计算为0.3999999999999996。至少这是随之而来的，因为那时8.0 // 0.4 * 0.4 + 8.0 % 0.4 评估为8.0

编辑：这不是 Is floating point math broken? 的重复，因为我在问为什么这个特定的操作会出现（也许可以避免的）舍入错误，以及为什么 a // b 没有被定义为 /等于floor(a / b)

REMARK：我想这不起作用的更深层原因是楼层划分是不连续的，因此有一个无限的condition number，使其成为一个不适定问题。楼层除法和浮点数根本不兼容，你永远不应该在浮点数上使用//。只需使用整数或分数即可。

Answer 1

正如您和 khelwood 已经注意到的，0.4 不能完全表示为浮点数。为什么？它是五分之二 (4/10 == 2/5)，它没有有限二进制分数表示。

试试这个：

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

然而

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

这里，0.4 被解释为浮点字面量（因此是浮点二进制数），需要（二进制）舍入，并且只有 then 转换为有理数 Fraction(3602879701896397, 9007199254740992)，这几乎但不完全是 4 / 10。然后执行地板除法，因为

19 * Fraction(3602879701896397, 9007199254740992) < 8.0

和

20 * Fraction(3602879701896397, 9007199254740992) > 8.0

结果是 19，而不是 20。

同样的情况也可能发生在

8.0 // 0.4

也就是说，似乎下限除法是原子确定的（但在解释的浮点文字的唯一近似浮点值上）。

为什么会这样

floor(8.0 / 0.4)

给出“正确”的结果？因为在那里，两个舍入误差相互抵消。 首先 ¹⁾ 执行除法，产生的结果略小于 20.0，但不能表示为浮点数。它被四舍五入到最接近的浮点数，恰好是20.0。只有然后，执行了floor 操作，但现在完全作用于20.0，因此不再更改数字。

---

¹⁾ 作为 Kyle Strand points out，确定了确切的结果然后四舍五入不是实际上 /em> 发生在²⁾ 级别（CPython 的 C 代码甚至 CPU 指令）。但是，它可能是确定预期 ³⁾ 结果的有用模型。

²⁾ 但是，在最低 ⁴⁾ 级别上，这可能不会太远。一些芯片组通过首先计算更精确（但仍然不精确，只是有更多二进制数字）的内部浮点结果然后四舍五入到 IEEE 双精度来确定浮点结果。

³⁾ Python 规范“预期”的，不一定是我们的直觉。

⁴⁾ 嗯，最低级别高于逻辑门。我们不必考虑使半导体能够理解这一点的量子力学。

Answer 2

在 github (https://github.com/python/cpython/blob/966b24071af1b320a1c7646d33474eeae057c20f/Objects/floatobject.c) 上查看了 cpython 中 float 对象的半官方来源后，可以理解这里发生了什么。

对于普通除法float_div 被调用（第 560 行），它在内部将 python floats 转换为 c-doubles，进行除法，然后将生成的 double 转换回 python float .如果您只是在 c 中使用 8.0/0.4 来执行此操作，您会得到：

#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

对于楼层划分，发生了其他事情。在内部，float_floor_div（第 654 行）被调用，然后它调用float_divmod，该函数应该返回一个 python 元组 floats 包含地板除法，以及 mod/remainder，即使后者被PyTuple_GET_ITEM(t, 0) 扔掉了。这些值按以下方式计算（转换为 c-doubles 后）：

1. 使用double mod = fmod(numerator, denominator)计算余数。
2. 在进行除法时，分子会减少 mod 以获得整数值。
3. 通过有效计算floor((numerator - mod) / denominator) 来计算下除法的结果
4. 之后，@Kasramvd 的回答中已经提到的检查完成。但这只会将(numerator - mod) / denominator 的结果捕捉到最接近的整数值。

这给出不同结果的原因是，fmod(8.0, 0.4) 由于浮点运算给出了0.4 而不是0.0。因此，计算的结果实际上是floor((8.0 - 0.4) / 0.4) = 19，将(8.0 - 0.4) / 0.4) = 19 捕捉到最接近的整数值并不能解决由fmod 的“错误”结果引入的错误。你也可以在 c 中轻松地破解它：

#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

我猜，他们选择了这种计算底除法的方式来保持(numerator//divisor)*divisor + fmod(numerator, divisor) = numerator 的有效性（如@0x539 答案中的链接中所述），即使现在这会导致@987654348 的一些意外行为@。

Answer 3

@jotasi 解释了背后的真正原因。

但是，如果您想阻止它，您可以使用 decimal 模块，该模块基本上旨在表示十进制浮点数，与二进制浮点表示完全相反。

因此，在您的情况下，您可以执行以下操作：

>>> from decimal import *
>>> Decimal('8.0')//Decimal('0.4')
Decimal('20')

参考：https://docs.python.org/2/library/decimal.html

Answer 4

好的，经过一些研究，我发现了这个issue。 似乎正在发生的事情是，正如@khelwood 建议的那样，0.4 在内部评估为0.40000000000000002220，当除以8.0 时，产生的结果略小于20.0。然后/ 运算符舍入到最接近的浮点数，即20.0，但// 运算符立即截断结果，产生19.0。

这应该更快，我认为它“接近处理器”，但我仍然不是用户想要/期望的。

Answer 5

那是因为 python 中没有 0.4（浮点有限表示），它实际上是一个像 0.4000000000000001 这样的浮点数，它使得除法的下限为 19。

>>> floor(8//0.4000000000000001)
19.0

但真正的除法 (/) returns a reasonable approximation of the division result if the arguments are floats or complex. 这就是为什么 8.0/0.4 的结果是 20。它实际上取决于参数的大小（在 C 双参数中）。 （不四舍五入到最近的浮点数）

阅读 Guido 本人关于 pythons integer division floors 的更多信息。

有关浮点数的完整信息，您可以阅读这篇文章https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

有兴趣的朋友可以看一下Cpython源码中的float_div，它对浮点数做真正的除法任务：

float_div(PyObject *v, PyObject *w)
{
    double a,b;
    CONVERT_TO_DOUBLE(v, a);
    CONVERT_TO_DOUBLE(w, b);
    if (b == 0.0) {
        PyErr_SetString(PyExc_ZeroDivisionError,
                        "float division by zero");
        return NULL;
    }
    PyFPE_START_PROTECT("divide", return 0)
    a = a / b;
    PyFPE_END_PROTECT(a)
    return PyFloat_FromDouble(a);
}

最终结果将由函数PyFloat_FromDouble计算：

PyFloat_FromDouble(double fval)
{
    PyFloatObject *op = free_list;
    if (op != NULL) {
        free_list = (PyFloatObject *) Py_TYPE(op);
        numfree--;
    } else {
        op = (PyFloatObject*) PyObject_MALLOC(sizeof(PyFloatObject));
        if (!op)
            return PyErr_NoMemory();
    }
    /* Inline PyObject_New */
    (void)PyObject_INIT(op, &PyFloat_Type);
    op->ob_fval = fval;
    return (PyObject *) op;
}