Лекция 16

Алгоритм Шёнхаге-Штрассена

(Конспект: Т. Лазовская)

Сложение n-битовых чисел можно осуществить за O(n) (во-первых, у нас имеется такая инструкция RAM-машины; во-вторых, ясно, что это действительно можно реализовать физически). Рассмотрим умножение n-битовых чисел (пусть n — степень двойки; можно считать так в любом случае, так как иначе время работы увеличится не более, чем в константу раз). Умножение �в столбик� даст, очевидно, время Ω(n^{2}). Хотелось бы построить более быстрый алгоритм.

Почему мы не включили операцию умножения в RAM-машину? Именно потому, что непонятно, как физически реализовать ее так, чтобы она работала время O(n). Мы могли бы включить другие инструкции (и включим их, так как они нам понадобятся): умножение и деление на степень двойки, т.е. сдвиг двоичного представления:

LSHIFT a,
RSHIFT a

(содержимое нулевого регистра сдвигается влево или вправо на a битов, т.е. умножается или делится на 2^{a}; RSHIFT — обобщение операции HALF). Ясно, что эту инструкцию совсем просто реализовать.

16.1 Простой алгоритм

Начнем с простого алгоритма, время работы которого составляет O(n^{log_{2} 3}).

Имеем два n-битовых числа a и b. Разделим их в битовом представлении на n/2-битовые a_{1}, a_{2} и b_{1}, b_{2} (соответственно), а затем перемножим эти числа рекурсивно:

a = a_{1} · 2^{n/2} + a_{2},
b = b_{1} · 2^{n/2} + b_{2},
a · b = a_{1}b_{1} · 2^{n} + (a_{1}b_{2} + a_{2}b_{1}) · 2^{n/2} + a_{2}b_{2};

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

a_{1}b_{2} + a_{2}b_{1} = (a_{1} + a_{2})(b_{1} + b_{2}) − a_{1}b_{1} − a_{2}b_{2}.

Тем самым получаем, что нам достаточно трех умножений n/2-битовых чисел (ибо умножение на степень двойки — по существу, не умножение), т.е. рекуррентное уравнение для времени работы —

T (n) ⩽ 3T (n/2) + cn.

Тогда (в первой части курса была соответствующая теорема) T (n) = O(n^{log_{2} 3}).

Замечание 16.1. Нюанс: строго говоря, мы перемножали не n/2-битные, а (n/2 + 1)-битные числа x = a_{1} + a_{2} и y = b_{1} + b_{2}, но одно к другому сводится за линейное время. Действительно, пусть x = 2A + x′, y = 2B + y′, где A и B — n/2-битные числа, x′ и y′ — биты. Тогда xy = 4AB + 2Ay′ + 2Bx′ + x′y′. �Сложным� умножением здесь является только A · B; остальные �умножения� реализуются за линейное время, поскольку это умножения на степени двойки или 0.

O(n^{log_{2} 3}) — лучше, чем O(n^{2}), но все еще много. Поэтому будем строить другой алгоритм.

16.2 Дискретное преобразование Фурье (ДПФ)

Пусть R — кольцо (в нем имеются 0, 1, +, ·). Будем работать с векторами размерности n над R. Пусть в R имеется w — первообразный корень степени n из 1, и пусть существует n^{−1}.

Определение 16.1. Для вектора a = (a_{0}, a_{1}, ..., a_{n−1})^{t} (в виде столбца) определим преобразование Фурье от него:

F (a) = Aa,

где A_{ij}=w^{ij}. Обратное преобразование Фурье — преобразование с матрицей (A^{−1})_{ij} = n^{−1}w^{−ij}.

Задача 16.1. Показать, что оно действительно обратное (A^{−1}A = E).

Использовать его мы будем примерно так. У нас векторы a и b задают числа (а можно считать — что многочлены). Сделаем ДПФ: a ↦ F (a), b ↦ F (b). После этого перемножим поэлементно F (a) и F (b). Получим что-то типа произведения значений многочленов в соответствующих точках; затем применим обратное ДПФ и получим искомый результат ab.

