Проблема, о которой вы наверняка не задумывались: print(.1+.2) / Хабр

Как следует отображать на экране результат деления 3.0 на 10.0 ? Сколько цифр следует вывести, если пользователь не указал точность?

Скорее всего, вы даже не знали, что вывод на экран чисел с плавающей запятой — это сложная проблема, настолько сложная, что по ней написаны десятки научных статей, причём последний прорыв был относительно недавно, в 2016 году. На самом деле, это одна из самых сложных частей поддержки чисел с плавающей запятой в среде выполнения языка.

Давайте продолжим разговор о самой неоптимизированной в мире библиотеке эмуляции плавающей точки при помощи целочисленной арифметики.

Это вторая статья из цикла «Санпросвет о плавающей точке»:

Компьютеры и числа
Вывод чисел с плавающей точкой на экран <- вы тут

Три основные проблемы

Спойлер: в моей реализации я намеренно обойду самые сложные части,
но давайте посмотрим, что делает эту задачу такой сложной.

Проблема № 1: нужное количество цифр

Обратите внимание, что я написал числа 3.0 и 10.0 в десятичной системе счисления (база 10), потому что это система, которую используют люди. Но большинство десятичных дробей не имеют точного двоичного представления:

$0.3_{10} = 0.0100110011001100110011001100\dots_2.$

Эта двоичная дробь бесконечна, паттерн 0011 повторяется бесконечно.
Поэтому, когда мы пишем x = 0.3, переменная x на самом деле хранит ближайшее представимое число с плавающей запятой.

Например, в Python плавающая запятая имеет двойную точность (64 бита). Таким образом, x на самом деле хранит 0.299999999999999988897769753748434595763683319091796875.

Вы можете проверить это сами:

from decimal import Decimal
 x = 0.3
 print(Decimal(x))

Если мы создадим вторую переменную y = 0.1 + 0.2, то сохраненное значение будет: 0.3000000000000000444089209850062616169452667236328125.

from decimal import Decimal
y = 0.1 + 0.2
print(Decimal(y))

Неудивительно, что x и y различаются: y накопило три ошибки аппроксимации — по одной для каждого операнда и еще одну от сложения:

0.1 округляется до 0.1000000000000000055511151231257827021181583404541015625
0.2 округляется до 0.200000000000000011102230246251565404236316680908203125
Их сумма вносит еще одну ошибку округления.

Теперь вопрос: как мы должны вывести x и y? Давайте посмотрим, как это делает Python:

Первое требование — безопасность при обратном преобразовании:
какую бы десятичную строку мы ни вывели, при обратном преобразовании должно быть восстановлено точно то же двоичное число с плавающей запятой.

Если система выводит слишком много цифр, лишние цифры могут быть просто «мусором» от двоичного округления.
Даже если x равно 0.299999999999999988897769753748434595763683319091796875, нет необходимости печатать все 54 цифры после десятичной запятой. На самом деле мы хотели сохранить 0.3, поэтому 54 цифры являются фактически численным мусором. Простого 0.3 достаточно для определения значения x, даже если x не равно .
Если система выводит слишком мало цифр, результат будет неверным в более строгом смысле: обратное преобразование в двоичную систему может дать другое число с плавающей запятой.
Например, самым коротким десятичным числом, которое однозначно идентифицирует y, является 0.30000000000000004. Вывод только 0.3 будет неправильным, потому что 0.3 при обратном преобразовании дает другое двоичное значение.

Таким образом, правило таково:
Нужно вывести самую короткую десятичную строку, которая точно преобразуется в обратном направлении.

Проблема № 2: скорость

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

Именно поэтому существуют специализированные алгоритмы:

Dragon4 (классический, точный, но сложный)
Grisu3, Errol3, Ryu (более новые, более быстрые, доказательно минимальные)

Эти алгоритмы являются одними из самых сложных частей языковых сред выполнения (C, Java, Python, JavaScript и т. д.). Кроме того, производительность тесно связана со следующей сложностью.

Проблема № 3: только правильные цифры

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

