Алгоритм Карацубы - Алгоритмика

В 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})$ .

И в алгоритме Штрассена, и в алгоритме Карацубы можно достичь и лучшей асимптотики, если разбивать объекты на большее число частей. Однако, в реальности это не применяется, потому что в асимптотиках подобных алгоритмов скрыта непрактично большая константа.