Сформулируем строгие утверждения.

Определение 16.2. Пусть a = (a_{0}, …, a_{n−1})^{t}, b = (b_{0}, ..., b_{n−1})^{t}. Определим новую операцию свертка (⊙); a ⊙ b — это 2n-мерный вектор с элементами

(a ⊙ b)_{i} = ∑_{j=0}^{n−1} a_{j}b_{i−j}

для 0 ⩽ i ⩽ 2n − 1; здесь и далее элементы, которые не определены (например, b_{−1}) считаются равными 0.

Замечание 16.2. Коэффициенты свертки равны коэффициентам произведения многочленов ∑_{i} a_{i}x^{i} и ∑_{i} b_{i}x^{i}.

Лемма 16.1. Дополним a и b до 2n элементов: a′′ = (a_{0}, …, a_{n−1}, 0, …, 0)^{t}, b′′ = (b_{0}, …, b_{n−1}, 0, …, 0)^{t}. Теперь применим ДПФ к каждому из них и перемножим результаты покомпонентно. Затем применим обратное ДПФ. Так мы получим вектор длины 2n, равный свертке a ⊙ b.

Упражнение 16.1. Доказать лемму 16.1.

Определение 16.3 (отрицательно обернутая свертка). Это n-мерный вектор a ⊖ b с элементами (a ⊖ b)_{i} = ∑_{j=0}^{i} a_{j}b_{i−j} − ∑_{j=i+1}^{n−1} a_{j}b_{n+i−j}.

Лемма 16.2. Пусть у нас есть не только w, но и ψ — корень из него (ψ^{2} = w). Рассмотрим вектора

a′ = (a_{0}, a_{1}ψ, …, a_{n−1}ψ^{n−1}),
b′ = (b_{0}, b_{1}ψ, …, b_{n−1}ψ^{n−1}),
d′ = (d_{0}, d_{1}ψ, …, d_{n−1}ψ^{n−1}),

где d_{i} — коэффициенты отрицательно обернутой свертки (d = a ⊖ b). Утверждение: d′ = F^{−1}(F (a′) ∗ F (b′)), где ∗ — покомпонентное произведение двух векторов.

Доказательство. Достаточно показать, что F (a′) ∗ F (b′) = F (d′).

(F (a′))_{i} = ∑_{j=0}^{n−1} w^{ij}ψ^{j}a_{j},

(F (b′))_{i} = ∑_{j′=0}^{n−1} w^{ij′}ψ^{j′}b_{j′},

(F (d′))_{i} = ∑_{j′′=0}^{n−1} w^{ij′′}ψ^{j′′} (∑_{k=0}^{j′′} a_{k}b_{j′′−k} − ∑_{k=j′′+1}^{n−1} a_{k}b_{n+j′′−k}) .

Теперь перемножим:

(F (a′))_{i} ∗ (F (b′))_{i} = ∑_{j=0}^{n−1} ∑_{j′=0}^{n−1} a_{j}b_{j′}ψ^{j+j′}w^{i(j+j′)}.

Слагаемые, соответствующие j + j′ < n, имеются и в F (d′)_{i} (при этом слагаемое с ψ^{j+j′} соответствует слагаемому с ψ^{j′′}). Заметим, что ψ^{n} = −1, w^{n} = 1. Поэтому слагаемые, соответствующие j + j′ ⩾ n, могут быть преобразованы так:

a_{j}b_{j′}ψ^{j+j′}w^{i(j+j′)} = −a_{j}b_{j′}ψ^{j+j′−n}w^{i(j+j′−n)},

но это как раз оставшиеся слагаемые из F (d′)_{i} (теперь слагаемое с ψ^{j+j′−n} соответствует слагаемому с ψ^{j′′}).

16.3 Алгоритм Шёнхаге-Штрассенa^{1}

Достаточно научиться перемножать n-битные числа по mod (2^{n} + 1). (Чтобы перемножить точно, можно представить n-битные числа как 2nбитные.) Не умаляя общности, можно считать, что n = 2^{k}.

