Алгоритм Тоома — Кука

Алгоритм Тоома — Кука, иногда упоминаемый как Тоом-3 — это алгоритм умножения^[англ.] больших чисел, названный именами Андрея Леоновича Тоома^[англ.], предложившего новый алгоритм с низкой сложностью, и Стивена Кука, более ясно его описавшего.

Если даны два больших числа a и b, согласно алгоритму Тоома — Кука, нужно разбить a и b на k меньших частей каждое длиной l и осуществить операции над элементами. При росте k можно комбинировать часть операций умножения частей разбиения, сокращая тем самым общую сложность алгоритма. Произведение частей разбиения можно затем вычислить с помощью того же самого алгоритма Тоома — Кука рекурсивно. Термины «алгоритм Тоома-3» и «алгоритм Тоома — Кука» иногда ошибочно используются как синонимы, хотя Тоом-3 является лишь частным случаем алгоритма Тоома — Кука для k = 3.

Тоом-3 сокращает сложность с 9 умножений до 5 и работает за время $\Theta (n^{\log(5)/\log(3)})\approx \Theta (n^{1{,}46})$ . В общем случае Тоом-k работает за время $\Theta (c(k)n^{e})$ , где $e=\log {(2k-1)}/\log {(k)}$ , $n^{e}$ является временем, затрачиваемым на умножения подчастей, а c — время, затрачиваемое на сложения и умножение на малые константы ^[1]. Алгоритм Карацубы является частным случаем алгоритма Тоома — Кука, где число разбивается на две части. Он сокращает сложность с 4 умножений до 3, а потому работает за время $\Theta (n^{\log(3)/\log(2)})\approx \Theta (n^{1{,}58})$ . Обычное умножение в столбик эквивалентно Тоом-1 со сложностью $\Theta (n^{2})$ .

Хотя степень e может быть установлена как можно близкой к 1 путём увеличения k, функция c растёт очень быстро^[1]^[2]. Скорость роста для смешанных схем Тоома — Кука оставалась открытой проблемой к 2005-му году^[3]. Реализация, описанная Дональдом Кнутом, добивается сложности $\Theta (n2^{\sqrt {2\log n}}\log n)$ ^[4].

За счёт накладных расходов алгоритм Тоома-Кука для малых чисел медленнее умножения в столбик и потому обычно использовался для множителей среднего размера, пока не был обнаружен асимптотически более быстрый алгоритм Шёнхаге — Штрассена (со сложностью $\Theta (n\log n\log {\log n}))$ .

Тоом описал свой алгоритм в 1963 году, а Кук опубликовал улучшенный (асимптотически эквивалентный) алгоритм в тезисах своей диссертации в 1966-м году^[5].

Детали

В этой секции обсуждается работа алгоритма Тоома-k для любого заданного значения k и даётся упрощённое описание умножения Тоома — Кука многочленов, которое описал Марко Бодрато^[6]. Алгоритм имеет пять основных шагов:

Разбиение
Вычисление в точках
Поточечное умножение
Интерполяция
Рекомпозиция

В типичной интерпретации больших чисел каждое целое представляется как последовательность цифр в позиционной системе счисления, где основанием счисления берётся некоторое (обычно большое) значение b. В нашем примере мы используем $b=10000$ , так что каждая цифра соответствует группе из четырёх десятичных цифр (в реализации в компьютере в качестве b обычно берётся степень двойки). Скажем, нужно перемножить два числа:

m	=	12	3456	7890	1234	5678	9012
n	=	9	8765	4321	9876	5432	1098.

Они много меньше, чем обычно обрабатываются алгоритмом Тоома — Кука, но они здесь иллюстрируют алгоритм.

Разбиение

Первым шагом является выбор основания счисления $B=b^{i}$ , так что число цифр как у числа m, так и у числа n по основанию B не превосходит k (например, 3 в Тоом-3). Обычно i задаётся выражением:

i=\max \left\{\left\lfloor {\frac {\left\lfloor \log _{b}m\right\rfloor }{k}}\right\rfloor ,\left\lfloor {\frac {\left\lfloor \log _{b}n\right\rfloor }{k}}\right\rfloor \right\}+1.

В нашем примере мы будем использовать Тоом-3, так что мы выбираем $B=b^{2}=10^{8}$ . Теперь мы разбиваем m и n по их основанию счисления B на цифры $m_{i},n_{i}$ :

{\begin{aligned}m_{2}&{}=123456\\m_{1}&{}=78901234\\m_{0}&{}=56789012\\n_{2}&{}=98765\\n_{1}&{}=43219876\\n_{0}&{}=54321098\end{aligned}}

