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

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

${\displaystyle 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, так что мы выбираем ${\displaystyle B=b^{2}=10^{8}}$. Теперь мы разбиваем m и n по их основанию счисления B на цифры ${\displaystyle m_{i},n_{i}}$:

${\displaystyle {\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), со свойствами ${\displaystyle p(B)=m}$ и ${\displaystyle q(B)=n}$:

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

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

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

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

${\displaystyle 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\}.}$

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

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

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

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

${\displaystyle {\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. В нашем примере, значения, которые мы получаем, равны:

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

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

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

${\displaystyle \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).}$

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

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

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

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

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

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

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

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

${\displaystyle {\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. Мы можем решить это уравнение с помощью метода Гаусса (метода исключений), но это будет очень дорого. Вместо этого мы используем факт, что используемые для вычисления значений точки выбраны специальным образом, позволяющим легко вычислить обратную матрицу (см. также Матрица Вандермонда), а тогда

${\displaystyle {\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 и она для нашего примера следующая:

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

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

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

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

Наконец, мы вычисляем r(B) для получения конечного результата. Это выполняется напрямую, поскольку B является степенью b, а потому умножение на степени B является сдвигом всего числа на целое число знаков основания b. В нашем примере ${\displaystyle b=10^{4}}$ и ${\displaystyle 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 иd 987654321987654321098.

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

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

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

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

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

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

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

${\displaystyle \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).}$