Adaptive precision arithmetic

Эта статья находится в разработке!

Мотивация

Все вычисления, производимые компьютером во floating-point[1] модели, имеют погрешность. При большом количестве арифметических действий она возрастает. Во многих случаях результирующая погрешность уже не устраивает, и требуется либо абсолютно точное вычисление, либо меньшая погрешность. Одним из решений данной проблемы является хранение чисел в виде рациональных дробей, в которых числитель и знаменатель представляется в виде длинного целого числа. Но работать с такими числами довольно "дорого" по времени и тяжело в реализации: необходимо писать факторизацию чисел, эффективно сокращать дроби. Для улучшения работы нужны определенные оптимизации. Одной из них и является использование adaptive precision arithmetic.

Background

Большинство современных процессоров поддерживают числа с плавающей точкой в форме [math] \pm significand \times 2^{exponent}[/math]. Значащая часть числа (мантисса) представляет собой [math]p[/math]-битное двоичное число в форме [math]b.bbb \dots[/math], где каждое [math]b[/math] обозначает один бит. Также имеется один бит на знак.

Числа с плавающей точкой, как правило, нормализованы, то есть если число не равно нулю, то первый значимый бит равен единице, а экспонента устанавливается соответственно. Например, в [math]p[/math]-битной арифметике число 1101 (десятичное 13) будет выглядеть как [math]1.101 \times 2^3[/math].

Далее в этой статье фраза "что-либо представимо в [math]p[/math] битах" будет означать представимость в [math]p[/math] битах мантиссы, не учитывая знак и экспоненту.

Базовые понятия

Расширения

Определение:

Два числа [math]x[/math] и [math]y[/math] называются неперекрывающимися (англ. nonoverlapping), если номер наименьшего значимого бита числа [math]x[/math] (нумерация справа налево) больше, чем номер наибольшего значимого бита числа [math]y[/math], или наоборот.

Более формально, [math]x[/math] и [math]y[/math] не перекрываются, если существует такое целое число [math]r[/math] и [math]s[/math], что [math]x = r2^s[/math] и [math]|y| \lt 2^s[/math], или [math]y = r2^s[/math] и [math]|x| \lt 2^s[/math].

Ноль не пересекается ни с одним другим числом.

Например, числа 1100 и -10.1 не пересекаются, а 101 и 10 - пересекаются.

Иногда для использовании точной арифметики может понадобиться больше, чем [math]p[/math] бит для хранения величин. В связи с этим вводится одно из базовых форм хранения чисел для такой арифметики.

Определение:

Расширением числа [math]x[/math] называется такое его представление [math]x = x_n + x_{n-1} + \dots + x_1[/math], где каждое [math]x_i[/math] выражено [math]p[/math]-битным числом с плавающей точкой и называется компонентой этого расширения.

Определение:

Расширение называется неперекрывающимся, если все его компоненты взаимно не перекрываются.

Как правило, расширения должны быть неперекрывающимися, а их компоненты должны быть упорядочены от большей к меньшей по величине (то есть [math]x_n[/math] - большая). Далее будут рассматриваться именно такая их форма.

Стоит отметить, что число может быть представлено несколькими возможными неперекрывающимися расширениями: 1100 + -10.1 = 1001 + 0.1 = 1000 + 1 + 0.1.

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

Округление

Все алгоритмы, представленные в этой статье, предполагают, что сложение, вычитание и умножение производятся с точным округлением. Предполагается, что числа представляются в [math]p[/math] битах.

Определение:

Точное округление (англ. exact rounding) - такой вид округления, что: если точный результат может быть представлен в [math]p[/math] битах, то результатом округления будет точное значение числа; если точный результат не может быть представлен в [math]p[/math] битах, то результатом округления будет ближайшее [math]p[/math]-битное значение.

Например, в 4-битной арифметике произведение [math]111 \times 101 = 100011[/math] будет округлено в [math]1.001 \times 2^5[/math].

При вычислении результата может возникнуть ситуация, когда значение попадает в точности между двумя соседними [math]p[/math]-битными значениями. Тогда требуется определить правило поведения в таком случае. Рассмотрим некоторые из них.

Определение:

Округление до ближайшего четного (англ. round-to-even) - правило, при котором округление в вышеописанном случае производится к ближайшему [math]p[/math]-битному четному значению.

Определение:

Округление к нулю (англ. round-toward-zero) - правило, при котором округление в вышеописанном случае производится [math]p[/math]-битному значению, находящемуся между точным значением и нулем, а также ближайшему к точному значению.

