В 1960-м году Андрей Колмогоров и несколько других советских пионеров информатики собрались на научном семинаре и выдвинули «гипотезу \(n^2\)»: невозможно перемножить два \(n\)-значных числа, быстрее, чем за \(O(n^2)\). Это подразумневает, что умножение «в столбик», придуманное шумерами как минимум четыре тысячи лет назад и никем на тот момент не побитое, является асимптотически оптимальным алгоритмом умножения двух чисел.
Через неделю 23-летний аспирант Анатолий Карацуба предложил метод умножения с оценкой времени работы \(O(n^{\log_2 3})\) и тем самым опроверг гипотезу.
Историческое примечание: эта задача сейчас решается за \(O(n \log n)\) с помощью алгоритма быстрого преобразования Фурье, который к концу 50-х годов уже изобрели, но использовали только «по назначению»: для обработки сигналов, а не для умножения чисел.
Любое число можно можно представить в виде многочлена, если вместо \(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 2 + a_2 \cdot 2^2 + \dots + a_n \cdot 2^n \end{aligned} \]
Наивная формула для умножения многочленов, если раскрыть все скобки, выглядит так:
\[ \left(\sum_{i=0}^n a_i x^i\right)\cdot\left(\sum_{j=0}^m b_j x^j\right)=\sum_{k=0}^{n+m}x^k\sum_{i+j=k}a_i b_j \]
Единственный нюанс: в случае чисел, когда мы хотим получить их реальное представление, нужно предварительно провести «каррирование»:
const int base = 10;
void carry(int *a, int n) {
int d = 0;
for (int i = 0; i < n; i++) {
a[i] += d;
d = a[i] % base;
a[i] /= base;
}
}Основание системы счисления может быть выбрано произвольно, но из соображений производительности его имеет смысл выбрать его как можно большим, избегая при этом переполнений. Например, при выборе основания \(10^6\) и использовании массива long long размера \(10^6\) (т. е. в данном случае перемножаются \((6 \cdot 10^6)\)-значные числа), ничего переполниться после одного умножения не должно: максимальное значение по формуле будет \((10^6)^3 = 10^{18}\), что помещается в 64 бит.
Основная схема алгоритма очень простая: в нём умножение двух чисел длины \(n\) небольшим алгебраическим трюком сводится к трём умножениям чисел длины \(\frac{n}{2}\). Примечательно, что нужные формулы были ещё у Чарльза Бэббиджа, работавшего в начале XIX века над механическим компьютером для подсчета налогов для британской короны, однако он не обратил внимания на возможность использования лишь трёх рекурсивных умножений вместо четырех.
То, почему он работает за такую странную асимптотику, весьма неочевидно, и для этого мы сначала докажем более мощную теорему, которая говорит об асимптотике большого класса алгоритмов вида «разделяй-и-властвуй», заменяющих исходную задачу на \(a\) задач размера \(\frac{n}{b}\).
Мастер-теорема. Пусть имеется рекуррента:
\[ T(n) = \begin{cases} a T(\frac{n}{b}) + \Theta(n^c), & n > n_0 \\ \Theta(1), & n \leq n_0 \end{cases} \]
Тогда:
Доказательство. Рассмотрим «дерево рекурсии» этого соотношения. В нём будет \(log_b n\) уровней. На \(k\)-том уровне будет \(a^k\) вершин, каждая из которых будет стоить \((\frac{n}{b^k})^c\) операций. Просуммируем значения во всех вершинах по всем уровням:
\[ T(n) = \sum_{k=0}^{\log_b n} a^k (\frac{n}{b^k})^c = n^c \sum_{k=0}^{\log_b n} (\frac{a}{b^c})^k \]
A. Если \(c > \log_b a\), то \(\sum (\frac{a}{b^с})^k\) это сумма убывающей геометрической прогрессии, которая не зависит от \(n\) и просто равна какой-то константе. Значит, \(T(n) = \Theta(n^c)\).
B. Если \(c = \log_b a\), то
\[ T(n) = n^c \sum_{k=0}^{\log_b n} (\frac{a}{b^c})^k = n^c \sum_{k=0}^{\log_b n} 1^k = \Theta(n^c \log_b n) \]
C. Если \(c < \log_b a\), то так как сумма прогрессии асимптотически эквивалентна своему старшему элементу,
\[ T(n) = n^c \sum_{k=0}^{\log_b n} (\frac{a}{b^c})^k = \Theta(n^c (\frac{a}{b^c})^{\log_b n}) = \Theta(n^c \cdot \frac{a^{\log_b n}}{n^c}) = \Theta(a^{\log_b n}) = \Theta(n^{\log_b a}) \]
Примечание. Для более точных оценок асимптотики «мерджа» теорема ничего не говорит. Например, если мердж занимает \(\Theta(n \log n)\) и задача разбивается каждый раз на две части, то асимптотика будет равна:
\[ \sum_{k=0}^{\log n} n \log \frac{n}{2^k} = \sum_{k=0}^{\log n} n (\log n - k) = n \sum_{k=0}^{\log n} k = \Theta (n \log^2 n) \]
В то же время эта рекуррента под условия теоремы не попадает. Можно лишь получить неточные границы \(\Omega (n \log n)\) и \(O(n^{1+\epsilon})\), если подставить \(c = 1\) и \(c = 1 + \epsilon\) соответственно. Заметим, что \(n \log n\) и \(n \log^2 n\) асимптотически меньше \(n^{1+\epsilon}\), каким бы маленьким \(\epsilon\) ни был.
Пусть у нас есть два многочлена \(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)) \]
Результат умножения исходных многочленов теперь можно посчитать по следующей формуле — внимание, алгебра:
\[ 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) \]
Корректность формулы можно проверить, просто выполнив нужные подстановки.
Обратим внимание, что результат умножения — многочлен размера \(2 n\).
Если посчитать необходимые операции, то выясняется, что для перемножения двух многочленов размера \(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\) за секунду,
В олимпиадных задачах длинное умножение само по себе встречается крайне редко. Гораздо чаще оно используется для подсчета каких-либо комбинаторных объектов через умножение многочленов.
Для примера мы решим задачу «Вор в магазине» с Educational Codeforces Round:
Дано \(n \leq 1000\) типов предметов разных положительных стоимостей. Каждая стоимость не превосходит \(1000\). Предметов каждого типа доступно бесконечное количество.
Требуется определить всевозможные суммы предметов, которые можно набрать, взяв ровно \(k\) предметов.
Составим многочлен степени \(1000\), в котором коэффициент при \(x_i\) равен единице, если в наборе существует предмет со стоимостью \(i\), и \(0\) в противном случае.
Если возвести этот многочлен в \(k\)-тую степень, то его можно записать так:
\[ \sum_k (x^k \sum_{t_1+\ldots+t_k=i} a_{t_1} a_{t_2} \ldots a_{t_k}) \]
Если посмотреть на эту формулу комбинаторно, то в получившемся многочлене коэффициент при \(x_i\) равен количеству способов набрать сумму \(i\), взяв ровно \(k\) предметов — а в задаче нас по сути интересует, равно ли это число нулю или нет.
Для возведения многочлена в \(k\)-тую степень можно использовать бинарное возведение в степень, а для умножения многочленов — алгоритм Карацубы.
const int maxn = (1<<20);
int a[maxn], res[maxn];
// мы не всегда хотим модифицировать исходные массивы при перемножении
// поэтому напишем обертку, которая работает, как *=
void mul(int *x, int *y, int n) {
int tx[n], ty[n];
memcpy(tx, x, sizeof tx);
memcpy(ty, y, sizeof ty);
memset(x, 0, sizeof tx);
karatsuba(tx, ty, x, n);
}
void binpow(int k) {
res[0] = 1;
int len = 1024;
while (k > 1) {
if (k & 1)
mul(res, a, len);
mul(a, a, len);
k /= 2;
len *= 2;
}
mul(res, a, len);
}
// осталось скормить функции binpow бинарный массив a,
// и res будет содержать a в k-той степениВ асимптотике будет учитёно только последнее (самое больше) умножение.
Решение в условиях задачи работает за 3 секунды с ограничением в 5. Его можно значительно ускорить, если вместо int-ов использовать байты или даже биты по аналогии с битсетом, потому что нас не интересует количество способов набрать какую-то сумму, а только равно ли это число нулю.
Похожий метод можно применить к матричному умножению — это называется алгоритмом Штрассена. В нём две матрицы разбиваются на \(8 = 4 + 4\) частей, перемножаются блочно, и сложной алгеброй от одного из 8 умножений получается избавиться, что даёт асимптотику \(O(n^{\log_7 8}) \approx O(n^{2.81})\).
И в алгоритме Штрассена, и в алгоритме Карацубы можно достичь и лучшей асимптотики, если разбивать объекты на большее число частей. Однако, в реальности это не применяется, потому что в асимптотиках подобных улгоритмов скрыта непрактично большая константа.