Aрифметика произвольной точности в Erlang

@rawpixel

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

Если говорить про Erlang, то он, как и многие другие языки, реализует IEEE754 стандарт для float, в то время как стандартный тип Integer в Erlang реализован с использованием арифметики произвольной точности. Однако, хотелось бы иметь не только bigint, но и возможность оперирования рациональными, комплексными и числами с плавающей точкой с необходимой точностью.

В статье представлен минимальный обзор теории кодирования чисел с плавающей точкой и наиболее яркие примеры возникающих эффектов. Решение, обеспечивающее необходимую точность операций через переход в представление с фиксированной точкой, оформлено в виде библиотеки EAPA (Erlang Arbitrary Precision Arithmetic), призванной удовлетворить потребности финансовых приложений, разрабатываемых на Erlang / Elixir.

---

Стандарты, стандарты, стандарты...

На сегодняшний день основным стандартом двоичной арифметики с плавающей точкой является широко применяемый в технике и программировании IEEE754. Он определяет четыре формата представления:

• с одинарной точностью (single-precision) 32 бита
• с двойной точностью (double-precision) 64 бита
• с одинарной расширенной точностью (single-extended precision) >=43 бит (редко используемый)
• с двойной расширенной точностью (double-extended precision) >= 79 бит (обычно используют 80 бит)
  и четыре режима округления:
• Округление, стремящееся к ближайшему целому.
• Округление, стремящееся к нулю.
• Округление, стремящееся к +∞
• Округление, стремящееся к -∞

Большинство современных микропроцессоров изготовляются с аппаратной реализацией представления вещественных переменных в формате IEEE754. Форматы представления ограничивают предельный размер числа, а режимы округления влияют на точность. Программисты зачастую не могут изменить поведение аппаратуры и реализации языков программирования. Например, официальная реализация Erlang хранит float в 3 словах на 64 разрядной машине и в 4 словах на 32 разрядной.

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

Основная масса чисел при отображении на конечное множество имеет стабильную и небольшую относительную погрешность. Так, для float она составляет 11,920928955078125e-6 %, а для double — 2,2204460492503130808472633361816e-14 %. В жизни программистов большинство решаемых повседневных задач позволяют пренебречь этой погрешностью, хотя следует учесть, что и в простых задачах можно наступить на грабли, так как величина абсолютной ошибки может достигать 10³¹ и 10²⁹² для float и double соответственно, вызывая затруднения в вычислениях.

Иллюстрация эффектов

От общей информации к делу. Давайте попробуем воспроизвести возникающие эффекты в Erlang.

Все примеры, приведенные ниже, оформлены в виде ct-тестов.

Округление и потеря точности

Начнем с классики – сложения двух чисел: 0,1+0,2 = ?:

t30000000000000004(_)->
 ["0.30000000000000004"] = io_lib:format("~w", [0.1 + 0.2]).

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

t30000000000000004_eapa(_)->
 %% prec = 1 symbols after coma
 X = eapa_int:with_val(1, <<"0.1">>),
 Y = eapa_int:with_val(1, <<"0.2">>),
 <<"0.3">> = eapa_int:to_float(1, eapa_int:add(X, Y)).

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

tiny(_)->
 X = 1.0,
 Y = 0.0000000000000000000000001,
 1.0 = X + Y.

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

tiny_eapa(_)->
 X1 = eapa_int:with_val(1, <<"1.0">>),
 X2 = eapa_int:with_val(25, <<"0.0000000000000000000000001">>),
 <<"1.0000000000000000000000001">> = eapa_int:to_float(eapa_int:add(X1, X2)).

Переполнение разрядной сетки

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

float_overflow(_) ->
 1.0 = 9007199254740991.0 - 9007199254740990.0,
 1.0 = 9007199254740992.0 - 9007199254740991.0,
 0.0 = 9007199254740993.0 - 9007199254740992.0,
 2.0 = 9007199254740994.0 - 9007199254740993.0.

Как видно из теста, в какой-то момент разница перестает быть равна 1.0, что очевидно является проблемой. EAPA позволяет решить и эту проблему:

float_overflow_eapa(_)->
 X11 = eapa_int:with_val(1, <<"9007199254740992.0">>),
 X21 = eapa_int:with_val(1, <<"9007199254740991.0">>),
 <<"1.0">> = eapa_int:to_float(1, eapa_int:sub(X11, X21)),
 X12 = eapa_int:with_val(1, <<"9007199254740993.0">>),
 X22 = eapa_int:with_val(1, <<"9007199254740992.0">>),
 <<"1.0">> = eapa_int:to_float(1, eapa_int:sub(X12, X22)),
 X13 = eapa_int:with_val(1, <<"9007199254740994.0">>),
 X23 = eapa_int:with_val(1, <<"9007199254740993.0">>),
 <<"1.0">> = eapa_int:to_float(1, eapa_int:sub(X13, X23)). 

Опасная редукция

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

Покажем, что в Erlang эта проблема присутствует:

reduction(_)->
 X = float(87654321098765432),
 Y = float(87654321098765431),
 16.0 = X-Y. %% has to be 1.0