$\frac{\text{целочисленный числитель}}{\text{степень двойки}}$

Но число с плавающей запятой имеет огромный динамический диапазон. Наибольшее конечное число типа double превышает $10^{308}$ . Это число имеет 308 десятичных цифр, что намного превышает возможности даже 128-битного целого числа (максимум $10^{38}$ ). Поэтому нам нужна высокоточная арифметика (большие числа), чтобы выполнить точное деление без переполнения.

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

Печать точного значения

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

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

Более десяти лет назад Брюс Доусон опубликовал реализацию на C++, которая выводит точное сохраненное значение float. Его метод использовал только арифметику int32_t для одинарной точности, но полагался на частичную реализацию длинной арифметики для десятичного преобразования (путем повторяющегося деления на 10).

Технически это прекрасно, но я хочу что-то более простое и прозрачное.

Класс Float7

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

n_e = 3
n_m = 7-n_e

anchors = [ 0 ]
for e in range(-2**(n_e-1)+1, 2**(n_e-1)+1): # prepare 2**n_e intervals
    anchors.append(2**e)

numbers = []
for i in range(len(anchors)-1): # for each interval
    for m in range(2**n_m):     # populate with 2**n_m numbers
        v = anchors[i] + m/2**n_m * (anchors[i+1]-anchors[i])
        numbers.append(v)
print(numbers)

Это делит диапазон на 8 интервалов, каждый из которых заполнен 16 равномерно распределенными числами. Всего у нас есть 128 значений:

[0.0, 0.0078125, 0.015625, 0.0234375, 0.03125, 0.0390625, 0.046875, 0.0546875, 0.0625, 0.0703125, 0.078125, 0.0859375, 0.09375, 0.1015625, 0.109375, 0.1171875, 0.125, 0.1328125, 0.140625, 0.1484375, 0.15625, 0.1640625, 0.171875, 0.1796875, 0.1875, 0.1953125, 0.203125, 0.2109375, 0.21875, 0.2265625, 0.234375, 0.2421875, 0.25, 0.265625, 0.28125, 0.296875, 0.3125, 0.328125, 0.34375, 0.359375, 0.375, 0.390625, 0.40625, 0.421875, 0.4375, 0.453125, 0.46875, 0.484375, 0.5, 0.53125, 0.5625, 0.59375, 0.625, 0.65625, 0.6875, 0.71875, 0.75, 0.78125, 0.8125, 0.84375, 0.875, 0.90625, 0.9375, 0.96875, 1.0, 1.0625, 1.125, 1.1875, 1.25, 1.3125, 1.375, 1.4375, 1.5, 1.5625, 1.625, 1.6875, 1.75, 1.8125, 1.875, 1.9375, 2.0, 2.125, 2.25, 2.375, 2.5, 2.625, 2.75, 2.875, 3.0, 3.125, 3.25, 3.375, 3.5, 3.625, 3.75, 3.875, 4.0, 4.25, 4.5, 4.75, 5.0, 5.25, 5.5, 5.75, 6.0, 6.25, 6.5, 6.75, 7.0, 7.25, 7.5, 7.75, 8.0, 8.5, 9.0, 9.5, 10.0, 10.5, 11.0, 11.5, 12.0, 12.5, 13.0, 13.5, 14.0, 14.5, 15.0, 15.5]

У нас есть 128 чисел $v_{e,m}$ , индексированных 8 значениями $e\in[-4\dots 3]$ и 16 значениями $m\in[0\dots 15]$ . Вставим явно значения опорных точек в строку 11 вышеприведенного фрагмента кода на Python.

Таким образом, мы можем записать общую формулу:

