15. Арифметика с плавающей точкой: проблемы и ограничения

Числа с плавающей точкой представляются в компьютерном оборудовании в виде дробей по основанию 2 (двоичных). Например, десятичная дробь 0.625 имеет значение 6/10 + 2/100 + 5/1000, и точно так же двоичная дробь 0.101 имеет значение 1/2 + 0/4 + 1/8. Эти две дроби имеют одинаковые значения, единственное различие заключается в том, что первая записывается в дробной системе счисления по основанию 10, а вторая - по основанию 2.

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

Задачу проще понять сначала в базе 10. Рассмотрим дробь 1/3. Вы можете приблизительно представить ее в виде дроби по основанию 10:

0.3

или, лучше,

0.33

или, лучше,

0.333

и так далее. Сколько бы цифр вы ни записали, результат никогда не будет равен 1/3, но будет все более и более точным приближением к 1/3.

Таким же образом, независимо от того, сколько цифр основания 2 вы хотите использовать, десятичное значение 0,1 не может быть представлено в точности как дробь основания 2. В основании 2 1/10 - это бесконечно повторяющаяся дробь

0.0001100110011001100110011001100110011001100110011...

Остановитесь на любом конечном количестве битов, и вы получите приближение. На большинстве современных машин плавающие числа аппроксимируются с помощью двоичной дроби, в числителе которой используются первые 53 бита, начиная со старшего бита, а в знаменателе - сила двух. В случае 1/10 двоичная дробь равна 3602879701896397 / 2 ** 55, что близко к истинному значению 1/10, но не совсем равно ему.

Многие пользователи не знают об этом приближении из-за способа отображения значений. Python печатает только десятичное приближение к истинному десятичному значению двоичного приближения, хранящегося в памяти машины. На большинстве машин, если бы Python выводил истинное десятичное значение двоичного приближения, хранящегося для 0,1, он должен был бы отобразить:

>>> 0.1
0.1000000000000000055511151231257827021181583404541015625

Это больше цифр, чем большинство людей считает нужным, поэтому Python сохраняет количество цифр, выводя вместо них округленное значение:

>>> 1 / 10
0.1

Помните, что даже если распечатанный результат выглядит как точное значение 1/10, на самом деле хранящееся в памяти значение - это ближайшая представимая двоичная дробь.

Интересно, что существует множество различных десятичных чисел, которые имеют одну и ту же ближайшую приближенную двоичную дробь. Например, числа 0.1, 0.10000000000000001 и 0.1000000000000000055511151231257827021181583404541015625 аппроксимируются 3602879701896397 / 2 ** 55. Поскольку все эти десятичные значения имеют одно и то же приближение, можно отобразить любое из них, сохранив при этом инвариант eval(repr(x)) == x.

Исторически сложилось так, что подсказка Python и встроенная функция repr() выбирали число с 17 значащими цифрами, 0.10000000000000001. Начиная с Python 3.1, Python (на большинстве систем) теперь может выбирать самое короткое из них и просто отображать 0.1.

Обратите внимание, что это заложено в самой природе двоичной плавающей точки: это не ошибка в Python, и это не ошибка в вашем коде тоже. Вы увидите то же самое во всех языках, поддерживающих аппаратную арифметику с плавающей точкой (хотя некоторые языки могут не отображать разницу по умолчанию или во всех режимах вывода).

Для более приятного вывода вы можете использовать форматирование строк, чтобы получить ограниченное количество значащих цифр:

>>> format(math.pi, '.12g')  # give 12 significant digits
'3.14159265359'

>>> format(math.pi, '.2f')   # give 2 digits after the point
'3.14'

>>> repr(math.pi)
'3.141592653589793'

Важно понимать, что в действительности это иллюзия: вы просто округляете отображение истинного значения машины.

Одна иллюзия может породить другую. Например, поскольку 0,1 не является в точности 1/10, суммирование трех значений 0,1 также не может дать в точности 0,3:

>>> 0.1 + 0.1 + 0.1 == 0.3
False

