Длинная арифметика

Длинная арифметика — в вычислительной технике операции (сложение, вычитание, умножение, деление, возведение в степень, элементарные функции) над числами, разрядность которых превышает длину машинного слова данной вычислительной машины. Эти операции реализуются не аппаратно, а программно, используя базовые аппаратные средства работы с числами меньших порядков. Частный случай — арифметика произвольной точности — относится к арифметике, в которой длина чисел ограничена только объёмом доступной памяти.

Основные потребители

• Компьютеры низкой разрядности, микроконтроллеры. Например, микроконтроллеры серии AVR имеют 8-битный регистр и 10-битный АЦП — так что при обработке информации с АЦП без длинной арифметики не обойтись.
• Криптография. Большинство систем шифрования данных, а также систем цифровой подписи данных используют целочисленную арифметику по модулю m, где m — очень большое положительное натуральное чисто, не обязательно простое. Для примера RSA[1], Rabin[2] или Эль Гамаль[3] требуют наличие эффективных методов для операций перемножения и возведения в степень для чисел, имеющих порядки 10³⁰⁹.
• Математическое и финансовое ПО, требующее, чтобы результат вычисления на компьютере совпал до последнего разряда с результатом вычисления на бумаге. В частности, калькулятор Windows (начиная с Windows 95). Калькулятор Windows Vista и далее проводит четыре арифметических действия с намного большей точностью, чем позволяет процессор x86. Для обычных научных и инженерных расчётов длинная арифметика применяется редко: ошибки во входных данных обычно намного больше, чем ошибки округления.
• Одна из излюбленных тем спортивного программирования. С распространением стандартных библиотек для работы с длинными числами постепенно сходит на нет.

Необходимые аппаратные средства для работы с длинной арифметикой

Строго говоря, для реализации арифметики произвольной точности от процессора требуется лишь косвенная адресация; в арифметике фиксированной точности можно обойтись даже без неё. Тем не менее, определённые функции процессора ускоряют длинную арифметику, одновременно упрощая её программирование.

• Флаг переноса. Операции «сложить/вычесть с переносом», «циклический сдвиг через бит переноса».
• Автоинкрементные и автодекрементные операции доступа к памяти.
• Умножение w·w→2w, деление 2w /w→w.

Реализация в языках программирования

Языки программирования имеют встроенные типы данных, размер которых, в основном, не превышает 64 бита (около 10¹⁹). Десятичная длинная арифметика была реализована в языке программирования АЛМИР-65 на ЭВМ МИР-1 и в языке программирования АНАЛИТИК на ЭВМ МИР-2. Для работы с большими числами, однако, существует довольно много готовых оптимизированных библиотек для длинной арифметики.

• GMP
• LibTomMath

Встроенные библиотеки работы с большими числами есть в Python и Java.

Алгоритмы

Алгоритмы умножения

Базовый

Представляет собой алгоритм по школьному методу «в столбик». Занимает время O(N*M), где N, M — размеры перемножаемых чисел. Его алгоритм подробно описан в книге [1]. Секция 4.3.1.

Умножение Карацубы

Этот алгоритм также описан в [1]. Секция 4.3.3, часть А[4]. Данный алгоритм представляет собой наиболее простую реализацию идеи разделения входных данных, которая стала базисной для нижеперечисленных алгоритмов.

Toom 3-Way