$\scriptstyle v_{e,m} = \left\{\begin{array}{lll}0 + \frac{m}{2^{n_m}}\left(2^{e+1} - 0\right) & =\quad m \cdot 2^{e+1-n_m} \quad & \text{если}~e= -2^{n_e-1}\\ 2^e + \frac{m}{2^{n_m}}\left(2^{e+1} - 2^e\right) & = \quad (m+2^{n_m}) \cdot 2^{e-n_m} & \text{в противном случае}\end{array}\right.$

Это отражает различие между денормализованными числами (первый случай) и нормальными (второй случай).

Следовательно, мы можем переписать вышеуказанный фрагмент в следующей форме:

numbers = [] 
for e in range(-2**(n_e-1), 2**(n_e-1)):
    for m in range(2**n_m):
        if e == -2**(n_e-1):
            e += 1          # subnormal number
        else:
            m += 2**n_m     # normal number
        v = m * 2**(e-n_m)
        numbers.append(v)
print(numbers)

Обратите внимание на строки 5 и 7:

мы увеличиваем $e\leftarrow e+1$ , если число является денормализованным,
и добавляем $2^{n_m}$ к , если число является нормальным $m\leftarrow m+2^{n_m}$ .

Таким образом, мы получаем уникальную формулу для конечного числа, которое мы генерируем в строке 9. Фактически это означает, что, несмотря на то, что мы фактически храним битов числа в памяти, это -битное число без знака. Дополнительный 5-й бит нашего Float7 не хранится явно, он «скрыт» и может быть восстановлен из значения e.

Это также означает, что мантисса кодирует число с фиксированной запятой с 1 битом для целой части и $2^{n_m}$ битами для дробной части.
Следовательно, $\frac m {2^{n_m}} \in [0, 2)$ . Умножив его на , мы сдвигаем точку на двоичных разрядов.

Вместо перечисления всех возможных чисел с плавающей запятой, давайте создадим класс Float7 для представления одного значения:

class Float7:
    def __init__(self, uint7):
        self.e = uint7 // 16 - 4
        self.m = uint7 %  16
        if self.e == -4:
            self.e += 1  # if subnormal, increment the exponent, the hidden bit = 0
        else:
            self.m += 16 # if normal, recover the hidden bit = 1

    def __float__(self):
        return self.m * 2**(self.e-4)

print(float(Float7(93)))

Конструктор (строка 2) принимает 7-битный шаблон и разделяет его на экспоненту и мантиссу , восстанавливая скрытый бит. Оператор __float__() преобразует наш пользовательский класс в нативные числа с плавающей запятой Python. Наконец, код выводит 93е число с плавающей запятой:

3.625

На всякий случай проверим вычисления:

$93_{10} = 101\ 1101_2$

Конструктор получает битовые шаблоны 101 и 1101 для экспоненты и мантиссы. Знаковая интерпретация со смещением говорит нам, что 101 интерпретируется как $1_{10}$ :

$e = \left\lfloor\frac{93}{16}\right\rfloor-4 = 1$

Число является нормальным, поэтому мы должны восстановить скрытый бит 1 для мантиссы:

$m = (93 \mod 16) + 16 = 29$

Наконец, приведение типа __float__() дает нам 3.625, как и должно быть:

$29 \cdot 2^{1-4} = 3.625$

Сложение в столбик

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

Оказывается, что наш Float7 может быть представлен в виде числа с фиксированной запятой, имеющего 4 бита в целой части и 7 бит в дробной части, что можно обозначить как формат Q4.7. Мы можем определить это, учитывая, что мантисса числа с плавающей запятой является числом с фиксированной запятой Q1.4. Максимальная экспонента числа с плавающей запятой равна 3, что эквивалентно сдвигу мантиссы влево на 3 позиции. Минимальная экспонента числа с плавающей запятой равна -3, что эквивалентно сдвигу мантиссы вправо на 3 позиции. Эти величины сдвига нашей мантиссы Q1.4 означают, что все числа с плавающей запятой могут поместиться в число с фиксированной запятой Q4.7, что составляет в общей сложности 11 бит.

Скрытый текст

Мы можем представить все значения float32 без знака, включая денормализованные, с помощью 277-битного формата с фиксированной запятой Q128.149. 277 бит — это огромное число, и чаще всего для его обработки используется длинная арифметика. Я хочу полностью избежать арифметических операций с такими большими числами.

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

Рассмотрим пример из предыдущего фрагмента кода на Python.

$m = 29_{10} = 11101_2$

Мантисса хранит число с фиксированной запятой Q1.4 . Экспонента говорит нам, что нам нужно сдвинуть точку основания вправо на 1 разряд, поэтому представленное число равно . Мы можем скопировать и вставить этот шаблон в число с фиксированной запятой Q4.7 0011.1010000.

Назовем битовый шаблон $b = \{0,0,0,0,1,0,1,1,1,0,0\}$ , где обозначает i-й бит. Значение можно вычислить как сумму битового шаблона, умноженного на соответствующие степени двойки:

$v = \sum_{i=0}^{10} b_i\ 2^{i-7}$

Мы можем восстановить его десятичное представление, выполнив то, что мы так любили в начальной школе, — сложение столбцами:

     0,0078125  * 0
   + 0,0156250  * 0
   + 0,0312500  * 0
   + 0,0625000  * 0
   + 0,1250000  * 1
   + 0,2500000  * 0
   + 0,5000000  * 1
   + 1,0000000  * 1
   + 2,0000000  * 1
   + 4,0000000  * 0
     8,0000000  * 0
   = ---------
     3,6250000

Да, действительно равно $3,625_{10}$ !

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

Вот полный код:

class Float7:
    def __init__(self, uint7):
        self.e = uint7 // 16 - 4
        self.m = uint7 %  16
        if self.e == -4:
            self.e += 1  # if subnormal, increment the exponent, the hidden bit = 0
        else:
            self.m += 16 # if normal, recover the hidden bit = 1

    def __float__(self):
        return self.m * 2**(self.e-4)

    def __str__(self):
        return str(Q4_7(self.e + 3, self.m))

class Q4_7:  # 11-bit number with 4 bits in the integer part and 7 in the fraction
    digits = [
        [5,0,0,0,0,0,0,0,0,0,0], # constant array, 11 powers of 2 in base 10
        [2,5,0,0,0,0,0,0,0,0,0],
        [1,2,5,0,0,0,0,0,0,0,0],
        [8,6,2,5,0,0,0,0,0,0,0],
        [7,5,1,2,5,0,0,0,0,0,0],
        [0,1,3,6,2,5,0,0,0,0,0],
        [0,0,0,0,1,2,5,0,0,0,0], # the decimal dot is here
        [0,0,0,0,0,0,0,1,2,4,8],
        [0,0,0,0,0,0,0,0,0,0,0]  # zero padding to avoid extra logic in line 40
    ]

    def __init__(self, offset, uint5):
        self.number = [ 0 ]*11   # 11-bit number, binary expansion of uint5 * 2**offset
        while uint5 > 0:
            self.number[offset] = uint5 % 2
            uint5 //= 2
            offset += 1

    def __str__(self):
        string = ''
        carry = 0
        for position in range(9):             # loop from the least significant digit to the most significant
            total = sum(bit * dgt for bit, dgt in zip(self.number, self.digits[position])) # sum of 11 digits
            digit = str((total + carry) % 10) # current digit to output
            carry = (total + carry)//10
            string = digit + string
            if position==6:                   # decimal dot position
                string = '.' + string
        return string
print(Float7(93))

В строке 14 битовый шаблон вставляется в нужное место 11-разрядного числа с фиксированной запятой Q4_7. Затем в строках 36-46 функция преобразования строк __str__() просто суммирует столбцы над 11 степенями двойки, хранящимися в массиве digits. Все очень просто!

Вот результат:

03.6250000

Удаляем ведущие и конечные нули, и все готово.

Вывод

Мы увидели, что вывод чисел с плавающей запятой — это удивительно сложная проблема. В то время как все более-менее нормальные языки используют сложные алгоритмы (Ryu, Grisu, Dragon4), для отладки реализаций soft-float мы можем пойти на хитрость: представлять значения в виде фиксированной запятой, заранее вычислить десятичные разложения и использовать простое сложение в столбик.

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

Проблема, о которой вы наверняка не задумывались: print(.1+.2)