Кроме того, поскольку 0,1 не может быть ближе к точному значению 1/10, а 0,3 не может быть ближе к точному значению 3/10, то предварительное округление с помощью функции round() не поможет:

>>> round(0.1, 1) + round(0.1, 1) + round(0.1, 1) == round(0.3, 1)
False

Хотя числа нельзя приблизить к их предполагаемым точным значениям, функция math.isclose() может быть полезна для сравнения неточных значений:

>>> math.isclose(0.1 + 0.1 + 0.1, 0.3)
True

В качестве альтернативы можно использовать функцию round() для сравнения грубых приближений:

>>> round(math.pi, ndigits=2) == round(22 / 7, ndigits=2)
True

Двоичная арифметика с плавающей точкой таит в себе множество подобных сюрпризов. Проблема с «0.1» подробно описана ниже, в разделе «Ошибка представления». В разделе Examples of Floating Point Problems вы найдете краткое описание принципов работы двоичной арифметики с плавающей точкой и проблем, часто встречающихся на практике. Также смотрите The Perils of Floating Point для более полного описания других распространенных сюрпризов.

Как говорится в конце, «простых ответов нет». Тем не менее, не стоит излишне опасаться плавающей точки! Погрешности в операциях с плавающей запятой в Python наследуются от аппаратуры с плавающей запятой, и на большинстве машин составляют не более 1 части на 2**53 на операцию. Этого более чем достаточно для большинства задач, но нужно помнить, что это не десятичная арифметика и что каждая операция с плавающей запятой может содержать новую ошибку округления.

Хотя патологические случаи и существуют, для большинства повседневных применений арифметики с плавающей точкой вы получите ожидаемый результат, если просто округлите отображаемые результаты до ожидаемого количества десятичных цифр. Обычно достаточно str(), а для более тонкого контроля смотрите спецификаторы формата метода str.format() в Синтаксис строки форматирования.

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

Другая форма точной арифметики поддерживается модулем fractions, который реализует арифметику на основе рациональных чисел (поэтому числа типа 1/3 могут быть представлены точно).

Если вы активно используете операции с плавающей точкой, вам стоит обратить внимание на пакет NumPy и многие другие пакеты для математических и статистических операций, поставляемые проектом SciPy. См. <https://scipy.org>.

В Python есть инструменты, которые могут помочь в тех редких случаях, когда вам действительно нужно знать точное значение плавающей величины. Метод float.as_integer_ratio() выражает значение плавающей величины в виде дроби:

>>> x = 3.14159
>>> x.as_integer_ratio()
(3537115888337719, 1125899906842624)

Поскольку соотношение является точным, его можно использовать для воссоздания исходного значения без потерь:

>>> x == 3537115888337719 / 1125899906842624
True

Метод float.hex() выражает значение float в шестнадцатеричном виде (основание 16), снова давая точное значение, хранящееся в памяти вашего компьютера:

>>> x.hex()
'0x1.921f9f01b866ep+1'

Это точное шестнадцатеричное представление может быть использовано для точного восстановления плавающего значения:

>>> x == float.fromhex('0x1.921f9f01b866ep+1')
True

Поскольку представление является точным, оно полезно для надежного переноса значений между различными версиями Python (независимость от платформы) и обмена данными с другими языками, поддерживающими тот же формат (например, Java и C99).

Еще один полезный инструмент - функция sum(), которая помогает уменьшить потери точности при суммировании. Она использует расширенную точность для промежуточных шагов округления при добавлении значений к общему итогу. Это может повлиять на общую точность, так что ошибки не накапливаются до такой степени, что влияют на итоговый результат:

>>> 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 == 1.0
False
>>> sum([0.1] * 10) == 1.0
True

math.fsum() идет дальше и отслеживает все «потерянные цифры» при сложении значений с общим итогом, чтобы результат имел только одно округление. Это медленнее, чем sum(), но будет более точным в редких случаях, когда большие величины входов в основном аннулируют друг друга, оставляя конечную сумму около нуля:

