CC BY
271
МАТЕМАТИЧЕСКОЕ И ПРОГРАММНОЕ ОБЕСПЕЧЕНИЕ...
УДК 004.021
Е. Г. Качко
РАСПАРАЛЛЕЛИВАНИЕ АЛГОРИТМОВ УМНОЖЕНИЯ ЧИСЕЛ МНОГОКРАТНОЙ ТОЧНОСТИ
Операция умножения для длинных чисел остается объектом исследования математиков и программистов с точки зрения минимизации вычислительной сложности. Для оценки вычислительной сложности этой операции традиционно используется количество элементарных операций. При этом не учитываются свойства современных процессоров, такие как суперскалярность и многоядерность. В работе рассмотрены современные алгоритмы умножения с точки зрения возможности их распараллеливания, сделана теоретическая оценка вычислительной сложности и экспериментальная проверка с помощью Open MP и TBB. Выполнен анализ полученных результатов. Умножение; многократная точность; вычислительная сложность; параллелизм
ВВЕДЕНИЕ
Практически все современные несимметричные криптографические алгоритмы используют операцию умножения для длинных чисел. При этом длина сомножителя все время увеличивается. Так, для алгоритма Я8Л в качестве сомножителей могут использоваться числа длиной до 4096 бит, в дальнейшем и этот диапазон может быть расширен. Так как операцию умножения необходимо выполнять многократно при возведении в степень, а последняя операция используется и в алгоритмах Я8Л, и в алгоритмах на основе эллиптических кривых, то даже незначительное уменьшение вычислительной сложности для этой операции приведет к существенному уменьшению времени для формирования и проверки цифровой подписи.
Многие математики и программисты неоднократно возвращались к теме минимизации вычислительной сложности операции умножения. Из наиболее распространенных алгоритмов следует назвать алгоритмы: умножения в столбик (так называемый школьный метод), Кара-цубы [1], Шенхаге и Штрассена [2], Фюрера [3]. В табл. 1 указаны значения вычислительной сложности для каждого из вышеназванных алгоритмов.
Следует заметить, что последний алгоритм предполагает представление сомножителей как
М6
полиномов с системой счисления 2 для длин чисел до 65536 бит включительно. Для 16битных процессоров это удачное представление, для 32-битных и 64-битных такое представление приведет к существенной потере производительности. В дальнейшем будет рассмотрен во-
Контактная информация: +38050-400-54-94 Статья рекомендована к публикации программным комитетом Международной научной конференции "Параллельные вычислительные технологии 2011".
прос о возможности перехода на систему счисления 232, в данной работе этот алгоритм не рассматривается.
Т аблица 1
Вычислительная сложность основных алгоритмов умножения
Алгоритм Вычислительная сложность
Вычисление «в столбик» O (п")
Карацубы и Офмана O (nlog3)
алгоритм
Шенхаге и Штрассена алгоритм O(nlognloglogn)
Фюрера алгоритм O(nlogn2O(log*n))1
Для оценки вычислительной сложности алгоритмов все авторы используют число элементарных операций, которые необходимо выполнить. При этом не учитывается ни возможность выполнения самих операций параллельно за счет суперскалярности процессора, ни возможность выполнения отдельных функций параллельно, ни свойства кэша. Учет этих особенностей может существенно влиять на вычислительную сложность, как в сторону уменьшения, так и в сторону увеличения.
Вопросы параллельного выполнения операции умножения рассматривались и ранее. В связи со спецификой вычислительных средств этот вопрос рассматривался в условиях использования распределенных вычислений для кластерных систем [4]. В работе [5] предложен параллельный алгоритм умножения, основанный на предвычислениях, который приводит к значительному увеличению количества циклов.
1 Обозначение log*n означает loglog...log n раз по количеству рекурсий, которые необходимо выполнить
Цель данной работы - исследование основных методов умножения для чисел многократной точности.
1. МАТЕМАТИЧЕСКАЯ ПОСТАНОВКА ЗАДАЧИ
Даны 2 числа:
X = Уп-1 хВ., Y = Уп-1 у В.,
i=0 1 ’ /—I i=0У 1 ’
где B - основание системы счисления, xг■, у. -«цифры» сомножителей.
Для современных процессоров B = 232 или 264. Все экспериментальные данные получены при B = 232.
Необходимо вычислить Z = ХУ, обеспечив минимальную вычислительную сложность за счет параллельного выполнения на многоядерном процессоре.
2. ОПРЕДЕЛЕНИЕ ПОКАЗАТЕЛЕЙ ЭФФЕКТИВНОСТИ ДЛЯ АЛГОРИТМОВ УМНОЖЕНИЯ
2.1. Вычисление «в столбик»
Алгоритм 1 (Л1)
1. Обнуление произведения.
2. Для всех «цифр» второго сомножителя
• Вычисление произведения на первый сомножитель
• Сложение с текущим значением произведения.
Для простоты будем считать, что количество «цифр» для обоих сомножителей одинаково и равно п. В этом случае потребуется п операций умножения п-разядного числа на одно разрядное и сложение полученных результатов. Обозначим время выполнения операции умножения т, операции сложения (а. Время выполнения алгоритма вычисления в столбик в последовательном режиме составляет:
Тл\ = П2(^от + ?а). (1)
Алгоритм 2 (Л2). «Быстрый столбик» [6]
Для всех цифр произведения вычисляем значение цифры по формуле:
Z =У‘=х .у.; 5 =0..п-1;
_ (г'=п-1,к=п-г' _ л
Z = > х,у.; 5 = п..2п -1.
5 £ш^.=з-п+1;к=п-1 к^1
(2)
Фактически «быстрый» столбик не что иное, как вычисление значения свертки для всех цифр результата, если «цифры» - это коэффициенты многочлена.
Для этого алгоритма необходимо выполнить такое же количество операций сложения и ум-
ножения, как для алгоритма 1, но в другом порядке, поэтому Тл2 = ТЛ1.
Так как количество ядер для современных процессоров небольшое, а цель данной работы -дать практические рекомендации по использованию параллельных вычислений при умножении длинных чисел, в работе не рассматривается вариант неограниченного параллелизма, т. е. бесконечного числа параллельных ветвей.
Пусть число параллельных ветвей р, причем р < п/2 и п кратно р.
Параллельный алгоритм 1 (РЛ1)
Для обеспечения максимальной грануляции и равномерного распределения нагрузки разделим один из сомножителей на порции равной длины. В этом случае каждый поток может выполняться параллельно. Если не учитывать накладных расходов, связанных с использованием потоков, то общее время выполнения РЛ1 равно:
ТрЛ1 = п2(^ + ?а)/р. (3)
Значения ускорения и эффективности соответственно равны:
Т -
-„, = тЪ = р; ЕРЛ1 = -^ = 1. (4)
Параллельный алгоритм 2 (РЛ2)
Так как сложность вычисления «цифры» зависит от ее номера, будем использовать динамический режим распределения итераций цикла. В этом случае потребуется специальный механизм учета переносов.
2.2. Алгоритм Карацубы и Офмана (КО)
Последовательный алгоритм (КО)
Для определения числа операций запишем формулы для вычисления значения произведения по этому методу. Пусть ХЬ - младшая часть первого сомножителя, а ХН - его старшая часть, т. е. X = ХНВпП + ХЬ. Пусть ХЬ - младшая часть первого сомножителя, а ХН - его старшая часть, т. е. X = УНВп2 + ХЬ. Обозначим ZL, 2Н младшую и старшую часть произведения Z. Тогда:
Z = СВп + ЕВп/2 + Л, (5)
где
Л = ХЬ • УЬ; С = ХН • УН;
Б = (ХЬ + ХН )(УЬ + УН); (6)
Е = Б - (Л + С).
Заметим, что значение коэффициента Е не отрицательно. В работе [7] предложена формула вычисления значения произведения:
г = (вп + вп 2)хи • уи +
+Вп/2 (хи - ХЬ)(УЬ - УИ) + (7)
+(Вп/2 +1) • ХЬ • УЬ.
Формула (7) эквивалентна (5), но второе слагаемое в (7) может быть как положительным, так и отрицательным, первое и третье слагаемое вычисляются сложнее (дополнительные операции сложения), поэтому в дальнейшем рассматривается формула (5). Для исключения рекурсии для вычисления А, С, Е используем умножение «в столбик» (А1). Общее время выполнения в последовательном режиме равно:
Тко - (3«2 /4 + п + \)ґт + +(3п2/4 + 11п/2 + 6)ґа.
(8)
Параллельный алгоритм (РКО)
Схема параллельного вычисления для р > 4 представлена в табл.2.
Т аблица 2 Параллельное вычисление произведения. Алгоритм Карацубы и Хофмана
Ядро процессора
1 2 3 4
1Ь 1Н X ' = ХЬ + ХН У' = УЬ + УН
- - 1' = 1Ь + 1Н II
- - Е = В - 1 ' -
- - 1+=Е -
Операции, которые выполняются параллельно, определены в одной строке таблицы. В этом случае время выполнения параллельного алгоритма равно:
ТрКО = (п2 /4 + п + 1)?т + (п2 /4 + 5п + 5Уа. (9)
Ускорение и эффективность для алгоритма Карацубы и Офмана для р > 4:
4 _ (3п2 /4 + п + \)т + (3п2 /4 + 11п/2 + 6)а
(П /4 + п + 1)^т + (п2 /4 + 5п + 5)^а
Ерко — V4 /4.
Г'.Л.и' рко
(10)
Для двух ядер соответствующие формулы для времени выполнения и эффективности:
ТрКО = (п2 /2 + п + 1Ут + (п2 /2 + 4п + 4)^а,
(3п2 /4 + п + 1)ґт + (3п2 /4 + 11п/2 + 6)ґа
(п2 / 2 + п + 1)^т + (п2 /2 + 4п + 4)ґа
Е2 — ^рко
рко
2
2.3 Алгоритм Шенхаге и Штрассена (88). Последовательный и параллельный режимы
Этот метод основан на замене умножения больших чисел на умножение полиномов, для вычисления которых используется дискретное преобразование Фурье. Для вычисления каждой цифры полинома достаточно в качестве модуля использовать число, больше или равное 2пВ2. Для того чтобы не выполнять вычисления с числами, превосходящими 264, будем использовать простые числа, произведения которых превосходят заданный модуль. Пусть максимальное число бит сомножителя пВих = 163842. Такая длина выбрана из практических соображений. В криптографии числа большей длины не используются. Минимальный размер модуля,
т74
достаточный для таких сомножителей, равен 2 . Этот модуль гарантируется произведением трех простых чисел, меньших, чем 232, но близких к этому значению. Выбор простых чисел гарантирует наличие корней и обратного элемента. Так как значение модуля зависит только от длины сомножителей, эти числа могут быть заранее вычислены для каждой из используемых длин.
Алгоритм 88 (Последовательный режим) [7]
1. Предвычисления (вычисление простых составляющих модуля, их инверсий, корней и их степеней, и коэффициентов для обратного преобразования числа из смешанной системы счисления).
2. Для всех простых чисел
• приведение первого сомножителя по модулю;
• формирование преобразования Фурье для 1 сомножителя;
• приведение второго сомножителя по модулю;
• формирование преобразования Фурье для 2 сомножителя;
• покомпонентное умножение по модулю;
• формирование обратного преобразования Фурье.
3. Преобразование для получения результата.
Ввиду громоздкости общей формулы для определения времени выполнения, она не приведена. Рассмотрим способы параллельного вычисления для алгоритма 88. При наличии не менее трех ядер каждое ядро может выполнить полный цикл обработки для одного простого
2 Современные криптографические алгоритмы не используют больших длин
числа. Так как предварительная обработка данных не учитывается, а последовательная часть программы в этом случае невелика, то ожидаемое ускорение приблизительно равно 3. Для двухядерного процессора первое ядро выполняет преобразования для первого и второго простого числа, а второе - для третьего. В этом случае ожидаемое ускорение примерно равно 1.5.
3. РЕАЛИЗАЦИЯ АЛГОРИТМОВ УМНОЖЕНИЯ
Для реализации описанных выше алгоритмов использовался процессор Intel (R) Core (TM)2 Duo CPU E6850 @3.00GHz. Среда разработки: Microsoft Visual C++ 2008.Компиляторы: Intel(R) C++ Compiler Professional Edition 11.1, Microsoft Visual C++ 2008. Для обоих компиляторов использовался режим оптимизации по времени. Среды для разработки параллельных программ: Open MP, версия 3.0 (май 2008), TBB 3.0. Были реализованы функции для всех описанных выше алгоритмов в последовательном и параллельном режимах. Для сравнения временных характеристик в качестве эталонной библиотеки использовалась библиотека Miracl3 (версия 5.4.3 от 29.10.10). Обе библиотеки откомпилированы в режиме без использования ассемблерных вставок. Все функции реализованы для сомножителей одинаковой длины в диапазоне от 512 до 16384 бита. Длина сомножителей в таблицах задается в 32-битных словах. В таблицах результаты определяют время выполнения функций (секунды). Результаты приведены в табл. 3-6.
Как следует из табл. 3, разработанная библиотека (алгоритм А1) более эффективна, чем эталонная (Miracl). Ускорение составляет от 32 до 62% для компилятора Intel C++ и от 13 до 78% для компилятора Visual C++. Сравнение эффективности кодов, построенных с помощью двух компиляторов, показывает, что первый компилятор строит более эффективный код (ускорение составляет от 9 до 63%). Параллельное выполнение кода, реализованное с помощью Open MP и TBB, не дает выигрыша. Так как использование TBB не дает выигрыша по сравнению с Open MP для этого типа задач, в дальнейшем результаты использования этой технологии для рассматриваемого класса алгоритмов не приводятся.
Как следует из табл. 4, алгоритм А2 работает медленнее алгоритма А1 для всех длин. И, хотя параллельный режим (алгоритм PA2),
позволяет уменьшить время выполнения, этот метод продолжает быть медленнее, чем А1.
Метод Карацубы и Офмана (табл. 5) опережает обычный метод на длине 64 слова (2048 бит). Параллельный режим позволяет получить ускорение до 50%.
Алгоритм Шенхаге и Штрассена (табл. 6) проигрывает даже простому столбику на всех рассмотренных длинах. Наилучшим с точки зрения вычислительной сложности является первый алгоритм для длин сомножителей меньше 2048 и параллельный алгоритм Карацу-бы и Офмана для всех остальных длин.
Для всех алгоритмов компилятор 1Пе1 С++ создает более эффективный код. На рисунке представлены графики зависимостей времени вычислений от длины сомножителей для компилятора Ш:е1 С++ для всех рассмотренных выше алгоритмов. Ожидаемые зависимости скорости вычислений от длин не совпадают с полученными. Так, для последовательного «столбика» ожидается зависимость (1), а получается близкая к (3). Это связано с суперска-лярностью современных процессоров, когда одновременно выполняется несколько команд. То же для других алгоритмов.
Ряд 1 соответствует функции умножения («в столбик») базовой библиотеки (Миракл). Эта функция менее эффективна по сравнению со всеми функциями, кроме варианта использования функции для алгоритма Шенхаге и Штрассена для длины сомножителей не более 8192.
Ряды 2 и 3 соответствуют функции умножения («в столбик»), последовательному и параллельному вариантам соответственно. Эффективности обоих методов практически совпадают и выше эффективности функции умножения (ряд 1).
Ряды 4 и 5 соответствуют функции умножения («быстрый столбик»), для последовательного и параллельного режимов. Эти графики находятся ниже графиков для рядов 2, 3 для всех длин.
Ряды 6 и 7 соответствуют алгоритму Кара-цубы и Офмана (последовательный и параллельный режим). Этому алгоритму соответствует минимальная вычислительная сложность для длин сомножителей не менее 2048 бит.
Ряды 8 и 9 - графики для алгоритма Шенха-ге и Штрассена. Эти графики пересекают график для функции умножения (Миракл) для длины сомножителей 16384 для последовательного режима и для длины 8192 для параллельного режима.
3 http://www.shamus.ie/index.php?page=Downloads
Т аблица 3
Вычисление «в столбик»
Len Intel(R) C++ Compiler Visual C++
Miracl A1 PA1, Open MP PA1, TBB Miracl A1 PA1, Open MP PA1, TBB
16 2.79E-06 2.23E-06 2.23E-06 6.12E-05 3.91E-06 3.07E-06 3.35E-06 7.29E-05
32 6.94E-06 4.75E-06 5.31E-06 7.57E-05 8.94E-06 7.82E-06 7.26E-06 6.03E-05
64 2.43E-05 1.53E-05 1.68E-05 8.41E-05 3.10E-05 2.74E-05 2.23E-05 9.92E-05
128 9.33E-05 5.64E-05 6.22E-05 0.00012 0.000119 0.000109 8.18E-05 0.000151
256 0.00037 0.00022 0.00024 0.00033 0.000472 0.000415 0.00032 0.000318
512 0.00154 0.00085 0.00098 0.00113 0.00187 0.0017 0.00126 0.001037
Т аблица 4
«Быстрый столбик»
Len Intel(R) C++ Compiler Visual C++
A1 A2 PA2,Open MP A1 A2 PA2,Open MP
16 2.23E-06 4.75E-06 1.54E-05 3.07E-06 2.79E-06 8.10E-06
32 4.75E-06 5.59E-06 1.62E-05 7.82E-06 6.42E-06 1.03E-05
64 1.53E-05 5.03E-06 1.20E-05 2.74E-05 1.87E-05 1.93E-05
128 5.64E-05 5.59E-05 3.88E-05 0.000109 6.42E-05 4.55 E-05
256 0.00022 0.000209 0.000142 0.000415 0.00024 0.00015
512 0.00085 0.000814 0.00055 0.0017 0.00094 0.00057
Т аблица 5
Алгоритм Карацубы и Офмана
Len Intel(R) C++ Compiler Visual C++
A1 KO PKO,Open MP A1 KO PKO,Open MP
16 2.23E-06 2.79E-06 4.47E-06 3.07E-06 3.63E-06 9.77E-06
32 4.75E-06 5.59E-05 5.59E-06 7.82E-06 6.42E-06 1.42E-05
64 1.53E-05 1.54E-05 1.23E-05 2.74E-05 1.65E-05 1.98E-05
128 5.64E-05 5.17E-05 3.72 E-05 0.000109 5.59 E-05 4.61E-05
256 0.00022 0.000208 0.000137 0.000415 0.000208 0.000150
512 0.00085 0.000751 0.000510 0.0017 0.000802 0.000542
Т аблица 6
Алгоритм Шенхаге и Штрассена
Len Intel(R) C++ Compiler Visual C++
Miracl (SS) SS PSS, Open MP Miracl (SS) SS PSS, Open MP
16 3.88E-05 3.57E-05 3.24E-05 3.94E-05 3.6E-05 3.18E-05
32 8.46E-05 7.23E-05 5.11E-05 8.57E-05 7.96E-06 5.89E-05
64 0.000185 0.000160 0.000109 0.000194 0.000176 0.000125
128 0.000409 0.000357 0.000241 0.000433 0.000392 0.000267
256 0.000915 0.000786 0.000531 0.000957 0.000866 0.000587
512 0.00202 0.00174 0.00116 0.00210 0.0014 0.00127
Умножение длинных чисел
1.80Е-03
1 . 60Е-03
1 40Е-03
„ 1.20Е-03
и
1.00Е-03 | 8 00Е-04
“ 6.00Е-04
4.00Е-04 2.00Е-04 0.00Е + 00
1 2 3 4 5 6
Код длины сомножителей
Графики зависимостей времени вычислений от длины сомножителей
Ожидаемые зависимости скорости вычислений от длин, вычисленные по формулам (1-11), не совпадают с полученными экспериментально. Так, для последовательного «столбика» ожидается зависимость (1), а получается близкая к (3). Это связано с суперскалярностью современных процессоров, когда одновременно выполняется несколько команд. То же для других алгоритмов.
СПИСОК ЛИТЕРАТУРЫ
1. Карацуба А., Офман Ю. Умножение многозначных чисел на автоматах // Доклады Академии Наук СССР. 1962. Т. 145. № 2/
2. Arnold Schonhage and Volker Strassen, Schnel-le Multiplikation grosser Zahlen, Computing 7 (1971). P. 281-292.
3. Forer M. Faster integer multiplication // Proceedings of the 39th ACM Symposium on Theory of Computing. 2007. P. 57-66.
4. Fagin B. S. Large integers multiplication on massively parallel processors. (http://barryfagin.name/ Papers/FMPC90.htm).
5. Efficient Parallel Multiplication Algorithm for Large Integers. (http://www.springerlink.com/content/ hqup5uxp5np2dhmy/fulltext.pdf).
6. Василенко О. Н. Теоретико-числовые алгоритмы в криптографии. (http://window.edu.ru/ win-dow_catalog/files/r23845/book.pdf).
7. Нуссбаумер Г. Быстрое преобразование Фурье и алгоритмы вычисления сверток (http://www. toroid.ru/nussbaymerG.html).
8. Кнут Д. Искусство программирования для ЭВМ. Т. 2. Получисленные алгоритмы. М.: Мир, 1977. 314 с.
ОБ АВТОРЕ
Качко Елена Григорьевна, проф. каф. про-граммн. обеспечения ЭВМ Харьковск. нац. ун-та радиоэлектроники. Дипл. инженер-электрик (Харьковск. ин-т радиоэлектроники, 1967). Канд. техн. наук по техн. кибернетике (1977). Иссл. в обл. параллельных вычислений в криптографии.
