Ошибки округления в Python floor division

Ответ 1

Как вы и хелвуд уже заметили, 0.4 не может быть точно представлен как float. Зачем? Это две пятые (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 интерпретируется как плавающий литерал (и, следовательно, двоичное число с плавающей запятой), который требует (двоичного) округления и только затем преобразуется в рациональное число Fraction(3602879701896397, 9007199254740992), которое почти, но не точно 4/10. Затем выполняется разделение полов, и потому что

19 * Fraction(3602879701896397, 9007199254740992) < 8.0

и

20 * Fraction(3602879701896397, 9007199254740992) > 8.0

результат равен 19, а не 20.

То же самое происходит и для

8.0 // 0.4

I.e., по-видимому, разделение пополам определяется атомарно (но только по приблизительным значениям поплавка интерпретируемых плавающих литералов).

Так почему же

floor(8.0 / 0.4)

дать "правильный" результат? Потому что там две ошибки округления компенсируют друг друга. Сначала ¹⁾ выполняется деление, принося что-то немного меньше 20.0, но не представляемое как float. Он округляется до ближайшего поплавка, что бывает 20.0. Только тогда выполняется операция floor, но теперь она действует именно на 20.0, тем самым больше не изменяя номер.

¹⁾ Как указывает Кайл Странд , точный результат определяется, а округленный не что на самом деле происходит с низким уровнем ²⁾ -level (код CPython C или даже инструкции CPU). Однако это может быть полезной моделью для определения ожидаемого результата ³⁾.

²⁾ На самом низком уровне ⁴⁾ это может быть не слишком далеко. Некоторые чипсеты определяют результаты поплавка, сначала вычисляя более точный (но все же не точный, просто несколько более двоичных цифр) результат внутренней с плавающей запятой, а затем округляя до двойной точности IEEE.

³⁾ "ожидаемый" по спецификации Python, не обязательно нашей интуицией.

⁴⁾ Ну, самый низкий уровень выше логических ворот. Нам не нужно рассматривать квантовую механику, позволяющую понять полупроводники.

Ответ 2

@jotasi объяснил истинную причину этого.

Однако, если вы хотите его предотвратить, вы можете использовать модуль decimal, который был в основном предназначен для представления десятичных чисел с плавающей запятой точно в отличие от двоичного представления с плавающей запятой.

Итак, в вашем случае вы можете сделать что-то вроде:

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

Ссылка: https://docs.python.org/2/library/decimal.html

Ответ 3

Хорошо после небольшого исследования я нашел это issue. Кажется, что происходит, так как @hhelwood предположил, что 0.4 внутренне оценивается до 0.40000000000000002220, который при делении 8.0 дает нечто немного меньшее, чем 20.0. Оператор / затем округляется до ближайшего числа с плавающей запятой, которое равно 20.0, но оператор // немедленно обрезает результат, получая 19.0.

Это должно быть быстрее, и я полагаю, что это "близко к процессору", но я все еще не то, чего хочет/ожидает пользователь.

Ответ 4

Это потому, что в питоне нет (0.4-конечное представление с плавающей запятой), на самом деле это float как 0.4000000000000001, что делает пол деления равным 19.

>>> floor(8//0.4000000000000001)
19.0

Но истинное деление (/) возвращает разумную аппроксимацию результата деления, если аргументы являются float или сложными. И вот почему результат 8.0/0.4 равен 20. Фактически это зависит от размера аргументов (в двойных аргументах C). ( не округление до ближайшего плавающего)

Подробнее о питонах целых делений поэта самим Гвидо.

Также для полной информации о номерах с плавающей запятой вы можете прочитать эту статью https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

Для тех, у кого есть интерес, следующая функция - это float_div, которая выполняет истинную задачу разделения для чисел с плавающей запятой, в исходном коде Cpython:

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

Ответ 5

После проверки полуофициальных источников объекта float в cpython на github (https://github.com/python/cpython/blob/966b24071af1b320a1c7646d33474eeae057c20f/Objects/floatobject.c) можно понять, что здесь происходит.

Для нормального деления float_div вызывается (строка 560), который внутренне преобразует python float в c- double s, выполняет деление и затем преобразует полученный double обратно в python float. Если вы просто делаете это с 8.0/0.4 в c, вы получаете:

#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 float, содержащий разделенное на пол разделение, а также mod/остаток, хотя последний просто выбрасывается на PyTuple_GET_ITEM(t, 0). Эти значения вычисляются следующим образом (после преобразования в c- double s):

• Остальная часть вычисляется с помощью double mod = fmod(numerator, denominator).
• Числитель уменьшается на mod, чтобы получить интегральное значение, когда вы затем выполняете деление.
• Результат для разделенного полов вычисляется путем эффективного вычисления floor((numerator - mod) / denominator)
• Затем выполняется проверка, уже упомянутая в ответе @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), хотя это теперь приводит к несколько неожиданному поведению floor(8.0/0.4) != 8.0//0.4.