Например, в 4-битной арифметике число [math]10011[/math] будет округлено до [math]1.010 \times 2^4[/math] по первому правилу, и до [math]1.001 \times 2^4[/math] по второму.

Стоит отметить, что стандарт IEEE 754 использует округление до ближайшего четного по умолчанию.

Далее в этой статье символами [math]\oplus, \ominus[/math] и [math]\otimes[/math] будут обозначаться [math]p[/math]-битные сложение, вычитание и умножение с точным округлением соответственно.

Из-за округления данные операции теряют некоторые важные свойства, например, ассициативность: [math](1000 \oplus 0.011) \oplus 0.011 = 1000[/math], но [math]1000 \oplus(0.011 \oplus 0.011) = 1001[/math].

При анализе округления часто используют так называемый ulp.

Определение:

ULP (англ. units in the last place) - эффективная величина самого младшего ([math]p[/math]-ого) бита.

Например, [math]ulp(-1100) = 1, ulp(1) = 0.001[/math] в [math]p[/math]-битной арифметике.

Так же полезной нотацией является [math] err(a \circledast b) [/math], которая обозначает ошибку округлении результата выполнения операции [math]\circledast[/math]. Отметим, что если ulp всегда беззнаковая величина, то err может иметь знак. Для базовых операций (сложение, вычитание, умножение) [math] a \circledast b = a \ast b + err(a \circledast b)[/math], и точное округление гарантирует, что [math] |err(a \circledast b)| \leqslant \frac{1}{2}ulp(a \circledast b)[/math].

Свойства

Иногда есть возможность найти более точные границы ошибки округления, что будет видно далее из лемм. Первая лемма используется, когда один операнд много меньше другого, а вторая - когда сумма близка к степени двойки. Для лемм 1 - 4 пусть [math]a, b[/math] - [math]p[/math]-битные числа с плавающей точкой. Леммы приводятся без доказательств, их можно найти в статье Джонатана Шевчука "Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates".

Лемма:

Пусть [math]a \oplus b = a + b + err(a \oplus b) [/math]. Ошибка округления [math]err(a \oplus b)[/math] не превзойдет [math]max(|a|, |b|)[/math]. (Аналогично для вычитания).

TODO: add picture

Как следствие, верна следующая лемма:

Лемма:

Ошибка округления [math]err(a \oplus b)[/math] может быть представлена в [math]p[/math] битах.

Лемма:

Пусть [math]|a + b| \leqslant |b|[/math] и [math]|a + b| \leqslant |a|[/math]. Тогда [math]a \oplus b = a + b[/math]. (Аналогично для вычитания).

Лемма:

Пусть [math]b \in \left [ \frac{a}{2}, 2a \right ][/math]. Тогда [math]a \ominus b = a - b[/math].

Алгоритмы суммирования

Простое суммирование

Важной базовой операцией во всех алгоритмах, основанных на представлении чисел в виде расширений, является сумма двух [math]p[/math]-битных величин, результатом которой является расширение длины два. Есть два алгоритма для выполнения этой задачи - алгоритмы Деккера и Кнута.

Теорема (Dekker):

Пусть [math]a[/math] и [math]b[/math] есть [math]p[/math]-битные числа, причем [math]|a| \geqslant |b|[/math]. Тогда следующий алгоритм вернет непересекающееся расширение [math]x + y[/math] такое, что [math]a + b = x + y[/math], где [math]x[/math] - приближение (аппроксимация) суммы [math]a + b[/math], а [math]y[/math] представляет собой ошибку округления при вычислении [math]x[/math]. (Аналогично для вычитания).

TODO: add pictures

Псевдокод для суммы:

[math]FastTwoSum(a, b)[/math]
[math]1[/math]   [math]x \Leftarrow a \oplus b[/math]
[math]2[/math]   [math]b_{virtual} \Leftarrow x \ominus a[/math]
[math]3[/math]   [math]y \Leftarrow b \ominus b_{virtual}[/math]
[math]4[/math]   [math]return (x, y)[/math]

Псевдокод для суммы:

[math]FastTwoSum(a, b)[/math]
[math]1[/math]   [math]x \Leftarrow a \ominus b[/math]
[math]2[/math]   [math]b_{virtual} \Leftarrow x \ominus a[/math]
[math]3[/math]   [math]y \Leftarrow b \ominus b_{virtual}[/math]
[math]4[/math]   [math]return (x, y)[/math]

Отметим, что величины на выходе этой процедуры вовсе не должны иметь один знак.