Метод умножения Шёнхаге — Штрассена
Содержание
Теория
Свертка
Допустим, что мы перемножаем числа 123 и 456 «в столбик», однако без выполнения переноса. Результат будет выглядеть так:
| 1 | 2 | 3 | ||
| × | 4 | 5 | 6 | |
| 6 | 12 | 18 | ||
| 5 | 10 | 15 | ||
| 4 | 8 | 12 | ||
| 4 | 13 | 28 | 27 | 18 |
Эта последовательность (4, 13, 28, 27, 18) называется ациклической или линейной свёрткой от последовательностей (1,2,3) и (4,5,6). Зная ациклическую свёртку двух последовательностей, рассчитать произведение несложно: достаточно выполнить перенос (например, в самом правом столбце, мы оставляем 8 и добавляем 1 к столбцу, содержащему 27). В нашем примере это приводит к результату 56088.
Есть и другие типы свёрток, которые могут быть полезны. Допустим, что входящие последовательности содержат n элементов (в примере — 3). Тогда результирующая линейная свёртка содержит n + n − 1 элементов; если мы возьмём самый правый столбец n элементов и добавим самый левый стоблец n − 1 ', в результате мы получим циклическую свёртку:
| 28 | 27 | 18 | ||
| + | 4 | 13 | ||
| 28 | 31 | 31 | ||
Если мы произведём перенос при циклическом свёртывании, результат будет тот же, что и при произведении чисел по модулю Bn − 1 (в данном примере это 103 − 1 = 999). Выполним перенос по (28, 31, 31) и получим 3141, при этом 3141 ≡ 56088 (mod 999).
Наоборот, если мы возьмём самый правый столбец n элементов и вычтем самый левый столбец n−1 элементов, то в результате мы получим обратную свёртку:
| 28 | 27 | 18 | ||
| − | 4 | 13 | ||
| 28 | 23 | 5 | ||
Если мы произведём перенос при обратном свёртывании, то результат будет тот же, что и при произведении чисел по модулю Bn + 1. В данном примере, 103 + 1 = 1001, выполним перенос по (28, 23, 5) и получим 3035, при этом 3035 ≡ 56088 (mod 1001). Обратная свёртка может содержать отрицательные числа, которые могут быть убраны во время переноса, используя ту же технику, что и при длинных вычитаниях.
Теорема о свёртке
Как и другие методы, основанные на быстром преобразовании Фурье, алгоритм Фюрера в корне зависит от теоремы о свёртке, которая обеспечивает эффективный способ посчитать циклическую свёртку двух последовательностей. Её идея состоит в следующем:
- Циклическая свёртка двух векторов может быть найдена через дискретное преобразование Фурье (ДПФ) каждого из них, путём произведения результирующих векторов элемент за элементом, с последующим обратным преобразованием Фурье (ОДПФ).
Или через формулы:
- CyclicConvolution(X, Y) = IDFT(DFT(X) · DFT(Y)), где:
- CyclicConvolution — циклическая свертка,
- DFT — дискретное преобразование Фурье,
- IDFT — обратное дискретное преобразование Фурье.
Если мы посчитаем ДПФ и ОДПФ, используя быстрое преобразование Фурье, и вызовем наш алгоритм перемножения рекурсивно, чтобы перемножить входы(?) преобразованных векторов DFT(X) и DFT(Y), то в результате мы получим эффективный алгоритм для расчёта циклической свёртки.
В этом алгоритме, гораздо эффективней считать обратную циклическую свёртку; как оказывается, немного модифицированная версия теоремы о свёртке может позволить и это. Предположим, что вектора X и Y имеют длину n, и a — примитивный корень порядка 2n (это означает, что a2n = 1 и все меньшие степени a не равны 1). Таким образом мы можем определить третий вектор A, называемый вектор веса, обладающий следующими свойствами:
- A = (aj), 0 ≤ j < n
- A−1 = (a−j), 0 ≤ j < n
Теперь мы можем записать:
- NegacyclicConvolution(X, Y) = A−1 · IDFT(DFT(A · X) · DFT(A · Y)), где
- NegacyclicConvolution — Обратная циклическая сертка,
- DFT — дискретное преобразование Фурье,
- IDFT — обратное дискретное преобразование Фурье.
Другими словами, это то же самое за исключением того, что входящие векторы умножены на A, а результат умножен на A−1.
Выбор кольца
Дискретное преобразование Фурье — абстрактная операция, которая может быть выполнена в любом алгебраическом кольце; обычно оно берётся из поля комплексных чисел, но фактически использовать комплексную арифметику с достаточной точностью, чтобы обеспечить точные результаты, медленно и неэффективно. Вместо этого мы можем использовать теоретико-числовое преобразование, которое производит преобразование в поле целых чисел по модулю N для некоторого целого N.
Так же как есть примитивные корни единицы любого порядка на комплексной плоскости, при любом заданном n мы можем выбрать подходящее N такое, что b — примитивный корень единицы порядка n в поле целых чисел по модулю N (другими словами, bn ≡ 1 (mod N), и все меньшие степени b не равны 1 mod N).
Алгоритм тратит большую часть времени на рекурсивное выполнение произведения меньших чисел; в простом варианте алгоритма это происходит в ряде мест:
- Внутри алгоритма быстрого преобразования Фурье, примитивный корень единицы b неоднократно возводится в степень и умножается на другие числа.
- При возведении в степень примитивного корня единицы a для получения вектора веса A с последующим умножением векторов A или A−1 на другие вектора.
- При выполнении последовательного перемножения преобразованных векторов.
Ключевой момент — выбрать N, модуль, равный 2n + 1 для некоторого целого n. У этого способа есть ряд преимуществ в ряде стандартных систем, в которых большие целые числа представлены в двоичном виде:
- Любое число может быть быстро уменьшено по модулю 2n + 1 используя только сдвиг и сложение.
- Любые примитивные корни единицы в этом кольце могут быть записаны в форме 2k; соответственно мы можем умножать или делить любое число на корень из единицы используя сдвиг.
- Поэлементное рекурсивное перемножение преобразованный векторов может быть выполнено, используя обратную свёртку, которая работает быстрее, чем ациклическая свёртка, и в которой уже есть уменьшение результата
по модулю 2n + 1.
Пример
Положим задано два числа: X и Y, каждое длиной 512 бит. Требуется найти значение RES = (X*Y) mod 2N+1
- Этап 1. Основные параметры алгоритма:
- len1 = 512
- len2 = 512
- N = 1024
- 2k ~ sqrt(N) = 32, k = 5
- n >= 2*N/2k+k = (2048/32)+5 = 69. Поскольку n должно делиться на 2k, то можно выбрать n = 128.
- Этап 2 (возникает на умножении DFT(X)*DFT(Y) mod 2128+1):
- len1 = 129
- len2 = 129
- N = 128
- 2k ~ (sqrt(N) = 11.3). Выберем k = 4.
- n >= 2*N/2k+k = (256/16)+4 = 20. Поскольку n должно делиться на 2k, то можно выбрать n = 32.
- Этап 3 (возникает на умножении DFT(X)*DFT(Y) mod 232+1):
- len1 = 33
- len2 = 33
это умножение можно делать аппаратно, не углубляюсь в рекурсию.
Описание
- Документация GMP LIB
- Статья A GMP-BASED IMPLEMENTATION OF SCHONHAGE-STRASSEN’S LARGE INTEGER MULTIPLICATION ALGORITHM
- Описание в английской Wikipedia