Получилось 16.0 вместо ожидаемого 1.0. Попробуем исправить данную ситуацию:

reduction_eapa(_)->
 X = eapa_int:with_val(1, <<"87654321098765432">>),
 Y = eapa_int:with_val(1, <<"87654321098765431">>),
 <<"1.0">> = eapa_int:to_float(eapa_int:sub(X, Y)).

Другие особенности арифметики с плавающей точкой в Erlang

Начнем, пожалуй, с игнорирования отрицательного нуля.

eq(_)->
 true = list_to_float("0.0") =:= list_to_float("-0.0").

Сразу хочется сказать, что EAPA сохраняет данное поведение:

eq_eapa(_)->
 X = eapa_int:with_val(1, <<"0.0">>),
 Y = eapa_int:with_val(1, <<"-0.0">>),
 true = eapa_int:eq(X, Y).

так как оно допустимо. В Erlang нет внятного синтаксиса и обработки NaN и бесконечностей, что порождает целый ряд особенностей, например, вот таких:

1> math:sqrt(list_to_float("-0.0")).
0.0

Следующим моментом является особенность обработки больших и маленьких чисел. Попробуем воспроизвести для маленьких:

2> list_to_float("0."++lists:duplicate(322, $0)++"1").
1.0e-323
3> list_to_float("0."++lists:duplicate(323, $0)++"1").
0.0

и для больших чисел:

4> list_to_float("1"++lists:duplicate(308, $0)++".0").
1.0e308
5> list_to_float("1"++lists:duplicate(309, $0)++".0").
** exception error: bad argument

Приведем еще пару примеров для маленьких чисел:

6> list_to_float("0."++lists:duplicate(322, $0)++"123456789").
1.0e-323
7> list_to_float("0."++lists:duplicate(300, $0)++"123456789").
1.23456789e-301

8> 0.123456789e-100 * 0.123456789e-100.
1.524157875019052e-202
9> 0.123456789e-200 * 0.123456789e-200.
0.0

Приведенные выше примеры подтверждают истину и для Erlang проектов: деньги нельзя считать в IEEE754.

EAPA (Erlang Arbitrary-Precision Arithmetic)

EAPA – NIF расширение, написанное на Rust. На данный момент в репозитории EAPA представлен максимально простой и удобный интерфейс eapa_int для работы с числами с фиксированной точкой. К числу особенностей eapa_int можно отнести следующие:

1. Отсутствие эффектов IEEE754 кодирования
2. Поддержка больших чисел
3. Настраиваемая точность до 126 знаков после запятой. (в текущей реализации)
4. Автоскейлинг
5. Поддержка всех основных операций над числами
6. Более или менее полное тестирование, в том числе property based.

Интерфейс eapa_int:

• with_val/2 – перевод числа с плавающей точкой в фиксированное представление, которое можно в том числе безопасно использовать в json, xml.
• to_float/2 – перевод числа с фиксированной точкой в число с плавающей точкой с заданной точностью.
• to_float/1 – перевод числа с фиксированной точкой в число с плавающей точкой.
• add/2 – сумма двух чисел
• sub/2 – разность
• mul/2 – умножение
• divp/2 – деление
• min/2 – минимальное из чисел
• max/2 – максимальное из чисел
• eq/2 – проверка равенства чисел
• lt/2 – проверка, что число меньше
• lte/2 – проверка меньше равно
• gt/2 – проверка, что число больше
• gte/2 – проверка больше равно

Код EAPA можно найти в репозитории https://github.com/Vonmo/eapa

Когда следует использовать eapa_int? Например, если ваше приложение работает с деньгами или вам нужно удобно и точно производить вычислительные операции над числами вида 92233720368547758079223372036854775807.92233720368547758079223372036854775807, можно смело использовать EAPA.

Как и любое решение, EAPA это компромисс. Мы получаем необходимую точность, жертвуя памятью и скоростью вычислений, Тесты производительности и статистика, собранная на реальных системах, показывают, что большинство операций выполняются в диапазоне 3-30 мкс. Данный момент тоже необходимо учитывать при выборе интерфейса с фиксированной точкой EAPA.

Заключение

Конечно, далеко не всегда приходится решать подобные задачи на Erlang или Elixir, но когда возникает проблема и не находится подходящий инструмент, приходится изобретать решение.
Данная статья – это попытка поделиться с сообществом инструментом и опытом, в надежде, что для кого-то эта библиотека окажется полезной и поможет сохранить время.
А как вы считаете деньги в Erlang?

P.S. Работа с рациональными и комплексными числами, а также нативный доступ к Integer, Float, Complex, Rational типам произвольной точности будут освещены в следующих публикациях. Не переключайтесь!

---

Материалы по теме:

• What Every Computer Scientist Should Know About Floating-Point Arithmetic
• Что нужно знать про арифметику с плавающей запятой
• IEEE Standard 754 Floating Point Numbers
• IEEE 754 — стандарт двоичной арифметики
  с плавающей точкой
• 0.30000000000000004

Автор: Максим Молчанов

Источник