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

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

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

• Компьютеры низкой разрядности, микроконтроллеры. Например, микроконтроллеры серии AVR имеют 8-битный регистр и 10-битный АЦП — так что при обработке информации с АЦП без длинной арифметики не обойтись.
• Криптография. Большинство систем шифрования данных, а также систем цифрофой подписи данных используют целочисленную арифметику по модулю m, где m — очень большое положительное натуральное чисто, не обязательно простое. Для примера RSA[1], Rabin[2] или Эль Гамаль[3] требуют наличие эффективных методов для операций перемножения и возведения в степень для чисел, имеющих порядки 10³⁰⁹Шаблон:-1
• Математическое и финансовое ПО, требующее, чтобы результат вычисления на компьютере совпал до последнего разряда с результатом вычисления на бумаге. В частности, калькулятор 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_{0}x_{0}+a_{1}x_{1}+...+a_{n-1}x_{n-1}$ на два многочлена,
$A0(x)=a_{0}x_{0}+a_{2}x_{1}+a_{4}x_{2}+\dots +a_{n-2}x_{\frac {n-2}{2}}$
$A1(x)=a_{1}x_{0}+a_{3}x_{1}+a_{5}x_{2}+\dots +a_{n-1}x_{\frac {n-2}{2}}$
Тогда $A(x)=A0(x^{2})+x*A1(x^{2})$.
Заметим, что среди всех чисел $\left(w_{i}^{n}\right)^{2}$ (0 <= i < $n$), только ${\frac {n}{2}}$ различных. Поэтому, ДПФ $A0$ и $A1$ будут ${\frac {n}{2}}$-элементными. А так как ДПФ A0 и A1 состоит из ${\frac {n}{2}}$ элементов, то комплексным корнем из 1 там будет корень степени ${\frac {n}{2}}$.
Значит, $A(w_{k}^{n})=A0(w_{2k}^{n})+w_{k}^{n}*A1(w_{2k}^{n})=A0(w_{k}^{\frac {n}{2}})+w_{k}^{n}*A1(w_{k}^{\frac {n}{2}})$, где 0 <= k < n / 2 и
$A(w_{k+{\frac {n}{2}}}^{n})=A0(w_{2*{k+{\frac {n}{2}}}}^{n})+w_{k+{\frac {n}{2}}}^{n}*A1(w_{2*{k+{\frac {n}{2}}}}^{n})=A0(w_{k+{\frac {n}{2}}}^{n/2})+w_{k+{\frac {n}{2}}}^{n}*A1(w_{k+{\frac {n}{2}}}^{\frac {n}{2}})=A0(w_{k}^{\frac {n}{2}})-w_{k}^{n}*A1(w_{k}^{\frac {n}{2}})$.
Мы использовали свойства комплексных чисел: различных корней степени n всего n. $w_{k+{\frac {n}{2}}}^{n}=w_{k}^{n}*e^{i*{\frac {2\pi }{n}}*{\frac {n}{2}}}=w_{k}^{n}*e^{i*\pi }=-w_{k}^{n}$.

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

Литература

1. Donald E. Knuth, «The Art of Computer Programming», volume 2, «Seminumerical Algo-rithms», 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