Идея впервые была предложена А. Л. Тоомом в [3], и заключается в разделении входных данных (множителей) на 3 равные части(для простоты все части считаем равными).
Рассмотрим данный алгоритм на следующем примере. Пусть дано два числа Y и X. Пусть определены операции над числами, размер которых равен 1/3 от размеров исходных чисел (см. рисунок). (предположим также, что числа занимают равную память и делятся ровно на 3 части, в противном случае дополним оба числа нулями до необходимых размеров.

Затем происходит отображение (параметризация) этих частей на полиномы 2 степени.

$X(t) = x2*t^2 + x1*t + x0$ $Y(t) = y2*t^2 + y1*t + y0$

Введем $b$, по определению такую, что значение многочленов$X(b) Y(b)$ соответственно равно входным числам $x$ и $y$. Для битового представления чисел это число равно возведению 2 в степень, равную длине каждой из частей в битах.

Также введем полином:
$W(t)=X(t)*Y(t)~~~~~~~~~~~~~~~~~~~~$ (1).
$W(t) = w4*t^4 + w3*t^3 + w2*t^2 + w1*t + w0$

После того как будут определены элементы $w[i]$ — коэффициенты многочлена$W(t)$, они будут использованы в целях получения $w=W(b)$, так как $x*y=X(b)*Y(b)=W(b)$. Размер коэффициентов $W[i]$ в 2 раза больше (по памяти), чем разбиение $x[i]$ или $y[i]$. Конечное число, равное произведению $x*y$ можно найти, складывая $W[i]$ со сдвигом, как показано на рисунке ниже.

Коэффициенты $w[i]$ могут быть вычислены следующим образом: $w4=x2*y2, w3=x2*y1+x1*y2, w2=x2*y0+x1*y1+x0*y2$ и так далее, но это потребует все 9 перемножений: $x[i]*y[j]$ для i, j=0,1,2, и будет эквивалентен простому перемножению. Вместо этого был использован иной подход. $X(t),Y(t)$ вычисляются в (1) в 5 точках при разных $t$.
В таблице, представленной ниже, приведены значения полиномов в равенстве (1)
$~~~t=0~~~~~~~~~~x0 * y0~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~W(0)~~~~~~~~~=~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~w0$
$~~~t=1~~~~~~~~~~~(x2+x1+x0)*(y2+y1+y0)~~~~~~~~~~~~~~~~~~~~~~W(1)~~~~~~~=~~~~w4 + w3 + w2 + w1 + w0$
$~~~t=-1~~~~~~~~(x2-x1+x0) * (y2-y1+y0)~~~~~~~~~~~~~~~~~~~~~~W(-1)~~~~~~=~~~w4 - w3 + w2 - w1 + w0$
$~~~t=2~~~~~~~~~~~~(4*x2+2*x1+x0) * (4*y2+2*y1+y0)~~~~W(2)~~~= 16*w4 + 8*w3 + 4*w2 + 2*w1 + w0$
$~~~t=inf~~~~~~x2 * y2~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~W(inf)~=~~w4$
Параметр t=inf условный. Он означает банальное равенство коэффициентов при $t^4$, тем самым мы полум значение $w4$ сразу. Данная система линейна относительно 5 неизвестных. При её разрешении мы получаем неизвестные $w[i]$. Далее получаем значение многочлена, как было описано выше.

Алгоритм Тоом-а для произвольного разбиения

Алгоритм Тоома для разделения входных чисел на n операндов эквивалентен описанному выше. В общем случае разделение двух одинаково длинных операндов на r частей приводит к вычислению значений в 2*r-1 точках. В качестве точек для решения системы, обычно берут 0, inf, +1, −1 and ±2^i, ±2^-i.

Перемножение методом Фурье

Данный алгоритм перемножения используется для больших и очень больших чисел. Описание данного алгоритма можно найти в различных источниках, в частности [1] секция 4.3.3 часть С[5], или Lipson глава 9. Далее приводится описание алгоритма.
Данный вид алгоритма перемножения двух чисел, за время $O(n*\log n)$ , где $n$ — количество значащих цифр в перемножаемых числах (предполагая их равными), приписывается Кули (Coolet) и Таки (Tukey) — 1965 г. Также есть предположения, что данный алгоритм был известен и раньше, но не имел большой важности до того как изобрели первые компьютеры. Вероятными кандидатами изобретения данного алгоритма также называют Рунг (Runge) и Кёниг (Konig) — 1924 г., а также Гаусса — 1805 г.
Пусть имеется число, представим его в виде многочлена, как мы делали ранее. Назовем этот многочлен $A(x)$. Также введем дискретное Фурье [6]преобразование многочлена как вектор [7], с координатами $\left (b_1, b_2, b_3, b_4, \dots b_{n-1}\right )$. Где $b[i]$ определены как $w^n[i]$ — комплексный корень $n$-ой степени из 1, не равный 1[8]. Дело в том, что комплексный корень из 1 определен с точностью до фазового множителя, количество этих корней — $n$. [9] Фурье преобразование применено для того, чтобы свертку[10] коэффициентов многочленов А и В: $(A(x)*B(x))$ — заменить на произведение их Фурье образов.
$F(A(x) * B(x)) = F (A) * F (B)$
$A(x) * B(x) = F^{-1} (F (A) * F (B))$,
где под умножением F (A) * F (B) подразумевается скалярное произведение векторов.
Предположим, что n есть степень двойки.
Нахождение F(A) производится рекурсивным методом (разделяй и властвуй). Идея заключается в разбиении исходного многочлена $A (x) = a_0x_0 + a_1x_1 + ... + a_{n-1}x_{n-1}$ на два многочлена,
$A0 (x) = a_0x_0 + a_2x_1 + a_4x_2 + \dots + a_{n-2}x_{\frac{n-2}{2}}$
$A1 (x) = a_1x_0 + a_3x_1 + a_5x_2 + \dots + a_{n-1}x_{\frac{n-2}{2}}$
Тогда $A(x) = A0 (x^2) + x * A1 (x^2)$.
Заметим, что среди всех чисел $\left (w^n_i\right )^2$ (0 <= i < $n$), только $\frac{n}{2}$ различных. Поэтому, ДПФ $A0$ и $A1$ будут $\frac{n}{2}$-элементными. А так как ДПФ A0 и A1 состоит из $\frac{n}{2}$ элементов, то комплексным корнем из 1 там будет корень степени $\frac{n}{2}$.
Значит, $A(w^n_k) = A0(w^n_{2k}) + w^n_k * A1(w^n_{2k}) = A0(w^\frac{n}{2}_k) + w^n_k * A1(w^\frac{n}{2}_k)$, где 0 <= k < n / 2 и
$A(w^n_{k+\frac{n}{2}}) = A0(w^n_{2*{k+\frac{n}{2}}}) + w^n_{k+\frac{n}{2}} * A1(w^n_{2*{k+\frac{n}{2}}}) = A0(w^{n/2}_{k+\frac{n}{2}}) + w^n_{k+\frac{n}{2}} * A1(w^\frac{n}{2}_{k+\frac{n}{2}}) = A0(w^\frac{n}{2}_k) - w^n_k * A1(w^\frac{n}{2}_k)$.
Мы использовали свойства комплексных чисел: различных корней степени n всего n. $w^n_{k+\frac{n}{2}} = w^n_k * e^{i*\frac{2\pi}{n} * \frac{n}{2}} = w^n_k * e^{i*\pi} = - w^n_k$.

Получаем рекурсивный алгоритм:
ДПФ одного элемента равно этому элементу
для нахождения ДПФ A разбиваем коэффициенты A на чётные и нечётные, получаем два многочлена A0 и A1, находим ДПФ для них, находим ДПФ A:
$b_k = b^0_k + w^n_k * b^1_k$
$b_{k + \frac{n}{2}} = b^0_k - w^n_k * b^1_k$
для 0 <= k < n / 2.
Существует доказательство следующего утверждения: Обратное ДПФ находится аналогично прямому ДПФ, за исключением того, что в качестве точек берутся точки, симметричные относительно вещественной оси, вместо тех, что использовались в прямом ДПФ преобразовании. Также необходимо все коэффициенты многочлена поделить на n

Литература

1. Donald E. Knuth, «The Art of Computer Programming», volume 2, «Seminumerical Algorithms», 3rd edition, Addison-Wesley, 1998
2. Dan Zuras, «On Squaring and Multiplying Large Integers», ARITH-11: IEEE Symposium on Computer Arithmetic, 1993, pp. 260 to 271.
3. ДАН СССР 150 (1963), 496—498