Требуется перемножить два n-битных числа, u и v. Разделим их двоичные представления на блоки длиной l битов каждый, где n = lb, b = 2^{k/2} или b = 2^{(k−1)/2} (в зависимости от четности k). Очевидно, b — количество блоков.

u = (u_{0}, u_{1}, …, u_{b−1})^{t},
v = (v_{0}, v_{1}, …, v_{b−1})^{t},
w = (w_{0}, w_{1}, …, w_{b−1})^{t},
w = uv mod (2^{n} + 1) (как числа).

(Если задача тривиальна, т.е. один из сомножителей — 2^{k}, этот алгоритм не используют.) Нам надо найти w_{i}. Рассмотрим uv, еще не взятое по модулю 2^{n} + 1; оно состоит из блоков

y_{i} = ∑_{j=0}^{b−1} u_{j}v_{i−j};

именно,

uv = y_{0} + y_{1} · 2^{l} + … + y_{b−1} · 2^{l(b−1)} + y_{b} · 2^{lb} + … + y_{b+i} · 2^{l(b+i)} + …

Поскольку 2^{lb} ≡ −1 (mod 2^{n} + 1), получаем w_{i} = y_{i} − y_{b+i}, т.е. w_{i} — коэффициенты отрицательно обернутой свертки.

Мы знаем для i < b, что 0 ⩽ y_{i} ⩽ (1 + i)2^{2l}, т.к. в y_{i} = ∑_{j=0}^{b−1} u_{j}v_{i−j} всего 1 + i ненулевых членов. Также 0 ⩽ y_{i} ⩽ (b − 1 − i)2^{2l} для i ⩾ b. Cледовательно, w_{i} достаточно вычислить по модулю ⩾ b2^{2l}.

Пусть w′_{i} = w_{i} mod b, w′′_{i} = w_{i} mod (2^{2l} +1) (одно из них мы вычислим с помощью преобразования Фурье, другое — с помощью первого, �простого�, алгоритма). Тогда можно найти и w_{i} = (2^{2l} + 1)[(w′_{i} − w′′_{i})mod b] + w′′_{i}.

1A. Schönhage; V. Strassen

Найдем w′′_{i}. Рассмотрим кольцо остатков по mod (2^{2l} + 1). Пусть ψ = 2^{2l/b}. Убедимся, что он — корень из первообразного корня степени b: ψ^{b} = 2^{2l} ≡ −1 (mod (2^{2l} + 1)). Теперь, чтобы найти w′′_{i}, можно воспользоваться леммой 16.2. Для этого надо уметь делить и умножать на степени ψ, но это легко, поскольку ψ — степень двойки. Кроме того, надо быстро уметь находить преобразование Фурье и обратное к нему (перемножать матрицу на вектор — слишком долго); это мы научимся делать на следующей лекции. И, наконец, надо поэлементно перемножить два вектора (F (u′) и F (v′)), т.е. уметь умножать 2l-битные числа; для этого мы рекурсивно воспользуемся нашим алгоритмом.

Найдем w′_{i}. Воспользуемся первым (�простым�) алгоритмом, но — для всех w′_{i} сразу. Рассмотрим два числa из 3b log b битов:

|0 0 (u_{0} mod b)|0 0 (u_{1} mod b)| … |,
|0 0 (v_{0} mod b)|0 0 (v_{1} mod b)| … | —

дополнили наши числа нулями таким образом, чтобы получилось b блоков по 3 log b битов. Перемножим эти длинные числа �простым� алгоритмом и получим

|y_{0} mod b|y_{1} mod b| … |y_{b−1} mod b| —

нули мы записывали так, чтобы после перемножения блоки не пересекались. Осталось попарно вычесть полученные числа: w′_{i} = (y_{i} mod b − y_{b+i} mod b)mod b.

Замечание 16.3. Сложение по модулю степени двойки делается просто.

Итак, нам осталось научиться быстро вычислять ДПФ (и обратное к нему) и оценить время работы всего алгоритма в целом.