В 1960-м году Андрей Колмогоров вместе с другими пионерами советской информатики собрались на научном семинаре и выдвинули «гипотезу $n^2$»: невозможно перемножить два $n$-значных числа, быстрее, чем за $O(n^2)$. Это подразумевает, что умножение «в столбик», придуманное шумерами как минимум четыре тысячи лет назад и никем на тот момент не побитое, является асимптотически оптимальным алгоритмом умножения двух чисел.
Через неделю 23-летний аспирант Анатолий Карацуба предложил метод умножения с оценкой времени работы $O(n^{\log_2 3})$ и тем самым опроверг гипотезу. В этой статье мы и рассмотрим этот метод.
#Алгоритм
Основная идея алгоритма очень простая: в нём умножение двух чисел длины $n$ небольшим алгебраическим трюком сводится к трём умножениям чисел длины $\frac{n}{2}$. Число элементов на каждом уровне рекурсии будет расти, но на самом нижнем будет суммарно всего $O(n^{\log_2 3})$ элементов, чем и объясняется такая странная асимптотика.
В нашей версии алгоритма мы будем перемножать не числа, а многочлены. Это эквивалентная задача — если вместо $x$ подставить основание системы счисления, то в качестве коэффициентов можно взять последовательность цифр числа:
$$ \begin{aligned} A(x) &= a_0 + a_1\cdot x + a_2 \cdot x^2 + \dots + a_n \cdot x^n \\ &= a_0 + a_1\cdot 10 + a_2 \cdot 10^2 + \dots + a_n \cdot 10^n \end{aligned} $$ Итак, пусть у нас есть два многочлена $a(x)$ и $b(x)$ равной длины $n = 2k$ и мы хотим их перемножить. Разделим их коэффициенты на две равные части и представим как $$ a(x) = a_1 (x) + x^k a_2(x) \\b(x) = b_1 (x) + x^k b_2(x) $$ Теперь рекурсивно вычислим многочлены-произведения $p_1$ и $p_2$: $$ p_1(x) = a_1(x) \cdot b_1(x) \\ p_2(x) = a_2(x) \cdot b_2(x) $$ А также многочлен $t$: $$ t(x) = ( a_1(x) + a_2(x) ) \cdot (b_1(x) + b_2(x)) $$ Результат умножения исходных многочленов (многочлен размера $2n$) теперь можно посчитать по следующей формуле — внимание, алгебра: $$ c(x) = a(x) \cdot b(x) = p_1(x) + x^k \cdot (t(x) - p_1(x) - p_2(x)) + x^{2k} \cdot p_2(x) $$Корректность формулы можно проверить, просто выполнив нужные подстановки.
#Анализ
Если посчитать необходимые операции, то выясняется, что для перемножения двух многочленов размера $n$ нам нужно посчитать три произведения — $p_1$, $p_2$ и $t$ — размера $\frac{n}{2}$ и выполнить константное количество сложений, вычитаний и сдвигов (домножений на $x^k$), которые суммарно можно выполнить за $O(n)$.
Мастер-теорема утверждает, что в данном случае асимптотика всего алгоритма будет $\Theta (n^{\log_2 3}) \approx \Theta (n^{1.58})$: наша задача разбивается на $a = 3$ части в $b = 2$ раз меньшего размера, а объединение происходит за $O(n)$.
#Реализация
Для простоты будем предполагать, что $n$ это степень двойки. Если это не так, то в зависимости от обстоятельств это можно исправить одним из двух костылей:
- Можно дополнить коэффициенты многочлена нулями до ближайшей степени двойки — в худшем случае это будет работать в $2^{1.58} \approx 3$ раза дольше.
- Можно «отщепить» последний коэффициент от многочленов и свести задачу размера $(2k + 1)$ к задаче размера $2k$ и константному количество сложений.
Мы будем использовать первый метод, так как он проще в реализации.
Основные соображения по поводу эффективной реализации:
- Нужно выделять как можно меньше лишней памяти, для чего нужно переиспользовать имеющиеся массивы.
- Все арифметические операции нужно реализовать как простые линейные проходы по массивам, чтобы компилятор смог их векторизовать.
- Вместо использования базы вида
if (n == 1) c[0] = a[0] * b[0], имеет смысл, начиная с какого-то размера задачи, использовать более эффективное наивное умножение за квадрат.
Сначала приведем код основной рекурсивной процедуры, а потом объясним, почему он работает:
void karatsuba(int *a, int *b, int *c, int n) {
if (n <= 64) {
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
c[i + j] += a[i] * b[j];
} else {
int k = n / 2;
int l[k], r[k], t[n] = {0};
for (int i = 0; i < k; i++) {
l[i] = a[i] + a[k + i];
r[i] = b[i] + b[k + i];
}
karatsuba(l, r, t, k); // считает t
karatsuba(a, b, c, k); // считает p1
karatsuba(a + k, b + k, c + n, k); // считает p2
int *t1 = t, *t2 = t + k;
int *s1 = c, *s2 = c + k, *s3 = c + 2 * k, *s4 = c + 3 * k;
for (int i = 0; i < k; i++) {
int c1 = s2[i] + t1[i] - s1[i] - s3[i];
int c2 = s3[i] + t2[i] - s2[i] - s4[i];
c[k + i] = c1;
c[n + i] = c2;
}
}
}
После трёх рекурсивных вызовов массив $c$ — это конкатенация $p_1$ и $p_2$.
После этого, для подсчета самого многочлена $c$ проще всего мысленно разделить его на четыре равные части, а многочлен $t$ — на две половины $t_1$ и $t_2$, а затем посмотреть на формулу и подумать, как изменится каждая часть:
- $s_1$: не меняется — это первая половина $p_1$.
- $s_2$: выражается как $s_2 + t_1 - s_1 - s_3$, то есть изменяется на «первую» половину $t - p_1 - p_2$.
- $s_3$: выражается как $s_3 + t_2 - s_2 - s_4$, то есть изменяется на «вторую» половину $t - p_1 - p_2$.
- $s_4$: не меняется — это вторая половина $p_2$.
Из-за векторизации важно использовать максимально «лёгкий» тип данных и при возможности компилировать с AVX:
#pragma GCC optimize("O3")
#pragma GCC target("avx2")
Реализация достаточно эффективна: она может перемножить два многочлена размера $4 \cdot 10^5$ за секунду.
#Обобщение идеи
Похожий метод можно применить к матричному умножению — это называется алгоритмом Штрассена. В нём две матрицы разбиваются на $8 = 4 + 4$ частей, перемножаются блочно, и сложной алгеброй от одного из 8 умножений получается избавиться, что даёт асимптотику $O(n^{\log_2 7}) \approx O(n^{2.81})$.
И в алгоритме Штрассена, и в алгоритме Карацубы можно достичь и лучшей асимптотики, если разбивать объекты на большее число частей. Однако, в реальности это не применяется, потому что в асимптотиках подобных алгоритмов скрыта непрактично большая константа.