Алгоритм Шёнхаг-Strassen является асимптотически быстрого алгоритма умножения больших чисел . Он был разработан Арнольдом Шёнхаге и Volker Strassen в 1971. [1] Во время выполнения битовой сложности есть в Big O нотации ,для двух n- значных чисел. Алгоритм использует рекурсивные быстрые преобразования Фурье в кольцах с 2 n +1 элементами, особый тип теоретико-числового преобразования .
Алгоритм Шёнхаге – Штрассена был асимптотически самым быстрым методом умножения, известным с 1971 по 2007 год, когда был объявлен новый метод, алгоритм Фюрера , с более низкой асимптотической сложностью; [2] однако алгоритм Фюрера в настоящее время дает преимущество только для астрономически больших значений и не используется на практике (см. Галактические алгоритмы ).
На практике алгоритм Шёнхаге-Штрассена начинает превосходить более старые методы, такие как умножение Карацубы и Тоома-Кука, для чисел от 2 2 15 до 2 2 17 (от 10 000 до 40 000 десятичных цифр). [3] [4] [5] Библиотека GNU Multi-Precision использует его для значений от 1728 до 7808 64-битных слов (от 33 000 до 150 000 десятичных цифр), в зависимости от архитектуры. [6] Существует Java-реализация Шёнхаге-Штрассена, в которой число десятичных знаков превышает 74 000. [7]
Применения алгоритма Шёнхаге-Штрассена включают математический эмпиризм , такой как Великий поиск прайма Мерсенна в Интернете и вычисление приближений числа π , а также практические приложения, такие как подстановка Кронекера , в которой умножение многочленов с целыми коэффициентами может быть эффективно сведено к большому целому числу. умножение; на практике это используется GMP-ECM для факторизации эллиптической кривой Ленстры . [8]
Подробности
В этом разделе подробно объясняется, как реализуется Schönhage – Strassen. Он основан прежде всего на обзоре метода Крэндалла и Померанса в их простых числах: вычислительная перспектива . [9] Этот вариант несколько отличается от оригинального метода Шёнхаге тем, что он использует дискретное взвешенное преобразование для более эффективного выполнения неациклических сверток . Еще один источник подробной информации - это книга Кнута « Искусство компьютерного программирования» . [10] Раздел начинается с обсуждения строительных блоков и заканчивается пошаговым описанием самого алгоритма.
Свертки
Предположим, мы умножаем два числа, например, 123 и 456, используя длинное умножение на цифры с основанием B , но без какого-либо переноса. Результат может выглядеть примерно так:
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 год |
Если мы выполняем циклическую свертку, результат эквивалентен произведению входов по модулю B n - 1. В этом примере 10 3 - 1 = 999, выполнение продолжения (28, 31, 31) дает 3141, и 3141 ≡ 56088 (мод.999).
И наоборот, если мы возьмем крайние правый п элементы и вычесть левый п -1 элементы, это производит negacyclic свертки :
28 год | 27 | 18 | ||
- | 4 | 13 | ||
28 год | 23 | 5 |
Если мы выполним неациклическую свертку, результат будет эквивалентен произведению входов по модулю B n + 1. В этом примере 10 3 + 1 = 1001, выполнение продолжения (28, 23, 5) дает 3035, и 3035 ≡ 56088 (мод 1001). Неациклическая свертка может содержать отрицательные числа, которые можно исключить при переносе с помощью заимствования, как это делается при длинном вычитании.
Теорема свертки
Как и другие методы умножения, основанные на быстром преобразовании Фурье , Шёнхаге – Штрассен в основном зависит от теоремы о свертке , которая обеспечивает эффективный способ вычисления циклической свертки двух последовательностей. В нем говорится, что:
- Циклическую свертку двух векторов можно найти, взяв дискретное преобразование Фурье (ДПФ) каждого из них, умножив результирующие векторы элемент на элемент, а затем используя обратное дискретное преобразование Фурье (IDFT).
Или символами:
- Циклическая свертка ( X, Y ) = IDFT (DFT ( X ) · DFT ( Y ))
Если мы вычисляем DFT и IDFT с использованием алгоритма быстрого преобразования Фурье и рекурсивно вызываем наш алгоритм умножения для умножения элементов преобразованных векторов DFT ( X ) и DFT ( Y ), это дает эффективный алгоритм для вычисления циклической свертки.
В этом алгоритме будет более полезно вычислить неациклическую свертку; как оказалось, слегка измененная версия теоремы о свертке (см. дискретное взвешенное преобразование ) также может позволить это. Предположим , что векторы X и Y имеют длину п , и является примитивным корнем из единицы из порядка 2 п (то есть, 2 п = 1 и для всех меньших полномочий не равен 1). Затем мы можем определить третий вектор A , называемый вектором весов , как:
- A = ( a j ), 0 ≤ j < n
- A −1 = ( a - j ), 0 ≤ j < n
Теперь мы можем констатировать:
- NegacyclicConvolution ( X , Y ) = A −1 · IDFT (DFT ( A · X ) · DFT ( A · Y ))
Другими словами, это то же самое, что и раньше, за исключением того, что входные данные сначала умножаются на A , а результат умножается на A −1 .
Выбор кольца
Дискретное преобразование Фурье - это абстрактная операция, которая может выполняться в любом алгебраическом кольце ; обычно это выполняется в комплексных числах, но на самом деле выполнение сложных арифметических операций с достаточной точностью, чтобы гарантировать точные результаты для умножения, является медленным и подверженным ошибкам. Вместо этого мы будем использовать подход числа теоретико преобразование , который должен выполнить преобразование в целых модах N для некоторых целого числа N .
Точно так же, как есть примитивные корни из единицы любого порядка на комплексной плоскости, для любого порядка n мы можем выбрать подходящее N такое, что b является примитивным корнем из единицы порядка n в целых числах по модулю N (другими словами, b n 1 (mod N ), и никакая меньшая степень b не эквивалентна 1 mod N ).
Алгоритм будет тратить большую часть своего времени на рекурсивное умножение меньших чисел; с наивным алгоритмом это происходит в нескольких местах:
- Внутри алгоритма быстрого преобразования Фурье, где примитивный корень из единицы b многократно приводится в действие, возводится в квадрат и умножается на другие значения.
- При использовании степеней первообразного корня из единицы a для формирования весового вектора A и при умножении A или A −1 на другие векторы.
- При выполнении поэлементного умножения преобразованных векторов.
Ключевой вывод Шёнхаге-Штрассена состоит в том, чтобы выбрать N, модуль, равным 2 n + 1 для некоторого целого числа n , кратного количеству преобразуемых частей P = 2 p . Это дает ряд преимуществ в стандартных системах, которые представляют большие целые числа в двоичной форме:
- Любое значение можно быстро уменьшить по модулю 2 n + 1, используя только сдвиги и сложения , как объясняется в следующем разделе .
- Все корни из единицы, используемые преобразованием, можно записать в виде 2 k ; следовательно, мы можем умножить или разделить любое число на корень из единицы, используя сдвиг, и степень или возвести в квадрат корень из единицы, оперируя только его показателем.
- Поэлементное рекурсивное умножение преобразованных векторов может быть выполнено с использованием неациклической свертки, которая быстрее, чем ациклическая свертка, и уже имеет «бесплатно» эффект уменьшения своего результата по модулю 2 n + 1.
Для удобства рекурсивного умножения мы представим Шёнхаге – Штрассена специализированный алгоритм умножения для вычисления не просто произведения двух чисел, а произведения двух чисел по модулю 2 n + 1 для некоторого заданного n . Это не потеря общности, поскольку всегда можно выбрать достаточно большое n , чтобы произведение по модулю 2 n + 1 было просто произведением.
Оптимизация смены
В ходе работы алгоритма во многих случаях умножение или деление на степень двойки (включая все корни из единицы) можно с успехом заменить небольшим количеством сдвигов и сложений. При этом используется наблюдение, что:
- (2 n ) k ≡ (−1) k mod (2 n + 1)
Обратите внимание, что k- значное число в базе 2 n, записанное в позиционной записи, может быть выражено как ( d k −1 , ..., d 1 , d 0 ). Он представляет собой число. Также обратите внимание, что для каждого d i мы имеем 0 ≤ d i <2 n .
Это упрощает уменьшение числа, представленного в двоичном формате mod 2 n + 1 : взять крайние правые (наименее значимые) n бит, вычесть следующие n бит, добавить следующие n бит и так далее, пока биты не будут исчерпаны. Если результирующее значение все еще не находится между 0 и 2 n , нормализуйте его, добавляя или вычитая кратное модулю 2 n + 1 . Например, если n = 3 (и поэтому модуль равен 2 3 +1 = 9) и сокращаемое число равно 656, мы имеем:
- 656 = 1010010000 2 ≡ 000 2 - 010 2 + 010 2 - 1 2 = 0 - 2 + 2 - 1 = −1 8 (mod 2 3 + 1).
Более того, можно произвести очень большие сдвиги, даже не построив сдвинутый результат. Предположим, у нас есть число A от 0 до 2 n , и мы хотим умножить его на 2 k . Разделив k на n, находим k = qn + r, где r < n . Следует, что:
- A (2 k ) = A (2 qn + r ) = A [(2 n ) q (2 r )] ≡ (−1) q Shift_Left (A, r ) (mod 2 n + 1).
Поскольку A ≤ 2 n и r < n , левый сдвиг r имеет не более 2 n -1 бит, поэтому требуется только один сдвиг и вычитание (с последующей нормализацией).
Наконец, чтобы разделить на 2 k , обратите внимание, что возведение в квадрат первой эквивалентности выше дает:
- 2 2 n ≡ 1 (модуль 2 n + 1)
Следовательно,
- A / 2 k = A (2 - k ) ≡ A (2 2 n - k ) = Shift_Left (A, 2 n - k ) (mod 2 n + 1).
Алгоритм
Алгоритм следует разделению, оценке (прямое БПФ), точечному умножению, интерполяции (обратное БПФ) и объединению фаз аналогично методам Карацубы и Тоома-Кука.
Для заданных входных чисел x и y и целого числа N следующий алгоритм вычисляет произведение xy mod 2 N + 1. Если N достаточно велико, это просто произведение.
- Разделить каждое входное число в векторы Х и Y 2 K частей каждых, где- K делит N . (например, 12345678 → (12, 34, 56, 78)).
- Чтобы добиться прогресса, необходимо использовать меньшее N для рекурсивного умножения. Для этого выберите n как наименьшее целое число не менее 2 N / 2 k + k и делимое на 2 k .
- Вычислите произведение X и Y по модулю 2 n + 1, используя неациклическую свертку:
- Умножьте X и Y каждый на весовой вектор A, используя сдвиги (сдвиньте j- ю запись влево на jn / 2 k ).
- Вычислите ДПФ X и Y, используя теоретико-числовое БПФ (выполните все умножения с использованием сдвигов; для корня 2 k -й степени из единицы используйте 2 2 n / 2 k ).
- Рекурсивно примените этот алгоритм для умножения соответствующих элементов преобразованных X и Y.
- Вычислите IDFT результирующего вектора, чтобы получить вектор результата C (выполните все умножения с использованием сдвигов). Это соответствует фазе интерполяции.
- Умножьте результирующий вектор C на A −1, используя сдвиги.
- Отрегулируйте знаки: некоторые элементы результата могут быть отрицательными. Мы вычисляем максимально возможное положительное значение для j- го элемента C, ( j + 1) 2 2 N / 2 k , и если оно превышает это значение, мы вычитаем модуль 2 n + 1.
- Наконец, выполните перенос по модулю 2 N + 1, чтобы получить окончательный результат.
Оптимальное количество частей для разделения ввода пропорционально , где N - количество входных битов, и этот параметр обеспечивает время работы O ( N log N log log N ), [1] [9], поэтому параметр k должен быть установлен соответствующим образом. На практике он устанавливается эмпирически на основе размеров входных данных и архитектуры, обычно в диапазоне от 4 до 16. [8]
На шаге 2 используется наблюдение, что:
- Каждый элемент входных векторов имеет не более n / 2 k битов;
- Произведение любых двух элементов входного вектора имеет не более 2 n / 2 k битов;
- Каждый элемент свертки представляет собой сумму не более 2 k таких произведений и, следовательно, не может превышать 2 n / 2 k + k битов.
- n должно делиться на 2 k, чтобы гарантировать, что в рекурсивных вызовах выполняется условие «2 k делит N » на шаге 1.
Оптимизация
В этом разделе объясняется ряд важных практических оптимизаций, которые учитывались при реализации Шёнхаге – Штрассена в реальных системах. Он основан в основном на работе Годри, Круппы и Циммерманна 2007 года, в которой описаны усовершенствования библиотеки GNU Multi-Precision . [8]
Ниже определенной точки отсечки более эффективно выполнять рекурсивные умножения с использованием других алгоритмов, таких как умножение Тоома – Кука . Результаты должны быть уменьшены по модулю 2 n + 1, что может быть выполнено эффективно, как описано выше в разделе Оптимизация сдвига со сдвигами и добавлением / вычитанием.
Вычисление IDFT включает в себя деление каждой записи на примитивный корень из единицы 2 2 n / 2 k , операцию, которая часто комбинируется с последующим умножением вектора на A -1 , поскольку оба включают деление на степень двойки.
В системе, где большое число представлено как массив из 2 w -битных слов, полезно убедиться, что размер вектора 2 k также кратен битам на слово, выбрав k ≥ w (например, выберите k ≥ 5 на 32-битный компьютер и k ≥ 6 на 64-битном компьютере); это позволяет разбивать входы на части без битовых сдвигов и обеспечивает единообразное представление для значений mod 2 n + 1, где старшее слово может быть только нулем или единицей.
Нормализация включает добавление или вычитание модуля 2 n + 1; это значение имеет только два установленных бита, что означает, что это может быть сделано в среднем за постоянное время с помощью специальной операции.
Итерационные алгоритмы БПФ, такие как алгоритм Кули – Тьюки , хотя часто используются для БПФ векторов комплексных чисел, имеют тенденцию демонстрировать очень плохую локальность кэша с большими элементами вектора, используемыми в Шёнхаге – Штрассене. Прямая рекурсивная реализация БПФ не на месте более успешна, когда все операции помещаются в кэш за пределами определенной точки глубины вызова, но все же неоптимально использует кеш при более высоких глубинах вызовов. Годри, Круппа и Циммерман использовали технику, сочетающую четырехшаговый алгоритм Бейли с преобразованиями более высокого основания, которые объединяют несколько рекурсивных шагов. Они также смешивают фазы, максимально углубляясь в алгоритм на каждом элементе вектора, прежде чем перейти к следующему.
«Уловка квадратного корня из 2», впервые описанная Шёнхаге, заключается в том, чтобы отметить, что при условии k ≥ 2, 2 3 n / 4 −2 n / 4 является квадратным корнем из 2 mod 2 n +1, и поэтому 4 корень n-й степени из единицы (так как 2 2 n 1). Это позволяет увеличить длину преобразования с 2 k до 2 k + 1 .
Наконец, авторы тщательно выбирают правильное значение k для разных диапазонов входных чисел, отмечая, что оптимальное значение k может меняться между одними и теми же значениями несколько раз по мере увеличения размера входных данных.
Рекомендации
- ^ a b А. Шёнхаге и В. Штрассен, " Schnelle Multiplikation großer Zahlen ", Computing 7 (1971), стр. 281–292.
- ^ Martin фюрера, " Быстрее целочисленное умножение архивации 2013-04-25 в Wayback Machine ", КТВЗР 2007 Материалы, стр. 57-66.
- ^ Ван Метер, Родни; Ито, Кохей М. (2005). «Быстрое квантовое модульное возведение в степень». Физический обзор . А. 71 (5): 052320. arXiv : Quant -ph / 0408006 . DOI : 10.1103 / PhysRevA.71.052320 . S2CID 14983569 .
- ^ Обзор функций Magma V2.9, арифметический раздел Архивировано 20 августа 2006 г. на Wayback Machine : Обсуждаются практические точки пересечения между различными алгоритмами.
- ^ Луис Карлос Коронадо Гарсия, " Может ли умножение Шёнхаге ускорить шифрование или дешифрование RSA? Архивировано ", Технологический университет, Дармштадт (2005)
- ^ «MUL_FFT_THRESHOLD» . Уголок GMP разработчиков . Архивировано из оригинального 24 ноября 2010 года . Проверено 3 ноября 2011 года .
- ^ «Улучшенный класс BigInteger, использующий эффективные алгоритмы, в том числе Шёнхаге – Штрассен» . Oracle . Проверено 10 января 2014 .
- ^ a b c Пьеррик Годри, Александр Круппа и Пол Циммерманн. GMP на основе реализации крупных Integer Умножение алгоритм Шёнхаге-Штрассена Архивированные . Труды Международного симпозиума 2007 г. по символическим и алгебраическим вычислениям, стр. 167–174.
- ^ а б Р. Крэндалл и К. Померанс. Простые числа - вычислительная перспектива . Второе издание, Springer, 2005. Раздел 9.5.6: Метод Шёнхаге, с. 502. ISBN 0-387-94777-9
- ^ Дональд Э. Кнут, Искусство компьютерного программирования , Том 2: получисловые алгоритмы (3-е издание), 1997. Addison-Wesley Professional, ISBN 0-201-89684-2 . Раздел 4.3.3.C: Дискретные преобразования Фурье, стр. 305.