Мы будем использовать эти цифры как коэффициенты в многочленах p и q степени (k − 1), со свойствами $p(B)=m$ и $q(B)=n$ :

p(x)=m_{2}x^{2}+m_{1}x+m_{0}=123456x^{2}+78901234x+56789012\,

q(x)=n_{2}x^{2}+n_{1}x+n_{0}=98765x^{2}+43219876x+54321098\,

После введения этих многочленов получаем, что если мы вычислим произведение $r(x)=p(x)q(x)$ , то ответом задачи будет $r(B)=m\times n$ .

В случае, когда перемножаемые числа имеют разные размеры, полезно использовать разные значения k для m и n, которые мы будем обозначать $k_{m}$ и $k_{n}$ . Например, алгоритм «Тоом-2.5» относится к алгоритму Тоома-Кука с $k_{m}=3$ и $k_{n}=2$ . В этом случае i в $B=b^{'}i$ обычно выбирается по выражению

i=\max \left\{\left\lfloor {\frac {\left\lceil \log _{b}m\right\rceil }{k_{m}}}\right\rfloor ,\left\lfloor {\frac {\left\lceil \log _{b}n\right\rceil }{k_{n}}}\right\rfloor \right\}.

Вычисление в точках

Подход Тоома — Кука вычисления произведения многочленов $p(x)q(x)$ обычен. Заметим, что многочлен степени d единственным образом определяется $d+1$ точками (например, прямая — это многочлен степени 1 и определяется двумя точками). Идеей является вычисление p(•) и q(•) в различных точках. Затем осуществляется произведение значений в этих точках, чтобы получить значения произведения многочленов. Наконец, интерполируем полученные значения для получения коэффициентов многочлена.

Поскольку $\deg(pq)=\deg(p)+\deg(q)$ , нам нужно $\deg(p)+\deg(q)+1=k_{m}+k_{n}-1$ точек для получения конечного результата. Назовём его d. В случае Тоом-3 $d=5$ . Алгоритм будет работать независимо от того, какие точки были выбраны (с несколькими небольшими исключениями — см. требование обратимости матрицы в разделе Интерполяция), но для упрощения алгоритма лучше использовать небольшие значения типа 0, 1, −1 и −2.

Одна необычная точка, которая часто используется, это бесконечность, то есть $\infty$ . Чтобы «вычислить» многочлен p на бесконечности, берут предел $p(x)/x^{\deg p}$ при стремлении x к бесконечности. Следовательно, $p(\infty )$ всегда равно значению старшего коэффициента (в примере выше — коэффициента $m_{2}$ ).

В нашем примере Тоом-3 мы используем точки 0, 1, −1, −2 и $\infty$ . Этот выбор упрощает вычисление в точках и даёт формулы:

{\begin{array}{lrlrl}p(0)&=&m_{0}+m_{1}(0)+m_{2}(0)^{2}&=&m_{0}\\p(1)&=&m_{0}+m_{1}(1)+m_{2}(1)^{2}&=&m_{0}+m_{1}+m_{2}\\p(-1)&=&m_{0}+m_{1}(-1)+m_{2}(-1)^{2}&=&m_{0}-m_{1}+m_{2}\\p(-2)&=&m_{0}+m_{1}(-2)+m_{2}(-2)^{2}&=&m_{0}-2m_{1}+4m_{2}\\p(\infty )&=&m_{2}&&\end{array}}

И аналогично для q. В нашем примере, значения, которые мы получаем, равны:

$p(0)$	=	$m_{0}$	=	56789012	=	56789012
$p(1)$	=	$m_{0}+m_{1}+m_{2}$	=	$56789012+78901234+123456$	=	135813702
$p(-1)$	=	$m_{0}-m_{1}+m_{2}$	=	$56789012-78901234+123456$	=	−21988766
$p(-2)$	=	$m_{0}-2m_{1}+4m_{2}$	=	$56789012-2\times 78901234+4\times 123456$	=	−100519632
$p(\infty )$	=	$m_{2}$	=	123456	=	123456
$q(0)$	=	$n_{0}$	=	54321098	=	54321098
$q(1)$	=	$n_{0}+n_{1}+n_{2}$	=	$54321098+43219876+98765$	=	97639739
$q(-1)$	=	$n_{0}-n_{1}+n_{2}$	=	$54321098-43219876+98765$	=	11199987
$q(-2)$	=	$n_{0}-2n_{1}+4n_{2}$	=	$54321098-2\times 43219876+4\times 98765$	=	−31723594
$q(\infty )$	=	$n_{2}$	=	98765	=	98765.

Как видно, эти значения могут быть отрицательными.

Для дальнейших объяснений полезно рассматривать этот процесс вычисления в точках как умножение матрицы на вектор справа, где каждая строка матрицы содержит степени одной из выбранных точек, а вектор содержит коэффициенты многочлена:

\left({\begin{matrix}p(0)\\p(1)\\p(-1)\\p(-2)\\p(\infty )\end{matrix}}\right)=\left({\begin{matrix}0^{0}&0^{1}&0^{2}\\1^{0}&1^{1}&1^{2}\\(-1)^{0}&(-1)^{1}&(-1)^{2}\\(-2)^{0}&(-2)^{1}&(-2)^{2}\\0&0&1\end{matrix}}\right)\left({\begin{matrix}m_{0}\\m_{1}\\m_{2}\end{matrix}}\right)=\left({\begin{matrix}1&0&0\\1&1&1\\1&-1&1\\1&-2&4\\0&0&1\end{matrix}}\right)\left({\begin{matrix}m_{0}\\m_{1}\\m_{2}\end{matrix}}\right).

Размерности матриц равны $d\times k_{m}$ для p и $d\times k_{n}$ для q. Строка для значения бесконечность имеет нулевые значения за исключением последнего столбца, в котором стоит 1.

Более быстрое вычисление значений

Вычисление значений в точках может быть выполнено быстрее, чем указано в приведённых выше формулах. Число элементарных операций (сложение/вычитание) может быть сокращено. Последовательность операций, придуманная Бодрато^[6] для Тоом-3, показана здесь для первого операнда (многочлена p):

$p_{0}$	$\leftarrow$	$m_{0}+m_{2}$	=	56789012 + 123456	=	56912468
$p(0)$	=	$m_{0}$	=	56789012	=	56789012
$p(1)$	=	$p_{0}+m_{1}$	=	56912468 + 78901234	=	135813702
$p(-1)$	=	$p_{0}-m_{1}$	=	56912468 − 78901234	=	−21988766
$p(-2)$	=	$(p(-1)+m_{2})\times 2-m_{0}$	=	$(-21988766+123456)\times 2-56789012$	=	− 100519632
$p(\infty )$	=	$m_{2}$	=	123456	=	123456.

Эта последовательность требует пять операций сложения/вычитания, на единицу меньше, чем при прямом вычислении. Более того, нам не нужно умножать на 4 при вычислении p(−2).

Поточечное умножение

В отличие от умножения многочленов p(•) и q(•), умножение вычисленных значений p(a) и q(a) просто сводится к умножению целых — более простого варианта исходной задачи. Мы рекурсивно используем нашу процедуру умножения для вычисления каждой пары значений в точках. В практических реализациях, когда операнды становятся меньше, алгоритм сводится к умножению в столбик^[англ.]. Если r — произведение многочленов, в нашем примере мы имеем:

$r(0)$	=	$p(0)q(0)$	=	$56789012\times 54321098$	=	3084841486175176
$r(1)$	=	$p(1)q(1)$	=	$135813702\times 97639739$	=	13260814415903778
$r(-1)$	=	$p(-1)q(-1)$	=	$-21988766\times 11199987$	=	−246273893346042
$r(-2$ )	=	$p(-2)q(-2)$	=	$-100519632\times -31723594$	=	3188843994597408
$r(\infty )$	=	$p(\infty )q(\infty )$	=	$123456\times 98765$	=	12193131840.

Как видим, эти значения могут быть отрицательными. Для достаточно больших чисел это наиболее дорогой шаг, единственный шаг, не зависящий линейно от размеров m и n.

Интерполяция

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

{\begin{aligned}\left({\begin{matrix}r(0)\\r(1)\\r(-1)\\r(-2)\\r(\infty )\end{matrix}}\right)&{}=\left({\begin{matrix}0^{0}&0^{1}&0^{2}&0^{3}&0^{4}\\1^{0}&1^{1}&1^{2}&1^{3}&1^{4}\\(-1)^{0}&(-1)^{1}&(-1)^{2}&(-1)^{3}&(-1)^{4}\\(-2)^{0}&(-2)^{1}&(-2)^{2}&(-2)^{3}&(-2)^{4}\\0&0&0&0&1\end{matrix}}\right)\left({\begin{matrix}r_{0}\\r_{1}\\r_{2}\\r_{3}\\r_{4}\end{matrix}}\right)\\&{}=\left({\begin{matrix}1&0&0&0&0\\1&1&1&1&1\\1&-1&1&-1&1\\1&-2&4&-8&16\\0&0&0&0&1\end{matrix}}\right)\left({\begin{matrix}r_{0}\\r_{1}\\r_{2}\\r_{3}\\r_{4}\end{matrix}}\right).\end{aligned}}

Эта матрица строится так же, как и на шаге вычисления значений многочлена в точках, за исключением, что матрица имеет размер d × d. Мы можем решить это уравнение с помощью метода Гаусса (метода исключений), но это будет очень дорого. Вместо этого мы используем факт, что используемые для вычисления значений точки выбраны специальным образом, позволяющим легко вычислить обратную матрицу (см. также Матрица Вандермонда), а тогда:

{\begin{aligned}\left({\begin{matrix}r_{0}\\r_{1}\\r_{2}\\r_{3}\\r_{4}\end{matrix}}\right)&{}=\left({\begin{matrix}1&0&0&0&0\\1&1&1&1&1\\1&-1&1&-1&1\\1&-2&4&-8&16\\0&0&0&0&1\end{matrix}}\right)^{-1}\left({\begin{matrix}r(0)\\r(1)\\r(-1)\\r(-2)\\r(\infty )\end{matrix}}\right)\\&{}=\left({\begin{matrix}1&0&0&0&0\\{\tfrac {1}{2}}&{\tfrac {1}{3}}&-1&{\tfrac {1}{6}}&-2\\-1&{\tfrac {1}{2}}&{\tfrac {1}{2}}&0&-1\\-{\tfrac {1}{2}}&{\tfrac {1}{6}}&{\tfrac {1}{2}}&-{\tfrac {1}{6}}&2\\0&0&0&0&1\end{matrix}}\right)\left({\begin{matrix}r(0)\\r(1)\\r(-1)\\r(-2)\\r(\infty )\end{matrix}}\right).\end{aligned}}

Теперь остаётся лишь вычислить произведение матрицы на вектор. Хотя матрица содержит дроби, результирующие коэффициенты будут целыми числами, так что вычисления можно вести в целочисленной арифметике путём сложения, вычитания и умножения/деления на маленькие константы. Здесь основной задачей является поиск эффективной последовательности операций, позволяющей вычислить произведение матрицы на вектор. Одну такую последовательность дал Бодрато^[6] для Тоом-3 и она для нашего примера следующая:

$r_{0}$	$\leftarrow$	$r(0)$	=	3084841486175176
$r_{4}$	$\leftarrow$	$r(\infty )$	=	12193131840
$r_{3}$	$\leftarrow$	$(r(-2)-r(1))/3$	=	(3188843994597408 − 13260814415903778)/3
			=	−3357323473768790
$r_{1}$	$\leftarrow$	$(r(1)-r(-1))/2$	=	$(13260814415903778-(-246273893346042))/2$
			=	6753544154624910
$r_{2}$	$\leftarrow$	$r(-1)-r(0)$	=	$-246273893346042-3084841486175176$
			=	−3331115379521218
$r_{3}$	$\leftarrow$	$(r_{2}-r_{3})/2+2r(\infty )$	=	$(-3331115379521218-(-3357323473768790))/2+2\times 12193131840$
			=	13128433387466
$r_{2}$	$\leftarrow$	$r_{2}+r_{1}-r_{4}$	=	−3331115379521218 + 6753544154624910 − 12193131840
			=	3422416581971852
$r_{1}$	$\leftarrow$	$r_{1}-r_{3}$	=	6753544154624910 − 13128433387466
			=	6740415721237444.

Мы теперь знаем произведение r наших многочленов:

{\begin{array}{rrr}r(x)=&{}&3084841486175176\\&+&6740415721237444x\\&+&3422416581971852x^{2}\\&+&13128433387466x^{3}\\&+&12193131840x^{4}\end{array}}

Если бы мы использовали другие $k_{m},k_{n}$ или выбрали другие точки для вычисления значений, матрица, а тогда и стратегия интерполяции, изменилась бы, но это не зависит от входных данных, а потому алгоритм может быть зашит «в железе» для любых данных параметров.

Рекомпозиция

Наконец, мы вычисляем r(B) для получения конечного результата. Это выполняется напрямую, поскольку B является степенью b, а потому умножение на степени B является сдвигом всего числа на целое число знаков основания b. В нашем примере $b=10^{4}$ и $B=b^{2}=10^{8}$ .

								3084	8414	8617	5176
						6740	4157	2123	7444
				3422	4165	8197	1852
		13	1284	3338	7466
+	121	9313	1840

	121	9326	3124	6761	1632	4937	6009	5208	5858	8617	5176

Результатом является произведение 1234567890123456789012 и 987654321987654321098.

Матрицы интерполяции для различных значений k

Здесь мы приведём матрицы интерполяции для нескольких различных малых значений $k_{m}$ и $k_{n}$ .

Тоом −1

Тоом-1 ( $k_{m}=k_{n}=1$ ) требует вычисления в одной точке, здесь выбирается 0. Алгоритм вырождается в умножение в столбик с единичной матрицей интерполяции:

\left({\begin{matrix}1\end{matrix}}\right)^{-1}=\left({\begin{matrix}1\end{matrix}}\right).

Тоом-1.5

Тоом-1.5 ( $k_{m}=2,k_{n}=1$ ) требует две точки, здесь выбраны 0 и $\infty$ . Матрица интерполяции тогда равна единичной матрице:

\left({\begin{matrix}1&0\\0&1\end{matrix}}\right)^{-1}=\left({\begin{matrix}1&0\\0&1\end{matrix}}\right).

Алгоритм также вырождается в умножение в столбик — оба коэффициента одного множителя умножается на один коэффициент другого множителя.

Тоом-2

Тоом-2 ( $k_{m}=2,k_{n}=2$ ) требует трёх точек, здесь выбраны 0, 1 и $\infty$ . Алгоритм получается тем же, что и алгоритм Карацубы с матрицей интерполяции:

\left({\begin{matrix}1&0&0\\1&1&1\\0&0&1\end{matrix}}\right)^{-1}=\left({\begin{matrix}1&0&0\\-1&1&-1\\0&0&1\end{matrix}}\right).

Тоом-2.5

Тоом-2.5 ( $k_{m}=3,k_{n}=2)$ требует 4 точки, здесь выбраны 0, 1, −1 и $\infty$ . Алгоритм тогда имеет матрицу интерполяции:

\left({\begin{matrix}1&0&0&0\\1&1&1&1\\1&-1&1&-1\\0&0&0&1\end{matrix}}\right)^{-1}=\left({\begin{matrix}1&0&0&0\\0&{\tfrac {1}{2}}&-{\tfrac {1}{2}}&-1\\-1&{\tfrac {1}{2}}&{\tfrac {1}{2}}&0\\0&0&0&1\end{matrix}}\right).

Примечания

↑ ¹ ² Knuth, 1997, с. 296.
↑ Crandall, Pomerance, 2005, с. 474.
↑ Crandall, Pomerance, 2005, с. 536.
↑ Knuth, 1997, с. 302.
↑ Positive Results Архивная копия от 6 января 2013 на Wayback Machine, chapter III of Stephen A. Cook: On the Minimum Computation Time of Functions.
↑ ¹ ² ³ Bodrato, 2007, с. 116–133.

Литература

D. Knuth. Section 4.3.3.A: Digital methods // The Art of Computer Programming. — Third Edition. — Addison-Wesley, 1997. — Т. 2. — С. 294.
Кнут Д.Э. Искусство программирования. Получисленные алгоритмы. — 2019. — Т. 2. — ISBN 978-5-907144-15-6.
R. Crandall, C. Pomerance. Section 9.5.1: Karatsuba and Toom–Cook methods // Prime Numbers – A Computational Perspective. — Second Edition. — Springer, 2005. — С. 473.
- Ричард Крэндалл, Карл Померанс. Простые числа. Криптографические и вычислительные аспекты. — Москва: Либроком, 2011. — ISBN 978-5-397-02060-2.
M. Bodrato. Toward Optimal Toom–Cook Multiplication for Univariate and Multivariate Polynomials in Characteristic 2 and 0 // Arithmetic of Finite Fields. — Springer, 2007. — Т. 4547. — (LNCS).

Ссылки

Toom-Cook 3-way multiplication from GMP documentation

[_46a6f8d240e854c8-1] ¹ ² Knuth, 1997, с. 296.

[_2ff7a89d8e8959bc-2] Crandall, Pomerance, 2005, с. 474.

[_2ff42c9d8e8650cb-3] Crandall, Pomerance, 2005, с. 536.

[_46a3ebd240e60922-4] Knuth, 1997, с. 302.

[COOK-5] Positive Results Архивная копия от 6 января 2013 на Wayback Machine, chapter III of Stephen A. Cook: On the Minimum Computation Time of Functions.

[_6625a9529a6c54f9-6] ¹ ² ³ Bodrato, 2007, с. 116–133.

[1]

[2]

[3]

[4]

[5]

[6]