>>> arr = [-0.10430216751806065, -266310978.67179024, 143401161448607.16,
...        -143401161400469.7, 266262841.31058735, -0.003244936839808227]
>>> float(sum(map(Fraction, arr)))   # Exact summation with single rounding
8.042173697819788e-13
>>> math.fsum(arr)                   # Single rounding
8.042173697819788e-13
>>> sum(arr)                         # Multiple roundings in extended precision
8.042178034628478e-13
>>> total = 0.0
>>> for x in arr:
...     total += x                   # Multiple roundings in standard precision
...
>>> total                            # Straight addition has no correct digits!
-0.0051575902860057365

15.1. Ошибка репрезентативности

В этом разделе подробно объясняется пример «0.1» и показывается, как можно самостоятельно провести точный анализ подобных случаев. Предполагается базовое знакомство с двоичным представлением с плавающей точкой.

Representation error относится к тому факту, что некоторые (большинство, на самом деле) десятичные дроби не могут быть представлены в точности как двоичные (основание 2). Это главная причина, по которой Python (или Perl, C, C++, Java, Fortran и многие другие) часто не отображает точное десятичное число, которое вы ожидаете.

Почему? 1/10 не совсем точно можно представить в виде двоичной дроби. По крайней мере с 2000 года почти все машины используют двоичную арифметику IEEE 754 с плавающей точкой, и почти все платформы отображают плавающие числа Python в значения IEEE 754 binary64 «двойной точности». Значения IEEE 754 binary64 содержат 53 бита точности, поэтому при вводе компьютер стремится преобразовать 0,1 в ближайшую возможную дробь вида J/2**N, где J - целое число, содержащее ровно 53 бита. Переписывание

1 / 10 ~= J / (2**N)

как

J ~= 2**N / 10

И если вспомнить, что J имеет ровно 53 бита (есть >= 2**52, но есть < 2**53), то наилучшее значение для N равно 56:

>>> 2**52 <=  2**56 // 10  < 2**53
True

То есть 56 - единственное значение для N, при котором в J остается ровно 53 бита. Наилучшее возможное значение для J - это округленный коэффициент:

>>> q, r = divmod(2**56, 10)
>>> r
6

Поскольку остаток больше половины от 10, наилучшее приближение получается при округлении в большую сторону:

>>> q+1
7205759403792794

Поэтому наилучшим приближением к 1/10 в IEEE 754 двойной точности будет:

7205759403792794 / 2 ** 56

Деление числителя и знаменателя на два сокращает дробь до:

3602879701896397 / 2 ** 55

Обратите внимание, что, поскольку мы округлили значение в большую сторону, оно на самом деле немного больше, чем 1/10; если бы мы не округлили, коэффициент был бы немного меньше, чем 1/10. Но ни в коем случае не может быть точно 1/10!

Поэтому компьютер никогда не «видит» 1/10: он видит точную дробь, приведенную выше, и лучшее двойное приближение IEEE 754, которое он может получить:

>>> 0.1 * 2 ** 55
3602879701896397.0

Если мы умножим эту дробь на 10**55, то получим значение с точностью до 55 десятичных цифр:

>>> 3602879701896397 * 10 ** 55 // 2 ** 55
1000000000000000055511151231257827021181583404541015625

означает, что точное число, хранящееся в компьютере, равно десятичному значению 0.1000000000000000055511151231257827021181583404541015625. Вместо того чтобы выводить полное десятичное значение, многие языки (включая старые версии Python) округляют результат до 17 значащих цифр:

>>> format(0.1, '.17f')
'0.10000000000000001'

Модули fractions и decimal упрощают эти вычисления:

>>> from decimal import Decimal
>>> from fractions import Fraction

>>> Fraction.from_float(0.1)
Fraction(3602879701896397, 36028797018963968)

>>> (0.1).as_integer_ratio()
(3602879701896397, 36028797018963968)

>>> Decimal.from_float(0.1)
Decimal('0.1000000000000000055511151231257827021181583404541015625')

>>> format(Decimal.from_float(0.1), '.17')
'0.10000000000000001'