Adaptive precision arithmetic — различия между версиями

Материал из Викиконспекты
Перейти к: навигация, поиск
м (rollbackEdits.php mass rollback)
 
(не показана 41 промежуточная версия 4 участников)
Строка 1: Строка 1:
 
==Мотивация==
 
==Мотивация==
Все вычисления, производимые компьютером во <tex>floating-point</tex>[http://en.wikipedia.org/wiki/Floating_point]'' модели, имеют погрешность. При большом количестве арифметических действий она возрастает. Во многих случаях результирующая погрешность уже не устраивает, и требуется либо абсолютно точное вычисление, либо меньшая погрешность. Одним из решений данной проблемы является хранение чисел в виде рациональных дробей, в которых числитель и знаменатель представляется в виде длинного целого числа. Но работать с такими числами довольно "дорого" по времени и тяжело в реализации: необходимо писать факторизацию чисел, эффективно сокращать дроби. Для улучшения работы нужны определенные оптимизации. Одной из них и является использование adaptive precision arithmetic.
+
Все вычисления, производимые компьютером в модели с ''плавающей точкой'', имеют погрешность. При большом количестве арифметических действий она возрастает. Во многих случаях результирующая погрешность уже не устраивает, и требуется либо абсолютно точное вычисление, либо меньшая погрешность. Одним из решений данной проблемы является хранение чисел в виде рациональных дробей, в которых числитель и знаменатель представляется в виде длинного целого числа. Но работать с такими числами довольно "дорого" по времени и тяжело в реализации: необходимо писать факторизацию чисел, эффективно сокращать дроби. Для улучшения работы нужны определенные оптимизации. Одной из них и является использование ''адаптивной арифметики'' (англ. ''adaptive precision arithmetic'').
  
==Background==
+
==Базовые понятия==
Большинство современных процессоров поддерживают числа с плавающей точкой в форме <tex> \pm significand \times 2^{exponent}</tex>. Значащая часть числа (мантисса) представляет собой <tex>p</tex>-битное двоичное число в форме <tex>b.bbb \dots</tex>, где каждое <tex>b</tex>   
+
===Представление чисел===
 +
Большинство современных процессоров поддерживают числа с плавающей точкой в форме <tex> \pm \mathrm{significand} \times 2^{\mathrm{exponent}}</tex>. Значащая часть числа (мантисса) представляет собой <tex>p</tex>-битное двоичное число в форме <tex>b.bbb \dots</tex>, где каждое <tex>b</tex>   
 
обозначает один бит. Также имеется один бит на знак.
 
обозначает один бит. Также имеется один бит на знак.
  
Строка 10: Строка 11:
 
Далее в этой статье фраза "что-либо представимо в <tex>p</tex> битах" будет означать представимость в <tex>p</tex> битах мантиссы, ''не'' учитывая  знак и экспоненту.
 
Далее в этой статье фраза "что-либо представимо в <tex>p</tex> битах" будет означать представимость в <tex>p</tex> битах мантиссы, ''не'' учитывая  знак и экспоненту.
  
==Базовые понятия==
+
===Разложения===
===Расширения===
 
 
{{Определение
 
{{Определение
 
|definition=
 
|definition=
Строка 21: Строка 21:
 
Ноль не пересекается ни с одним другим числом.
 
Ноль не пересекается ни с одним другим числом.
  
'''Например''', числа <tex>1100</tex> и <tex>-10.1</tex> не пересекаются, а <tex>110</tex> и <tex>10</tex> - пересекаются.
+
''Например'', числа <tex>1100</tex> и <tex>-10.1</tex> не пересекаются, а <tex>110</tex> и <tex>10</tex> {{---}}  пересекаются.
  
 
{{Определение
 
{{Определение
Строка 28: Строка 28:
 
}}
 
}}
  
'''Например''', числа <tex>1100</tex> и <tex>11</tex> - смежные, а <tex>1100</tex> и <tex>1000</tex> - нет.
+
''Например'', числа <tex>1100</tex> и <tex>11</tex> {{---}}  смежные, а <tex>1100</tex> и <tex>1000</tex> {{---}}  нет.
  
 
Иногда для использовании точной арифметики может понадобиться больше, чем <tex>p</tex> бит для хранения величин. В связи с этим вводится одно из базовых форм хранения чисел для такой арифметики.
 
Иногда для использовании точной арифметики может понадобиться больше, чем <tex>p</tex> бит для хранения величин. В связи с этим вводится одно из базовых форм хранения чисел для такой арифметики.
 
{{Определение
 
{{Определение
 
|definition=
 
|definition=
'''Расширением''' числа <tex>x</tex> называется такое его представление <tex>x = x_n + x_{n-1} + \dots + x_1</tex>, где каждое <tex>x_i</tex> выражено <tex>p</tex>-битным числом с плавающей точкой и называется '''компонентой''' этого расширения.
+
'''Разложением''' числа <tex>x</tex> называется такое его представление <tex>x = x_n + x_{n-1} + \dots + x_1</tex>, где каждое <tex>x_i</tex> выражено <tex>p</tex>-битным числом с плавающей точкой и называется '''компонентой''' этого разложения.
 
}}
 
}}
  
 
{{Определение
 
{{Определение
 
|definition=
 
|definition=
Расширение называется '''неперекрывающимся (несмежным)''', если все его компоненты взаимно не перекрываются (не являются смежными).
+
Разложение называется '''неперекрывающимся (несмежным)''', если все его компоненты взаимно не перекрываются (не являются смежными).
 
}}
 
}}
  
Как правило, расширения должны быть неперекрывающимися, а их компоненты должны быть упорядочены от большей к меньшей по величине (то есть <tex>x_n</tex> - большая). Далее будут рассматриваться именно такая их форма.
+
Как правило, разложения должны быть неперекрывающимися, а их компоненты должны быть упорядочены от большей к меньшей по величине (то есть <tex>x_n</tex> {{---}} большая). Далее будут рассматриваться именно такая их форма.
  
<wikitex>
+
<wikitex>Стоит отметить, что число может быть представлено несколькими возможными неперекрывающимися разложениями: $1100 + -10.1 = 1001 + 0.1 = 1000 + 1 + 0.1$.
Стоит отметить, что число может быть представлено несколькими возможными неперекрывающимися расширениями: $1100 + -10.1 = 1001 + 0.1 = 1000 + 1 + 0.1$.
 
 
</wikitex>
 
</wikitex>
Неперекрывающиеся расширения нужны, например, для того, чтобы быстро вычислять знак выражения (смотрим знак большей по размеру компоненты), или для грубой оценки значения всего расширения (берем большую по величине компоненту).
+
Неперекрывающиеся разложения нужны, например, для того, чтобы быстро вычислять знак выражения (смотрим знак большей по размеру компоненты), или для грубой оценки значения всего разложения (берем большую по величине компоненту).
  
 
===Округление===
 
===Округление===
Строка 52: Строка 51:
 
{{Определение
 
{{Определение
 
|definition=
 
|definition=
'''Точное округление''' (англ. ''exact rounding'') - такой вид округления, что:
+
'''Точное округление''' (англ. ''exact rounding'') {{---}}  такой вид округления, что:
 
* если точный результат ''может'' быть представлен в <tex>p</tex> битах, то результатом округления будет точное значение числа;
 
* если точный результат ''может'' быть представлен в <tex>p</tex> битах, то результатом округления будет точное значение числа;
 
* если точный результат ''не может'' быть представлен в <tex>p</tex> битах, то результатом округления будет ближайшее <tex>p</tex>-битное значение.
 
* если точный результат ''не может'' быть представлен в <tex>p</tex> битах, то результатом округления будет ближайшее <tex>p</tex>-битное значение.
 
}}
 
}}
  
'''Например''', в 4-битной арифметике произведение <tex>111 \times 101 = 100011</tex> будет округлено в <tex>1.001 \times 2^5</tex>.
+
Например, в 4-битной арифметике произведение <tex>111 \times 101 = 100011</tex> будет округлено в <tex>1.001 \times 2^5</tex>.
  
 
При вычислении результата может возникнуть ситуация, когда значение попадает в точности между двумя соседними <tex>p</tex>-битными значениями. Тогда требуется определить правило поведения в таком случае. Рассмотрим некоторые из них.
 
При вычислении результата может возникнуть ситуация, когда значение попадает в точности между двумя соседними <tex>p</tex>-битными значениями. Тогда требуется определить правило поведения в таком случае. Рассмотрим некоторые из них.
Строка 63: Строка 62:
 
{{Определение
 
{{Определение
 
|definition=
 
|definition=
'''Округление до ближайшего четного''' (англ. ''round-to-even'') - правило, при котором округление в вышеописанном случае производится к ближайшему <tex>p</tex>-битному четному значению.
+
'''Округление до ближайшего четного''' (англ. ''round-to-even'') {{---}}  правило, при котором округление в вышеописанном случае производится к ближайшему <tex>p</tex>-битному четному значению.
 
}}
 
}}
  
 
{{Определение
 
{{Определение
 
|definition=
 
|definition=
'''Округление к нулю''' (англ. ''round-toward-zero'') - правило, при котором округление в вышеописанном случае производится <tex>p</tex>-битному значению, находящемуся между точным значением и нулем, а также ближайшему к точному значению.
+
'''Округление к нулю''' (англ. ''round-toward-zero'') {{---}}  правило, при котором округление в вышеописанном случае производится <tex>p</tex>-битному значению, находящемуся между точным значением и нулем, а также ближайшему к точному значению.
 
}}
 
}}
  
'''Например''', в 4-битной арифметике число <tex>10011</tex> будет округлено до <tex>1.010 \times 2^4</tex> по первому правилу, и до <tex>1.001 \times 2^4</tex> по второму.
+
Например, в 4-битной арифметике число <tex>10011</tex> будет округлено до <tex>1.010 \times 2^4</tex> по первому правилу, и до <tex>1.001 \times 2^4</tex> по второму.
  
 
Стоит отметить, что стандарт IEEE 754 использует округление до ближайшего четного по умолчанию.
 
Стоит отметить, что стандарт IEEE 754 использует округление до ближайшего четного по умолчанию.
Строка 77: Строка 76:
 
Далее в этой статье символами <tex>\oplus, \ominus</tex> и <tex>\otimes</tex> будут обозначаться <tex>p</tex>-битные сложение, вычитание и умножение с точным округлением соответственно.
 
Далее в этой статье символами <tex>\oplus, \ominus</tex> и <tex>\otimes</tex> будут обозначаться <tex>p</tex>-битные сложение, вычитание и умножение с точным округлением соответственно.
  
Из-за округления данные операции теряют некоторые важные свойства, например, ассициативность: <tex>(1000 \oplus 0.011) \oplus 0.011 = 1000</tex>, но <tex>1000 \oplus(0.011 \oplus 0.011) = 1001</tex>.
+
Из-за округления данные операции теряют некоторые важные свойства, например, ассоциативность: <tex>(1000 \oplus 0.011) \oplus 0.011 = 1000</tex>, но <tex>1000 \oplus(0.011 \oplus 0.011) = 1001</tex>.
  
При анализе округления часто используют так называемый <tex>ulp</tex>.
+
При анализе округления часто используют так называемый <tex>\operatorname{ulp}</tex>.
 
{{Определение
 
{{Определение
 
|definition=
 
|definition=
'''ULP''' (англ. ''units in the last place'') - эффективная величина самого младшего (<tex>p</tex>-ого) бита.
+
'''ULP''' (англ. ''units in the last place'') {{---}}  эффективная величина самого младшего (<tex>p</tex>-ого) бита.
 
}}
 
}}
'''Например''', <tex>ulp(-1100) = 1, ulp(1) = 0.001</tex> в <tex>p</tex>-битной арифметике.
+
Например, <tex>\operatorname{ulp}(-1100) = 1, \operatorname{ulp}(1) = 0.001</tex> в <tex>p</tex>-битной арифметике.
  
Так же полезной нотацией является <tex> err(a \circledast b) </tex>, которая обозначает ошибку округлении результата выполнения операции <tex>\circledast</tex>. Отметим, что если <tex>ulp</tex> всегда беззнаковая величина, то <tex>err</tex> может иметь знак. Для базовых операций (сложение, вычитание, умножение) <tex> a \circledast b = a \ast b + err(a \circledast b)</tex>, и точное округление гарантирует, что <tex> |err(a \circledast b)| \leqslant \frac{1}{2}ulp(a \circledast b)</tex>.
+
Так же полезной нотацией является <tex> \operatorname{err}(a \circledast b) </tex>, которая обозначает ошибку округлении результата выполнения операции <tex>\circledast</tex>. Отметим, что если <tex>\operatorname{ulp}</tex> всегда беззнаковая величина, то <tex>\operatorname{err}</tex> может иметь знак. Для базовых операций (сложение, вычитание, умножение) <tex> a \circledast b = a \ast b + \operatorname{err}(a \circledast b)</tex>, и точное округление гарантирует, что <tex> |\operatorname{err}(a \circledast b)| \leqslant \frac{1}{2}\operatorname{ulp}(a \circledast b)</tex>.
  
 
==Свойства==
 
==Свойства==
Строка 94: Строка 93:
  
 
{{Лемма
 
{{Лемма
 +
|id=
 +
lemma1
 
|statement=
 
|statement=
 
Пусть <tex>a \oplus b = a + b + err(a \oplus b) </tex>.
 
Пусть <tex>a \oplus b = a + b + err(a \oplus b) </tex>.
Строка 102: Строка 103:
 
Как следствие, верна следующая лемма:
 
Как следствие, верна следующая лемма:
 
{{Лемма
 
{{Лемма
 +
|id=
 +
lemma2
 
|statement=
 
|statement=
 
Ошибка округления <tex>err(a \oplus b)</tex> может быть представлена в <tex>p</tex> битах.
 
Ошибка округления <tex>err(a \oplus b)</tex> может быть представлена в <tex>p</tex> битах.
 
}}
 
}}
 +
  
  
 
{{Лемма
 
{{Лемма
 +
|id=
 +
lemma3
 
|statement=
 
|statement=
 
Пусть <tex>|a + b| \leqslant |b|</tex> и <tex>|a + b| \leqslant |a|</tex>. Тогда <tex>a \oplus b = a + b</tex>. (Аналогично для вычитания).
 
Пусть <tex>|a + b| \leqslant |b|</tex> и <tex>|a + b| \leqslant |a|</tex>. Тогда <tex>a \oplus b = a + b</tex>. (Аналогично для вычитания).
Строка 114: Строка 120:
  
 
{{Лемма
 
{{Лемма
 +
|id=
 +
lemma4
 
|statement=
 
|statement=
 
Пусть <tex>b \in \left [ \frac{a}{2}, 2a \right ]</tex>. Тогда  <tex>a \ominus b = a - b</tex>.
 
Пусть <tex>b \in \left [ \frac{a}{2}, 2a \right ]</tex>. Тогда  <tex>a \ominus b = a - b</tex>.
Строка 120: Строка 128:
 
==Алгоритмы суммирования==
 
==Алгоритмы суммирования==
 
===Простое суммирование===
 
===Простое суммирование===
Важной базовой операцией во всех алгоритмах, основанных на представлении чисел в виде расширений, является сумма двух <tex>p</tex>-битных величин, результатом которой является расширение длины два. Есть два алгоритма для выполнения этой задачи - алгоритмы Деккера и Кнута.
+
Важной базовой операцией во всех алгоритмах, основанных на представлении чисел в виде разложений, является сумма двух <tex>p</tex>-битных величин, результатом которой является разложение длины два. Есть два алгоритма для выполнения этой задачи - алгоритмы Деккера и Кнута.
 
====Сумма по Деккеру====
 
====Сумма по Деккеру====
 
{{Теорема
 
{{Теорема
Строка 126: Строка 134:
 
Dekker
 
Dekker
 
|statement=
 
|statement=
Пусть <tex>a</tex> и <tex>b</tex> есть <tex>p</tex>-битные числа, причем <tex>|a| \geqslant |b|</tex>. Тогда следующий алгоритм вернет неперекрывающееся расширение <tex>x + y</tex> такое, что <tex>a + b = x + y</tex>, где <tex>x</tex> - приближение (аппроксимация) суммы <tex>a + b</tex>, а <tex>y</tex> представляет собой ошибку округления при вычислении <tex>x</tex>. (Аналогично для вычитания).
+
Пусть <tex>a</tex> и <tex>b</tex> есть <tex>p</tex>-битные числа, причем <tex>|a| \geqslant |b|</tex>. Тогда следующий алгоритм вернет неперекрывающееся разложение <tex>x + y</tex> такое, что <tex>a + b = x + y</tex>, где <tex>x</tex> - приближение (аппроксимация) суммы <tex>a + b</tex>, а <tex>y</tex> представляет собой ошибку округления при вычислении <tex>x</tex>. (Аналогично для вычитания).
 
}}
 
}}
  
Строка 157: Строка 165:
 
<br clear="all" />
 
<br clear="all" />
  
Проблема с использованием этой процедуры заключается в требовании <tex>|a| \geqslant |b|</tex>. Если это заранее не известно, то необходимо выполнить сравнение перед ее вызовом. В большинстве С компиляторов, возможно, самым быстрым переносимым способом реализовать эту проверку является выражение <tex>if ((a > b) == (a > -b))</tex>. На эту проверку уйдет некоторое время, но увеличение времени может быть на удивление большим из-за современных процессоров с суперскалярными и конвейерными архитектурами, в которых вызов условного оператора может сбросить ветку предсказаний.
+
<wikitex>Проблема с использованием этой процедуры заключается в требовании <tex>|a| \geqslant |b|</tex>. Если это заранее не известно, то необходимо выполнить сравнение перед ее вызовом. В большинстве С компиляторов, возможно, самым быстрым переносимым способом реализовать эту проверку является выражение <tex>if ((a > b) == (a > -b))</tex>. На эту проверку уйдет некоторое время, но увеличение времени может быть на удивление большим из-за современных процессоров с суперскалярными и конвейерными архитектурами, в которых вызов условного оператора может сбросить ветку предсказаний.
Это объяснение лишь гипотетическое и зависит от машины, но алгоритм TwoSum, что будет описан ниже, избегает этого сравнения посредством  трех дополнительных операций, что обычно на практике даже быстрее. Конечно же, FastTwoSum все же быстрее, если результат сравнения известен ''априори''.
+
Это объяснение лишь гипотетическое и зависит от машины, но алгоритм $\mathrm {TwoSum}$, что будет описан ниже, избегает этого сравнения посредством  трех дополнительных операций, что обычно на практике даже быстрее. Конечно же, $\mathrm {FastTwoSum}$ все же быстрее, если результат сравнения известен ''априори''.
 +
</wikitex>
  
 
====Сумма по Кнуту====
 
====Сумма по Кнуту====
Строка 165: Строка 174:
 
Knuth
 
Knuth
 
|statement=
 
|statement=
Пусть <tex>a</tex> и <tex>b</tex> есть <tex>p</tex>-битные числа, причем <tex>p > 3</tex>. Тогда следующий алгоритм вернет неперекрывающееся расширение <tex>x + y</tex> такое, что <tex>a + b = x + y</tex>, где <tex>x</tex> - приближение (аппроксимация) суммы <tex>a + b</tex>, а <tex>y</tex> представляет собой ошибку округления при вычислении <tex>x</tex>.
+
Пусть <tex>a</tex> и <tex>b</tex> есть <tex>p</tex>-битные числа, причем <tex>p > 3</tex>. Тогда следующий алгоритм вернет неперекрывающееся разложение <tex>x + y</tex> такое, что <tex>a + b = x + y</tex>, где <tex>x</tex> {{---}}  приближение (аппроксимация) суммы <tex>a + b</tex>, а <tex>y</tex> представляет собой ошибку округления при вычислении <tex>x</tex>.
 
}}
 
}}
  
Строка 179: Строка 188:
 
  <tex>7</tex>  <tex>return (x, y)</tex>
 
  <tex>7</tex>  <tex>return (x, y)</tex>
  
<wikitex>
+
<wikitex>Пример работы последнего алгоритма, когда $|a| < |b|$ и $|a| < |x|$. Сумма $11.11 + 1101$ будет представлена в виде разложение $10000 + 0.11$.
Пример работы последнего алгоритма, когда $|a| < |b|$ и $|a| < |x|$. Сумма $11.11 + 1101$ будет представлена в виде расширения $10000 + 0.11$.
 
 
[[Файл:Adaptive_10.jpg|слева]]
 
[[Файл:Adaptive_10.jpg|слева]]
 
<br clear="all" />
 
<br clear="all" />
Строка 186: Строка 194:
 
{{Лемма
 
{{Лемма
 
|statement=
 
|statement=
Пусть $x$ и $y$ - величины, возращенные алгоритмом $FastTwoSum$ или $TwoSum$. На машине с округлением до ближайшего четного $x$ и $y$ - несмежные.
+
Пусть $x$ и $y$ {{---}}  величины, возращенные алгоритмом $\mathrm {FastTwoSum}$ или $\mathrm {TwoSum}$. На машине с округлением до ближайшего четного $x$ и $y$ {{---}}  несмежные.
 
}}
 
}}
 
</wikitex>
 
</wikitex>
  
===Суммирование расширений===
+
===Суммирование разложений===
Базируясь на алгоритмах сложения двух <tex>p</tex>-битных чисел, описанных выше, можно предложить алгоритмы суммирования расширений. В этой статье их будет предложено два: ExpansionSum и FastExpansionSum. Первый алгоритм прибавляет к m-элементному расширению n-элементное расширение за время <tex>O(mn)</tex>, в то время, как последний алгоритм делает это за <tex>O(n + m)</tex>.
+
<wikitex>Базируясь на алгоритмах сложения двух <tex>p</tex>-битных чисел, описанных выше, можно предложить алгоритмы суммирования разложений. В этой статье их будет предложено два: $\mathrm {ExpansionSum}$ и $\mathrm {FastExpansionSum}$. Первый алгоритм прибавляет к $m$-элементному разложению $n$-элементное разложение за время <tex>O(mn)</tex>, в то время, как последний алгоритм делает это за <tex>O(n + m)</tex>.
 +
 
 +
Несмотря на такое различие в асимптотике, первый алгоритм на практике может оказаться быстрее на  разложениях, чей размер мал и фиксирован, потому что программные циклы могут быть полностью развернуты, а косвенные расходы времени исчезают, так как можно отказаться от использования массива). Линейный же алгоритм имеет определенные условия, при которых подобные оптимизации невозможны.  
  
Несмотря на такое различие в асимптотике, первый алгоритм на практике может оказаться быстрее на расширениях, чей размер мал и фиксирован, потому что программные циклы могут быть полностью развернуты, а косвенные расходы времени исчезают, так как можно отказаться от использования массива). Линейный же алгоритм имеет определенные условия, при которых подобные оптимизации невозможны.  
+
ExpansionSum имеет свойство, что если входные данные были неперекрывающимися (несмежными), то и на выходе мы получим  разложения с соответствующим свойством.  
  
ExpansionSum имеет свойство, что если входные данные были неперекрывающимися (несмежными), то и на выходе мы получим расширения с соотв.свойством.  
+
FastExpansionSum имеет несколько серьезных недостатков. Первый из них - алгоритм не сохраняет свойств неперекрываемости/несмежности. Второй - алгоритм основывается на округлении до ближайшего четного, что делает его непереносимым.
  
FastExpansionSum имеет несколько серьезных недостатков.Первый из них - алгоритм не сохраняет свойств неперекрываемости/несмежности. Второй - алгоритм основывается на округлении до ближайшего четного, что делает его непереносимым.
+
Как правило, общим недостатком алгоритмов работы с разложениями является то, что в разложении на выходе могут быть нулевые элементы, даже если в исходном разложении их не было. Например, если подать на вход разложениея $1111 + 0.0101$ и $1100 + 0.11$, то результатом будет $11100 + 0 + 0 + 0.0001$. К счастью, алгоритмы, описанные в этой статье хорошо справляются с этой проблемой.
<wikitex>
 
Как правило, общим недостатком алгоритмов работы с расширениями является то, что в расширении на выходе могут быть нулевые элементы, даже если в исходном расширении их не было. Например, если подать на вход расширениея $1111 + 0.0101$ и $1100 + 0.11$, то результатом будет $11100 + 0 + 0 + 0.0001$. К счастью, алгоритмы, описанные в этой статье хорошо справляются с этой проблемой.
 
 
</wikitex>
 
</wikitex>
  
Перед описанием алгоритмов суммирования расширений, приведем алгоритм прибавления к расширению <tex>p</tex>-битного числа.
+
Перед описанием алгоритмов суммирования разложений, приведем алгоритм прибавления к разложению <tex>p</tex>-битного числа.
  
====Grow expansion====
+
====Сумма разложения и числа====
<wikitex>
+
<wikitex>{{Теорема
{{Теорема
 
 
|statement=
 
|statement=
Пусть $e = \sum^{m}_{i=1}e_i$ - неперекрывающееся $m$-компонентное расширение; $b$ - $p$-битное число, где $p \geqslant 3$. Предполагается, что $e_1, e_2, \dots, e_m$ отсортированы в '''возрастающем''' порядке, причем все компоненты ненулевые. Тогда следующий алгоритм вернет такое расширение $h$, что $h = \sum^{m + 1}_{i=1}h_i = e + b$, где компоненты $h$ также отсортированы в возрастающем порядке и не равны нулю. Алгоритм также сохраняет свойство несмежности/неперекрываемости.
+
Пусть $e = \sum^{m}_{i=1}e_i$ - неперекрывающееся $m$-компонентное разложение; $b$ - $p$-битное число, где $p \geqslant 3$. Предполагается, что $e_1, e_2, \dots, e_m$ отсортированы в '''возрастающем''' порядке, причем все компоненты ненулевые. Тогда следующий алгоритм вернет такое разложение $h$, что $h = \sum^{m + 1}_{i=1}h_i = e + b$, где компоненты $h$ также отсортированы в возрастающем порядке и не равны нулю. Алгоритм также сохраняет свойство несмежности / неперекрываемости.
 
}}
 
}}
  
Строка 218: Строка 225:
 
  $3$      $(Q_i, h_i) \Leftarrow TwoSum(Q_{i-1}, e_i)$
 
  $3$      $(Q_i, h_i) \Leftarrow TwoSum(Q_{i-1}, e_i)$
 
  $4$  $h_{m+1} \Leftarrow Q_m$
 
  $4$  $h_{m+1} \Leftarrow Q_m$
  $5$  $return h$
+
  $5$  $return$ $h$
  
 
$Q_i$ - приближенное значение суммы $b$ и первых $i$ компонент $e$.
 
$Q_i$ - приближенное значение суммы $b$ и первых $i$ компонент $e$.
Строка 227: Строка 234:
 
{{Лемма
 
{{Лемма
 
|statement=
 
|statement=
Каждая их первых $m$ компонент расширения $h$ не больше чем соответствующая компонента из $e$, то есть $|h_i| \leqslant |e_i|$. Более того, $|h_1| \leqslant |b|$.
+
Каждая их первых $m$ компонент разложения $h$ не больше чем соответствующая компонента из $e$, то есть $|h_i| \leqslant |e_i|$. Более того, $|h_1| \leqslant |b|$.
 
}}
 
}}
 
</wikitex>
 
</wikitex>
  
 
+
====Простое суммирование====
====Expansion sum====
+
<wikitex>{{Теорема
<wikitex>
 
{{Теорема
 
 
|statement=
 
|statement=
Пусть $e = \sum^{m}_{i=1}e_i$ и $f = \sum^{n}_{i=1}f_i$ - неперекрывающиеся $m$- и $n$-компонентные расширения, компоненты которых - $p$-битные числа, где $p \geqslant 3$. Предполагается, что компоненты обоих расширений отсортированы в '''возрастающем''' порядке, причем все компоненты ненулевые. Тогда следующий алгоритм вернет такое расширение $h$, что $h = \sum^{m + n}_{i=1}h_i = e + f$, где компоненты $h$ также отсортированы в возрастающем порядке и не равны нулю. Алгоритм также сохраняет свойство несмежности/неперекрываемости.
+
Пусть $e = \sum^{m}_{i=1}e_i$ и $f = \sum^{n}_{i=1}f_i$ - неперекрывающиеся $m$- и $n$-компонентные разложения, компоненты которых - $p$-битные числа, где $p \geqslant 3$. Предполагается, что компоненты обоих разложений отсортированы в '''возрастающем''' порядке, причем все компоненты ненулевые. Тогда следующий алгоритм вернет такое разложение $h$, что $h = \sum^{m + n}_{i=1}h_i = e + f$, где компоненты $h$ также отсортированы в возрастающем порядке и не равны нулю. Алгоритм также сохраняет свойство несмежности/неперекрываемости.
 
}}
 
}}
  
Строка 251: Строка 256:
 
<br clear="all" />
 
<br clear="all" />
  
====Fast expansion sum====
+
====Быстрое суммирование====
<wikitex>
+
<wikitex>В отличие от $\mathrm {ExpansionSum}$, $\mathrm {FastExpansionSum}$ не сохраняет свойства неперекрываемости и несмежности, но гарантирует, что если если входное разложение было ''строго неперекрывающимся'', то на выходе получится разложение с таким же свойством.
В отличие от $ExpansionSum$, $FastExpansionSum$ не сохраняет свойства неперекрываемости и несмежности, но гарантирует, что если если входное расширение было ''строго неперекрывающимся'', то на выходе получится расширение с таким же свойством.
 
  
 
{{Определение
 
{{Определение
 
|definition=
 
|definition=
Расширение называется '''строго неперекрывающимся''', если если его компоненты попарно неперекрываются, ни одна компонента не смежна никаким двум другим, а также любая пара смежных компонент состоит из степеней двойки.
+
Разложение называется '''строго неперекрывающимся''', если если его компоненты попарно неперекрываются, ни одна компонента не смежна никаким двум другим, а также любая пара смежных компонент состоит из степеней двойки.
 
}}
 
}}
'''Например''', $11000 + 11$ и $10000 + 1000 + 10 + 1$ - строго неперекрывающиеся расширения, а $11100 + 11$ и $100 + 10 + 1$ - нет.
+
''Например'', $11000 + 11$ и $10000 + 1000 + 10 + 1$ {{---}}  строго неперекрывающиеся разложения, а $11100 + 11$ и $100 + 10 + 1$ {{---}}  нет.
  
Для расширения эта характеристика означает, что ноль в записи расширения должен появляться ''как минимум'' каждые $p + 1$ бит. Например, расширение, каждая компонента которого есть $p$-битное число, и максимальная величина которой равна $1111$, может быть не больше $1111.01111011110\dots$.
+
Для разложения эта характеристика означает, что ноль в записи разложения должен появляться ''как минимум'' каждые $p + 1$ бит. Например, разложение, каждая компонента которого есть $p$-битное число, и максимальная величина которой равна $1111$, может быть не больше $1111.01111011110\dots$.
  
Каждое несмежное расширение является строго неперекрывающимся, а также каждое строго неперекрывающееся расширение есть неперекрывающееся. В общем случае, обратное не верно.
+
Каждое несмежное разложение является строго неперекрывающимся, а также каждое строго неперекрывающееся разложение есть неперекрывающееся. В общем случае, обратное не верно.
  
Утверждается, что после предположения о том, что все расширения строго неперекрывающиеся, алгоритм, приведенный ниже, корректно работает на машинах с округлением до ближайшего четного.
+
Утверждается, что после предположения о том, что все разложения строго неперекрывающиеся, алгоритм, приведенный ниже, корректно работает на машинах с округлением до ближайшего четного.
  
 
{{Теорема
 
{{Теорема
 
|statement=
 
|statement=
Пусть $e = \sum^{m}_{i=1}e_i$ и $f = \sum^{n}_{i=1}f_i$ - строго неперекрывающиеся $m$- и $n$-компонентные расширения, компоненты которых - $p$-битные числа, где $p \geqslant 4$. Предполагается, что компоненты обоих расширений отсортированы в '''возрастающем''' порядке, причем все компоненты ненулевые. Тогда следующий алгоритм,запущенный на машине с правилом округления до ближайшего четного, вернет такое расширение $h$, что $h = \sum^{m + n}_{i=1}h_i = e + f$, где компоненты $h$ также отсортированы в возрастающем порядке и не равны нулю. Алгоритм также сохраняет свойство несмежности/неперекрываемости.
+
Пусть $e = \sum^{m}_{i=1}e_i$ и $f = \sum^{n}_{i=1}f_i$ {{---}}  строго неперекрывающиеся $m$- и $n$-компонентные разложения, компоненты которых {{---}}  $p$-битные числа, где $p \geqslant 4$. Предполагается, что компоненты обоих разложений отсортированы в '''возрастающем''' порядке, причем все компоненты ненулевые. Тогда следующий алгоритм,запущенный на машине с правилом округления до ближайшего четного, вернет такое разложение $h$, что $h = \sum^{m + n}_{i=1}h_i = e + f$, где компоненты $h$ также отсортированы в возрастающем порядке и не равны нулю. Алгоритм также сохраняет свойство несмежности / неперекрываемости.
 
}}
 
}}
  
Строка 294: Строка 298:
 
* [http://en.wikipedia.org/wiki/Rounding Rounding]
 
* [http://en.wikipedia.org/wiki/Rounding Rounding]
 
* [http://www.google.ru/url?sa=t&rct=j&q=Adaptive%2Bprecision%2Barithmetic&source=web&cd=4&ved=0CDgQFjAD&url=http%3A%2F%2Fwww.cs.berkeley.edu%2F~jrs%2Fpapers%2Frobustr.pdf&ei=fZigTp6xHoWa-waGxKSOBQ&usg=AFQjCNGl9T4V-Y02Wi99TjgSDLFotO5xeg&sig2=pQ5SPC_lRGCtwoBZhMhGhA  Jonathan Richard Shewchuk, "Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates"]
 
* [http://www.google.ru/url?sa=t&rct=j&q=Adaptive%2Bprecision%2Barithmetic&source=web&cd=4&ved=0CDgQFjAD&url=http%3A%2F%2Fwww.cs.berkeley.edu%2F~jrs%2Fpapers%2Frobustr.pdf&ei=fZigTp6xHoWa-waGxKSOBQ&usg=AFQjCNGl9T4V-Y02Wi99TjgSDLFotO5xeg&sig2=pQ5SPC_lRGCtwoBZhMhGhA  Jonathan Richard Shewchuk, "Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates"]
 +
 +
 +
[[Категория: Вычислительная геометрия]]

Текущая версия на 19:17, 4 сентября 2022

Мотивация

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

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

Представление чисел

Большинство современных процессоров поддерживают числа с плавающей точкой в форме [math] \pm \mathrm{significand} \times 2^{\mathrm{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].

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

Например, числа [math]1100[/math] и [math]-10.1[/math] не пересекаются, а [math]110[/math] и [math]10[/math] — пересекаются.


Определение:
Два числа [math]x[/math] и [math]y[/math] называются смежными (англ. adjacent), если они перекрываются, если [math]x[/math] и [math]2y[/math] перекрываются, или [math]2x[/math] и [math]y[/math] перекрываются.


Например, числа [math]1100[/math] и [math]11[/math] — смежные, а [math]1100[/math] и [math]1000[/math] — нет.

Иногда для использовании точной арифметики может понадобиться больше, чем [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] — большая). Далее будут рассматриваться именно такая их форма.

<wikitex>Стоит отметить, что число может быть представлено несколькими возможными неперекрывающимися разложениями: $1100 + -10.1 = 1001 + 0.1 = 1000 + 1 + 0.1$. </wikitex> Неперекрывающиеся разложения нужны, например, для того, чтобы быстро вычислять знак выражения (смотрим знак большей по размеру компоненты), или для грубой оценки значения всего разложения (берем большую по величине компоненту).

Округление

Все алгоритмы, представленные в этой статье, предполагают, что сложение, вычитание и умножение производятся с точным округлением. Предполагается, что числа представляются в [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].

При анализе округления часто используют так называемый [math]\operatorname{ulp}[/math].

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

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

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


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

Лемма:
Ошибка округления [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]. (Аналогично для вычитания).


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

[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]FastTwoDiff(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]

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

В первом примере [math]a[/math] и [math]b[/math] имеют один знак.

Рисунок к теореме Деккера


Во втором примере [math]a[/math] и [math]b[/math] имеют противоположные знаки, и [math]|b| \gt \frac{|a|}{2}[/math].

Adaptive 3.jpg


<wikitex>Проблема с использованием этой процедуры заключается в требовании [math]|a| \geqslant |b|[/math]. Если это заранее не известно, то необходимо выполнить сравнение перед ее вызовом. В большинстве С компиляторов, возможно, самым быстрым переносимым способом реализовать эту проверку является выражение [math]if ((a \gt b) == (a \gt -b))[/math]. На эту проверку уйдет некоторое время, но увеличение времени может быть на удивление большим из-за современных процессоров с суперскалярными и конвейерными архитектурами, в которых вызов условного оператора может сбросить ветку предсказаний. Это объяснение лишь гипотетическое и зависит от машины, но алгоритм $\mathrm {TwoSum}$, что будет описан ниже, избегает этого сравнения посредством трех дополнительных операций, что обычно на практике даже быстрее. Конечно же, $\mathrm {FastTwoSum}$ все же быстрее, если результат сравнения известен априори. </wikitex>

Сумма по Кнуту

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

Псевдокод:

[math]TwoSum(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]a_{virtual} \Leftarrow x \ominus b_{virtual}[/math]
[math]4[/math]   [math]b_{roundoff} \Leftarrow b \ominus b_{virtual}[/math]
[math]5[/math]   [math]a_{roundoff} \Leftarrow a \ominus a_{virtual}[/math]
[math]6[/math]   [math]y \Leftarrow a_{roundoff} \oplus b_{roundoff}[/math]
[math]7[/math]   [math]return (x, y)[/math]

<wikitex>Пример работы последнего алгоритма, когда $|a| < |b|$ и $|a| < |x|$. Сумма $11.11 + 1101$ будет представлена в виде разложение $10000 + 0.11$.

Adaptive 10.jpg


Лемма:
Пусть $x$ и $y$ — величины, возращенные алгоритмом $\mathrm {FastTwoSum}$ или $\mathrm {TwoSum}$. На машине с округлением до ближайшего четного $x$ и $y$ — несмежные.

</wikitex>

Суммирование разложений

<wikitex>Базируясь на алгоритмах сложения двух [math]p[/math]-битных чисел, описанных выше, можно предложить алгоритмы суммирования разложений. В этой статье их будет предложено два: $\mathrm {ExpansionSum}$ и $\mathrm {FastExpansionSum}$. Первый алгоритм прибавляет к $m$-элементному разложению $n$-элементное разложение за время [math]O(mn)[/math], в то время, как последний алгоритм делает это за [math]O(n + m)[/math].

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

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

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

Как правило, общим недостатком алгоритмов работы с разложениями является то, что в разложении на выходе могут быть нулевые элементы, даже если в исходном разложении их не было. Например, если подать на вход разложениея $1111 + 0.0101$ и $1100 + 0.11$, то результатом будет $11100 + 0 + 0 + 0.0001$. К счастью, алгоритмы, описанные в этой статье хорошо справляются с этой проблемой. </wikitex>

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

Сумма разложения и числа

<wikitex>
Теорема:
Пусть $e = \sum^{m}_{i=1}e_i$ - неперекрывающееся $m$-компонентное разложение; $b$ - $p$-битное число, где $p \geqslant 3$. Предполагается, что $e_1, e_2, \dots, e_m$ отсортированы в возрастающем порядке, причем все компоненты ненулевые. Тогда следующий алгоритм вернет такое разложение $h$, что $h = \sum^{m + 1}_{i=1}h_i = e + b$, где компоненты $h$ также отсортированы в возрастающем порядке и не равны нулю. Алгоритм также сохраняет свойство несмежности / неперекрываемости.

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

$GrowExpansion(e, b)$
$1$  $Q_0 \Leftarrow b$
$2$  $for$ $i = 1 \dots m$
$3$       $(Q_i, h_i) \Leftarrow TwoSum(Q_{i-1}, e_i)$
$4$  $h_{m+1} \Leftarrow Q_m$
$5$  $return$ $h$

$Q_i$ - приближенное значение суммы $b$ и первых $i$ компонент $e$.

Рисунок к теореме Деккера


Лемма:
Каждая их первых $m$ компонент разложения $h$ не больше чем соответствующая компонента из $e$, то есть $

</wikitex>

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

<wikitex>
Теорема:
Пусть $e = \sum^{m}_{i=1}e_i$ и $f = \sum^{n}_{i=1}f_i$ - неперекрывающиеся $m$- и $n$-компонентные разложения, компоненты которых - $p$-битные числа, где $p \geqslant 3$. Предполагается, что компоненты обоих разложений отсортированы в возрастающем порядке, причем все компоненты ненулевые. Тогда следующий алгоритм вернет такое разложение $h$, что $h = \sum^{m + n}_{i=1}h_i = e + f$, где компоненты $h$ также отсортированы в возрастающем порядке и не равны нулю. Алгоритм также сохраняет свойство несмежности/неперекрываемости.

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

$ExpansionSum(e, f)$
$1$  $h \Leftarrow e$
$2$  $for$ $i = 1 \dots n$
$3$       $\langle h_i, h_{i+1}, \dots, h_{i+m} \rangle \Leftarrow GrowExpansion(\langle h_i, h_{i+1}, \dots, h_{i+m-1} \rangle, f_i)$
$4$  $return$ $h$

</wikitex>

Рисунок к теореме Деккера


Быстрое суммирование

<wikitex>В отличие от $\mathrm {ExpansionSum}$, $\mathrm {FastExpansionSum}$ не сохраняет свойства неперекрываемости и несмежности, но гарантирует, что если если входное разложение было строго неперекрывающимся, то на выходе получится разложение с таким же свойством.


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

Например, $11000 + 11$ и $10000 + 1000 + 10 + 1$ — строго неперекрывающиеся разложения, а $11100 + 11$ и $100 + 10 + 1$ — нет.

Для разложения эта характеристика означает, что ноль в записи разложения должен появляться как минимум каждые $p + 1$ бит. Например, разложение, каждая компонента которого есть $p$-битное число, и максимальная величина которой равна $1111$, может быть не больше $1111.01111011110\dots$.

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

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

Теорема:
Пусть $e = \sum^{m}_{i=1}e_i$ и $f = \sum^{n}_{i=1}f_i$ — строго неперекрывающиеся $m$- и $n$-компонентные разложения, компоненты которых — $p$-битные числа, где $p \geqslant 4$. Предполагается, что компоненты обоих разложений отсортированы в возрастающем порядке, причем все компоненты ненулевые. Тогда следующий алгоритм,запущенный на машине с правилом округления до ближайшего четного, вернет такое разложение $h$, что $h = \sum^{m + n}_{i=1}h_i = e + f$, где компоненты $h$ также отсортированы в возрастающем порядке и не равны нулю. Алгоритм также сохраняет свойство несмежности / неперекрываемости.

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

$ExpansionSum(e, f)$
$1$  $Merge$ $e$ $and$ $f$ $into$ $a$ $single$ $sequence$ $g$, $in$ $order$ $of$ $nondecreasing$ $magnitude$ $(possibly$ $with$ $interspersed$ $zeros)$
$2$  $(Q_i, h_1) \Leftarrow FastTwoSum(g_2, g_1)$
$3$  $for$ $i = 3 \dots (n+m)$
$4$      $(Q_i, h_{i-1}) \Leftarrow TwoSum(Q_{i-1}, g_i)$
$5$  $h_{n+m} \Leftarrow Q_{m+n}$
$6$  $return$ $h$
Adaptive 8.jpg



На рисунке ниже показан конкретный пример работы алгоритма

Adaptive 9.jpg


</wikitex>

Ссылки