CC BY
49
2013 Вычислительные методы в дискретной математике №4(22)
УДК 519.7
О РЕАЛИЗАЦИИ ОСНОВНЫХ ЭТАПОВ БЛОЧНОГО АЛГОРИТМА ВИДЕМАНА — КОППЕРСМИТА ДЛЯ ДВОИЧНЫХ СИСТЕМ ЛИНЕЙНЫХ УРАВНЕНИЙ НА ВЫЧИСЛИТЕЛЯХ КЛАСТЕРНОГО ТИПА
А. С. Рыжов Лаборатория ТВП, г. Москва, Россия E-mail: ryzhovalexander@mail.ru
Рассматриваются вопросы реализации наиболее трудоёмких этапов алгоритма Видемана — Копперсмита поиска решений сильно разреженных систем линейных уравнений на современных ЭВМ. Исследуются вопросы эффективной реализации отдельных операций алгоритма. Отдельно рассмотрены проблемы, возникающие при реализации алгоритма на вычислителях кластерного типа.
Ключевые слова: системы линейных уравнений, алгоритм Видемана — Копперсмита.
Введение
Необходимость нахождения решений системы однородных линейных уравнений (СОЛУ) с сильно разреженной матрицей возникает в задачах вычислительной алгебры достаточно часто. Не являются исключением и теоретико-числовые задачи, имеющие непосредственное отношение к криптографии. Например, при нахождении дискретного логарифма в поле GF(2n) (при простом p = 2n — 1) методом Копперсмита [1] необходимо найти решения системы линейных уравнений над простым полем Zp. При факторизации чисел методом решета числового поля одним из этапов является нахождение нескольких решений большой разреженной системы линейных уравнений над полем GF(2).
В настоящий момент для решения данной задачи обычно используется один из двух алгоритмов: Ланцоша — Монтгомери [2] и Видемана — Копперсмита [3]. Оба алгоритма имеют свои преимущества и недостатки, а также ограничения на используемые вычислители. Алгоритм Ланцоша — Монтгомери реализован, например, в пакете GGNFS. Он ориентирован на выполнение на одном вычислителе кластерного типа с высокоскоростным внутренним соединением узлов. Алгоритм Видемана — Копперсмита, напротив, позволяет наиболее трудоёмкие этапы выполнять на нескольких не связанных друг с другом кластерах. Именно этот алгоритм был выбран авторами работы [4] для реализации линейного этапа при факторизации 768-битного числа.
При реализации математических алгоритмов на современных вычислительных системах возникает ряд вопросов, связанных со способом представления данных, выполнением отдельных шагов и координацией работы узлов вычислителя. Знание и грамотное использование особенностей вычислителя, на котором предполагается выполнение алгоритма, помогают существенно снизить время его работы по сравнению с «лобовой» реализацией. В то же время в работах, содержащих результаты экспериментальных исследований математических алгоритмов, эти вопросы, как правило, не освещаются. Например, в [4] приводится время выполнения отдельных этапов факторизации 768-битного числа и характеристики используемых вычислительных кла-
стеров, однако не затрагиваются вопросы, связанные с реализацией этих этапов на указанных вычислителях.
В данной работе рассматриваются вопросы реализации наиболее трудоёмких этапов алгоритма Видемана — Копперсмита для двоичных систем линейных уравнений на современных вычислителях. В п. 1 излагается сам алгоритм Видемана — Копперсмита; п. 2 посвящён вопросам реализации ряда операций, выполняющихся на отдельных шагах алгоритма. В п. 3 рассмотрены нюансы выполнения алгоритма Видемана — Копперсмита на вычислителях кластерного типа; п. 4 посвящён вопросам применения алгоритма в рамках факторизации больших целых чисел методом решета числового поля.
Работу можно рассматривать как обзор подходов к программированию отдельных процедур алгоритма Видемана — Копперсмита и реализации его в целом на вычислителях кластерного типа. Автору не удалось найти публикаций, в которых бы детально анализировались вопросы эффективной реализации приводимых методов на современных многопроцессорных системах, и настоящая работа преследует цель восполнить этот пробел.
1. Алгоритм Видемана — Копперсмита
Блочный алгоритм Видемана — Копперсмита является обобщением алгоритма Видемана поиска решения разреженной СОЛУ [5], ориентированным на вычислительную технику, работающую со словами в несколько байтов. Он подробно описан в [3,6]. Приведём краткое описание основных шагов алгоритма.
Пусть А — матрица размера М х N над конечным полем Р и необходимо найти ненулевое решение СОЛУ
Ах = 0. (1)
Будем считать, что Р = ОЕ(2), параметры алгоритма т, п Є N (т ^ п) кратны длине машинного слова в битах: т = 64т/, п = 64п;, и что выполнено неравенство N ^ М + т.
Выберем случайную матрицу Z Є Рт,м размера т х М над полем Р. Пусть В Є Рм,м — матрица из первых М столбцов матрицы А СОЛУ (1), У Є Рм,п — матрица, составленная из столбцов матрицы А с номерами М+1, М+2,..., М+п. Выполняем следующие три этапа:
1) Вычислим матрицы аг = ZBгУ, і = 0,1,... , Ь (здесь Ь — параметр метода; см. далее).
і
2) Найдём такие векторы /0,..., / Є Р(п), что £ аг+ / = 0 Є Р(т) для всех і ^ 0.
3=0
1
Формальная сумма £ /3-X3 называется векторным аннулирующим многочле-
3=о
ном последовательности (аг : і = 0,...).
1-і
3) Вычислим вектор ш, составленный из векторов £ В3У/3+1 и /0 и дополненный
3=о
нулями до длины N. Можно показать [3, 6], что ш с большой вероятностью является решением системы (1).
Для нахождения векторного аннулирующего многочлена Копперсмитом предложен отдельный алгоритм [3, 6], основанный на последовательном вычислении приближений Паде. При этом в результате работы алгоритма может быть получено до п различных линейно независимых аннулирующих многочленов, что в результате при-
водит к получению до n линейно независимых решений системы (1). Из теоретиковероятностных соображений [3] достаточно иметь L = _M/mj + _M/nj + А элементов последовательности (a¿), где А œ 100.
Заметим, что на этапах 1 и 3 алгоритма при вычислении матриц W¿ = B*Y и a = ZBгУ вместо вычисления степеней матрицы B достаточно выполнять последовательное умножение матрицы Y на матрицу B: W¿+i = BW¿. При этом трудоёмкость этого умножения, равно как и вся трудоёмкость этапов 1 и 3, полностью определяется размером матрицы B и её весом, т. е. числом ненулевых элементов. Если ш — среднее число ненулевых элементов в строке матрицы B, то трудоёмкость вычисления очередной матрицы W¿+i равна О(Мш); соответственно трудоёмкость этапов 1 и 3 равна O(M 2ш).
Важно отметить, что на 64-разрядных ЭВМ умножение одного двоичного вектора на разреженную матрицу B выполняется с такой же трудоёмкостью, что и умножение на матрицу B одновременно 64 векторов. За счёт этого в блочном алгоритме Видемана — Копперсмита достигается существенное снижение трудоёмкости, так как сокращается число шагов на 1-м и 3-м этапах по сравнению с исходным алгоритмом Видемана (частный случай при m = n = 1). Например, при m = n = 64 число шагов L уменьшится в 64 раза! Если n = 64n;, то вычисления матриц W¿ можно выполнять параллельно на n1 несвязанных вычислителях, умножая по 64 столбца матрицы W¿ на матрицу B независимо на каждом вычислителе.
Замечание 1. Не стоит забывать о том, что вместо алгоритма Берлекэмпа — Мэсси, используемого на втором этапе при m = n = 1, возникает необходимость применения алгоритма Копперсмита для построения векторного аннулирующего многочлена. Трудоёмкость этого алгоритма равна O(m3L2) = O (mM2) (при m ^ n), и при сокращении числа шагов на первом и третьем этапах алгоритма растет трудоёмкость второго этапа. Однако при небольших значениях m и n трудоёмкость второго этапа в разы меньше, чем трудоёмкость этапов 1 и 3. Кроме того, в работе [7] предложена рекурсивная версия алгоритма Копперсмита, имеющая трудоёмкость порядка Оm3L log2 L. При факторизации 768-битного числа [4] для нахождения решений большой разреженной СОЛУ с помощью алгоритма Видемана — Копперсмита были выбраны параметры m =16 • 64, n = 8 • 64. Размер матрицы системы был равен N œ 193 • 106. Время работы алгоритма на трёх мощных вычислительных кластерах составило 119 сут, причём на выполнение второго этапа с использованием лишь одного из кластеров ушло менее 18 ч.
Замечание 2. Использование параметров m и n, больших 64, позволяет выполнять этапы 1 и 3 на независимых вычислителях и снижает число шагов на этих этапах, за счёт чего снижается общее время выполнения первого и третьего этапов алгоритма. Но трудоёмкость этих этапов не снижается: при m = n = 128 число шагов L в 2 раза меньше, чем при m = n = 64, но эти шаги необходимо выполнить для каждой из двух частей матрицы Y (хотя их и можно выполнять параллельно на двух независимых вычислителях).
2. Вопросы реализации первого и третьего этапов алгоритма
на современных ЭВМ
Итак, на каждом шаге первого этапа алгоритма Видемана — Копперсмита выполняются следующие операции:
— умножение матрицы W¿ размера M х n на сильно разреженную матрицу B размера
M х M: Wi+i = BW¿;
— умножение матрицы Wi+1 на матрицу Z размера m х M: ai+1 = ZWi+1.
d—1
На третьем этапе выполняется вычисление вектора v = £ Bj Y fj+1, где fj —векто-
j=0
ры длины n (коэффициенты аннулирующего многочлена). Данные вычисления можно также выполнять параллельно на независимых вычислителях, разбивая BjY на П блоков по 64 столбца, а каждый вектор fj — на П составляющих длины 64. Поскольку на втором этапе может быть получено до n различных векторных аннулирующих многочленов /¿(ж) = £ f xj, i = 1,... , n, на третьем этапе можно одновременно вычислять j=0 j di—1 (i)
n векторов vi = £ BjYfj(+i1. Соответственно на каждом шаге третьего этапа выпол-j=0 j
няются следующие операции:
— умножение матрицы Wt—1 = Bt—1Y размера M х n на сильно разреженную матрицу B размера M х M: Wt = BWt—1;
— умножение матрицы Wt на набор векторов-коэффициентов ft+1, i = 1,... , n, каждый длины n, и побитное сложение с результатами предыдущих шагов, т. е. с сум-
t—1 . /Л
мами £ Bj Yfj(i+1.
j=0 j
Очевидно, что наиболее трудоёмкой операцией каждого шага первого и третьего этапов алгоритма является умножение на разреженную матрицу B: Wi+1 = BWi. Тем не менее неэффективная реализация других операций может существенно повлиять на время работы алгоритма, поэтому необходимо оптимизировать каждую операцию, выполняемую на шагах различных этапов алгоритма.
2.1. Умножение вектора на разреженную матрицу Двоичная разреженная матрица обычно представляется в виде последовательности строк (или столбцов); каждая строка определяется числом ненулевых элементов и их номерами. Например, матрица
10 1 B = I 0 10
1 1 0
записывается в виде массива (2,1, 3,1, 2, 2,1, 2) (жирным выделены веса строк, т. е. число ненулевых элементов). При построчном хранении матрицы B (размера M х M) умножение Wi+1 = BWi, где Wi — матрица M х 64, выполняется следующим образом (здесь и далее приводится пример реализации на языке C+—+).
Алгоритм 1. Умножение (слева) на разреженную матрицу, представленную построчно в виде списка ненулевых элементов.
Вход: массив B (представление разреженной матрицы), массив W (матрица Wi, представленная в виде массива 64-разрядных чисел), число M.
Выход: массив W2 (представление матрицы Wi+1 = BWi).
1 int st, i, len; //st - номер строки матрицы B
2 //len - вес очередной строки
3 __int64 tmp; //результат умножения на очередную строку
4 //(64-разрядное число)
5 i = 0; //индекс в массиве B
6 for (st=0; st<M; st++) //проходим все строки матрицы
7 {
8 len=B[i]; //вес очередной строки
9 i ++;
10 tmp = 0;
11 while (len--) //проходим все ненулевые элементы
12 {
13 tmp = tmp ~ W[B[i]]; // ~ - побитное сложение
14 i ++;
15 }
16 W2 [st] = tmp; //очередной элемент результата
17 }
Конец алгоритма.
Хранение матрицы В можно реализовать и в виде набора весов и номеров ненулевых элементов её столбцов. Для приведённой выше матрицы В её представление по столбцам — (2,1, 3, 2, 2, 3,1,1). В этом случае умножение Жг+1 = ВЖг можно выполнить с помощью следующего алгоритма.
Алгоритм 2. Умножение (слева) на разреженную матрицу, представленную по столбцам в виде списка ненулевых элементов.
Вход: массив В (представление разреженной матрицы по столбцам), массив Ч (матрица Жг, представленная в виде массива 64-разрядных чисел), число М.
Выход: массив Ч2 (представление матрицы Жг+1 = ВЖг).
1 i nt st ,i,len ; //st - номер столбца матрицы B
2 //len - вес очередного столбца
3 __int64 tmp; //очередной элемент массива W
4 for (st=0; st<M; st++) W2[st] = 0; //обнуляем результат
5 i= 0; //индекс в массиве B
6 for (st=0; st<M; st++) //проходим все столбцы матрицы
7 {
8 len=B [i]; //вес очередного столбца
9 tmp = W[st]; //следующий элемент массива Ч
10 i ++;
11 while (len--) // проходим все ненулевые элементы
12
13 W2[B[i]] = W2[B[i]] ~ tmp; //побитное сложение
14 i ++;
15 }
16 W2 [str] = tmp ; //очередной элемент результата
17 }
Конец алгоритма.
Если ш — среднее число ненулевых элементов в строке (столбце) матрицы В, то трудоёмкость обоих алгоритмов равна О(Мш). Проведены исследования по сравнению времени работы алгоритмов 1 и 2 для различных размеров входных данных. Результаты экспериментов приведены в табл. 1 (время работы указано в миллисекундах).
Из табл. 1 видно, что алгоритм 1 выигрывает по скорости работы. Это объясняется тем, что внутри основного цикла алгоритма 2 выполняется долгая (по сравнению с остальными) операция записи в память (прибавление элемента массива Ч к соответствующему элементу результата). В алгоритме 1 промежуточное значение очередного элемента выходного массива Ч2 хранится в переменной Шр, и если в машинном коде
Таблица 1 Время работы алгоритмов 1 и 2
M ш Алгоритм 1 Алгоритм 2
100 000 100 24 29
200 000 100 51 62
500 000 100 137 164
1000 000 100 484 648
2 000 000 100 1445 2269
5 000 000 100 4537 6949
10 000 000 100 11677 17135
200 000 500 251 303
500 000 200 272 321
1000 000 100 484 648
2 000 000 50 712 1050
эта переменная реализована с помощью одного из регистров процессора, то запись в неё происходит во много раз быстрее, чем запись в ячейку оперативной памяти.
Следует также отметить разницу во времени работы алгоритмов в четырёх последних экспериментах. Несмотря на то, что теоретическая трудоёмкость алгоритмов O(Mw) для выбранных параметров одинакова, умножение на более плотные матрицы работает быстрее за счёт того, что реже выполняется операция записи в память. Этим также объясняется нелинейный рост времени работы алгоритмов относительно размера матрицы при одинаковом среднем числе ненулевых элементов в строке (столбце).
Алгоритм 1 допускает эффективное распараллеливание по строкам матрицы B: умножение на строки с номерами n + 1, n + 2,..., n + k можно выполнять независимо от умножения на другие строки; при этом будут получены k элементов выходного массива W2 с номерами n + 1,n + 2,...,ni + k. Тем не менее из-за особенностей работы потоков с общей оперативной памятью не рекомендуется распараллеливать данный алгоритм на большое число потоков.
2.2. Умножение двух векторов В исходном алгоритме Видемана (т. е. при m = n = 1) на каждом шаге первого этапа после умножения вектора на матрицу вычисляется скалярное произведение двух векторов длины M: a = ZWj. Трудоёмкость вычисления скалярного произведения равна 2M операций в поле GF(2). В общем случае в алгоритме Видемана — Копперсмита требуется перемножить две матрицы размеров m х M и M х n. При фиксированных m и n это также можно сделать за O(M) операций, используя тривиальный алгоритм. Однако с увеличением параметров m и n трудоёмкость тривиального алгоритма быстро растет, поэтому для ускорения данной процедуры целесообразно использовать более быстрые алгоритмы, например LUT-алгоритм на основе таблиц поиска (Look-Up Tables). Описание LUT-алгоритма можно найти, например, в [8, гл. 6].
Идея метода заключается в следующем. Пусть, для начала, m = n = 8, а матрицы Z и W размеров 8 х M и M х 8 представлены в виде массивов байтов Z и W длины M. Результат умножения a = ZW есть матрица 8 х 8, представляемая массивом из 8 байтов. Каждая строка результата есть линейная комбинация строк второго множителя: a^ = £ ZjWj. Выделение отдельных элементов Zj, т. е. отдельных битов j
элементов массива Z, является достаточно долгой операцией. Вместо этого удобнее сохранять промежуточный результат вычисления линейных комбинаций £ Zj Wj сра-
j
зу для всех i = 1,... , 8, а именно: в дополнительный массив LUT из 28 = 256 байтов в ячейку с номером k G {0,..., 255} запишем сумму строк второго множителя с такими номерами j G {1,... , M}, для которых значение столбца первого множителя (т. е. элемента массива Z) равно k:
LUTfc = £ Wj.
je{1,...,M }:
Zj =k
Создание массива LUT выполняется за один проход по массивам Z и W: для j = = 1,... , M выполняем побитное сложение LUTzj := LUTzj Ф Wj. Получение результата умножения из массива LUT тривиально. Например, первая строка результата есть линейная комбинация элементов массива LUT, первый (младший) бит номеров которых равен 1.
Для случая m = n = 64 невозможно выделить дополнительную память LUT объёмом 264 64-разрядных элементов. Вместо этого будем использовать 8 различных таблиц LUT, разбив вычисления по номерам байтов в строке результата. При больших значениях m и n удобнее решать задачу, разбивая на подзадачи с m = n = 64.
Итоговый алгоритм на языке C+—+ можно записать следующим образом.
Алгоритм 3. Умножение матриц размеров 64 х M и M х 64: а = ZW.
Вход: массивы Z и W 64-разрядных чисел, представляющие матрицы Z и W соответственно, число M.
Выход: массив A из 64-х 64-разрядных чисел, представляющий результат (матрицу а размера 64 х 64).
1 __int64 LUT [256 * 8] ;
2 __int64 zj , wj;
3 int i, j, k;
4 for (j = 0; j < 256*8; j++) LUT [j] = 0;
5 for (j =0; j <M; j ++)
6 {
7 zj = Z[j];
8 wj = W[j];
9 for (i=0; i<8; i++)
10 LUT [i *256 + (zj >> 8*i) & 0xFF ] ~= wj;
11 //при m=n=8 здесь вместо цикла по i было бы только
12 //одно побитное сложение: LUT [Z[j]] ~= W[j];
13 }
14 for (i = 0; i < 64; i++) a[i] = 0;
15 for (i = 0; i < 8; i++)
16 {
17 for (j = 0; j < 256; j++)
18 if ((j >> i) & 1)
19 for (k = 0; k < 8; k++)
20 a [8*k + i] ~= LUT [k*256 + j];
21 }
Конец алгоритма.
Алгоритм 3 работает на ЭВМ существенно быстрее тривиального алгоритма (см. табл. 2). Выигрыш достигается за счёт того, что результатом является небольшая двоичная матрица размера 64 х 64, а также за счёт эффективного использования дополнительной памяти и многоразрядности ЭВМ.
Таблица 2 Время работы алгоритма умножения матриц 64 х М и М х 64
М Тривиальный алгоритм Алгоритм 3
100 000 4,69 0,42
1000 000 47,02 4,27
10 000 000 470,03 42,59
100 000 000 4715,9 436,8
2.3. Умножение на коэффициенты аннулирующего
многочлена
^-1
При вычислении вектора V = £ В3У/3+1 на каждом шаге третьего этапа необхо-
3=0
димо, помимо вычисления матрицы Ж, = В3 У, реализовать умножение этой матрицы размера М х п на вектор /3+1 длины п. Если п = 64 и Ж, представляется в виде массива длины М, состоящего из 64-разрядных чисел, а вектор /3+1 записан как одно 64-разрядное число, то умножение Ж, на /3+1 сводится к операциям побитного умножения и вычисления чётности двоичного веса Хэмминга 64-разрядных чисел. С помощью сдвига и побитного сложения в 2 раза уменьшается размер задачи вычисления чётности двоичного веса вектора, а для 16-разрядных чисел соответствующие значения можно вычислить заранее и хранить в памяти. В результате умножение Ж, на /3+1 и сложение с результатом предыдущих шагов выполняется за 6М операций с 64-разрядными числами и М операций обращения к массиву двоичных весов 16-разрядных чисел. Однако при наличии достаточного объёма оперативной памяти можно выполнить требуемое умножение и сложение более эффективно. Поскольку отображение х : Уб4 ^ VI, реализующее вычисление чётности двоичного веса Хэмминга 64-разрядного вектора, является линейным, можно на каждом шаге вычислять только побитное произведение /3+1 и элементов матрицы Ж, и складывать их (побитно) с соответствующими результатами предыдущих шагов, а отображение х вычислить после выполнения всех d шагов третьего этапа. При этом необходимый объём памяти для хранения промежуточных сумм £ В3У /3+1 возрастёт в 64 раза.
,
Заметим, что требуемый объём оперативной памяти можно уменьшить в 2 раза с помощью М сдвигов и М сложений на каждом шаге.
3. Реализация первого и третьего этапов алгоритма на вычислителях кластерного типа
Рассмотрим особенности, возникающие при реализации алгоритма Видемана — Копперсмита на кластерах. Под кластером, или вычислителем кластерного типа, будем понимать некоторый набор ЭВМ (узлов), соединённых между собой высокоскоростной сетью. Поскольку в этом случае отсутствует возможность использования быстрой общей памяти, необходимо решать вопросы о распределении обрабатываемых данных по узлам кластера, координации работы узлов, а также о способе обмена данными между узлами.
3.1. Распределение данных и распараллеливание вычислений
При выборе способа представления данных и распределения их между узлами кластера необходимо исходить из соображений эффективной реализации наиболее трудо-
ёмких операций алгоритма. В случае алгоритма Видемана — Копперсмита наиболее трудоёмкой операцией, выполняемой на каждом шаге первого и третьего этапов, является умножение векторов на разреженную матрицу. Если использовать алгоритм 1, то данная операция легко распараллеливается между узлами кластера путём разбиения разреженной матрицы на блоки строк примерно одинакового веса. Каждый узел хранит свою порцию строк разреженной матрицы В, а также всю матрицу Ж,. После умножения на каждом узле находится соответствующая часть результата, т. е. матрицы Ж,+1 = ВЖ,. Перед следующим умножением необходимо выполнить обмен между узлами, чтобы на каждом из них вновь была полная копия правого множителя.
Строки матрицы В необходимо распределять по узлам кластера таким образом, чтобы каждый узел тратил примерно одинаковое время на умножение матрицы Ж, на свою порцию строк матрицы В. Как правило, используются узлы кластера с одинаковой производительностью, поэтому время умножения определяется весом строк разреженной матрицы В, доставшихся данному узлу, и их количеством. Важно учитывать время записи результата умножения Ж, на одну строку матрицы В, поэтому умножение на несколько «тяжелых» строк матрицы В выполняется быстрее, чем умножение на большее количество более «легких» строк с таким же суммарным количеством ненулевых элементов (см. п. 2.1).
Умножение двух матриц Z и Ж, для получения очередного элемента последовательности Крылова а, = ZЖ:,• = ZB3 У, особенно с помощью алгоритма 3, выполняется намного быстрее, чем умножение Ж, на разреженную матрицу В, поэтому данную операцию можно не распараллеливать и выполнять на одном из узлов после формирования на нем полной копии матрицы Ж,. При этом на выбранном узле необходимо хранить всю матрицу Z.
Вычисление сумм £ В, У /,+1 на шагах третьего этапа также является достаточно ,
быстрой операцией, к тому же её легко распараллелить, так как одновременно может вычисляться до п таких сумм.
3.2. П е р е д а ч а д а н н ы х и к о о р д и н а ц и я р а б о т ы у з л о в
Для выполнения алгоритма на кластере необходимо выделить один из узлов для координации действий. На этом управляющем узле хранится результат вычислений первого этапа — матричная последовательность (а^). Этот же узел выполняет умножение матриц Z и Ж, на каждом шаге первого этапа.
Управляющий узел должен следить за ходом выполнения шагов алгоритма. Он даёт команду остальным узлам на выполнение умножения на разреженную матрицу, ждёт отчёта узлов о завершении умножения, затем инициирует процедуру обмена данными между узлами. Как отмечалось выше, после умножения матрицы Ж, на разреженную матрицу В на каждом узле будет получена часть матрицы Ж,+і = ВЖ,, соответствующая доли матрицы В на этом узле. Перед следующим умножением, а также перед вычислением матрицы а,+1 необходимо получить полную копию Ж?+1 на каждом узле кластера.
Наиболее эффективным способом обмена данными между узлами в этом случае представляется передача данных по кольцу. Кратко изложим идею этого способа. Пусть имеется к узлов и на г-м узле хранится соответствующая часть X некоторого массива X, г = 1,... , к. Заметим, что размер части массива X на некоторых узлах может быть равен нулю. Передача осуществляется следующим образом. Сначала каждый узел передаёт предыдущему узлу свою часть массива X, т. е. узел с номером г передает узлу с номером г — 1 часть Х^ (первый узел передаёт Х1 на узел с номером к). По-
сле этого каждый узел передает предыдущему только что полученную от следующего узла часть: i-й узел передает часть Xi+1. И так далее. На j-м шаге узел с номером i передает часть Xi+j-1 и получает часть Xi+j, j = 1,..., k. Через k — 1 шагов каждый узел будет иметь весь массив X. Применительно к алгоритму Видемана — Копперсмита, при равномерном распределении данных между k узлами каждый узел отправляет k1
и получает —-— 64M битов. Заметим, что нет необходимости координировать данную k
процедуру на различных шагах: каждый узел самостоятельно продолжает передачу и прием данных в необходимом объеме, если перед началом процедуры ему известен размер части на предыдущем узле и общий размер массива X.
При использовании технологии MPI выполнить подобный обмен между процессами одной группы можно с помощью готовой процедуры MPI_ALLGATHER (см., например, [9]). Реализация обмена была подготовлена путем использования системных функций для установления соединения и пересылки данных между двумя узлами (библиотека Winsock под ОС Windows и sys/socket под ОС Linux). Заметим, что если сетевые карты узлов поддерживают режим полного дуплекса (что типично для современных вычислителей), то удается достичь максимально возможной скорости передачи для данной сети. Например, при M = 10 000 000 и равномерном распределении матрицы между k =10 узлами, соединенными сетью gigabit Ethernet (передача до 109 бит/с), сетевой обмен между узлами на каждом шаге выполняется примерно за 585мс (каждый узел отправляет и получает 576 • 106 битов).
В случае, когда умножение узлами матрицы Wj на их порции разреженной матрицы B выполняется за время, сравнимое с временем обмена результатами умножения, имеет смысл выполнять умножение и передачу одновременно. При этом каждый узел начинает выполнять умножение и сразу же передает либо вычисленную им часть матрицы Wj+1, либо полученную от соседнего узла. Из-за нагрузки на оперативную память снижение времени работы алгоритма в 2 раза невозможно, однако при оптимальном балансе между временем умножения и передачи удается получить выигрыш примерно в 1,5 раза. Заметим, что подобный подход исключает возможность обмена результатами вычисления Wj+i между узлами с помощью готовых решений технологии MPI.
Результаты работы третьего этапа алгоритма Видемана — Копперсмита могут занимать большой объем памяти, и в этом случае их можно хранить на отдельно выделенных для этого узлах кластера. По окончании выполнения третьего этапа эти узлы передают результаты вычислений на управляющий узел.
3.3. Реализация проверки и сохранения п р о м е ж у т о ч н ы х р е з у л ь т а т о в
При решении задач больших размеров важными аспектами реализации алгоритма являются проверка и сохранение промежуточных результатов. Наибольшей вычислительной и емкостной трудоемкостью обладает операция умножения на разреженную матрицу. Наибольший объем данных, передаваемых между узлами кластера, также связан с этой операцией. Поэтому проверка правильности вычисления матриц Wj, i = 1,2,... , представляет наиболее важную задачу с точки зрения контроля правильности работы алгоритма.
Проверка правильности вычисления матриц Wj = BjY может быть реализована следующим образом. Пусть itest G N — параметр метода, Ztest — случайная матрица размера mtest х M. Вычислим значение ZtestBitest. Далее выполняем шаги первого или
?
третьего этапов алгоритма и на шаге i = itest выполняем проверку: (ZtestBitest ) Wo =
= . Если равенство не выполняется, то в процессе выполнения данного этапа
алгоритма была допущена ошибка и необходимо начать выполнение этого этапа сначала. В случае выполнения сохраняем значение и выполняем следующие ша-
?
гов, после чего (на шаге % = 2%^^ снова выполняем проверку: (^е81В=
?
= . Если равенство не выполняется, возвращаемся к шагу % = иначе
сохраняем и продолжаем выполнение шагов алгоритма. И так далее. На ша-
?
ге % = к ■ гtest, к = 1, 2,..., проверяем равенство ^) Ж^-!)^ = и
сохраняем значение либо возвращаемся на шаг % = (к — 1^е^.
Для обнаружения ошибок достаточно использовать матрицу ^е^ размеров 64 х М. В качестве ^е^ можно взять первые 64 столбца случайной матрицы Z. Однако иногда вместо случайной матрицы Z выбирают матрицу, состоящую из первых т столбцов единичной матрицы (см. далее), и в этом случае Ztest лучше выбирать случайно. Параметр %^(| выбирают исходя из размеров задачи и надёжности вычислителя. При факторизации 768-битного числа, описанной в [4], проверка выполнялась каждые %^(| = 2!4 шагов (всего на первом этапе было вычислено Ь ^ 565000 элементов последовательности (аг), причём в работе не сказано о наличии сбоев при вычислении).
При выполнении первого этапа алгоритма Видемана — Копперсмита в качестве промежуточных результатов выступают элементы последовательности (аг) и матрицы = ВгУ, поэтому их сохранение реализуется автоматически при выполнении проверки правильности вычислений. В случае сбоя на шаге % (например, при отключении питания вычислителя и т.п.) достаточно вернуться к шагу к ■ %^^, где к = |_%^е^_|, и продолжить выполнение алгоритма. На третьем же этапе необходимо выполнять сохранение промежуточных сумм £ В3У/3+!. Эта процедура также выполняется каж-
3
дые шагов, и промежуточные суммы можно хранить либо на узлах, которые эти суммы вычисляют, либо на управляющем узле кластера. В последнем случае необходимо выполнить передачу этих значений с соответствующих узлов кластера на управляющий узел. Заметим, что если применяется метод ускорения вычислений сумм £ В3У/3+!, описанный в п. 2.3, то объём передаваемых и сохраняемых данных не воз-
з
растает, так как перед передачей и сохранением промежуточных результатов можно выполнить вычисление чётности двоичных весов Хэмминга и хранить собственно значения £ В3 У/з+1.
3
3.4. Н е к о т о р ы е с п о с о б ы с н и ж е н и я в ы ч и с л и т е л ь н о й и ёмкостной сложности алгоритма Приведём некоторые незначительные изменения алгоритма Видемана — Коппер-смита, которые позволяют снизить требуемый объём оперативной памяти вычислителя и сократить трудоёмкость отдельных операций.
Первое изменение связано с выбором матрицы Z. В исходном описании алгоритма сказано, что матрица Z должна выбираться случайным образом. При этом появляются затраты на её хранение и умножение на матрицы Жг. Действительно, если матрицу У можно разбить на подматрицы по 64 столбца и выполнять первый и третий этапы алгоритма для каждой из этих подматриц на несвязанных между собой вычислителях, то матрица Z должна храниться на каждом вычислителе целиком, что при большом значении т приводит к существенным затратам оперативной памяти. Чтобы этого избежать, в качестве Z можно использовать первые т строк единичной матрицы. В этом случае матрицу Z хранить не нужно, а результат умножения аг = ZWi — это первые
т строк матрицы Жг. Следует отметить, что в этом случае нарушается обоснование надёжности алгоритма [3], основанное на вероятностных соображениях, однако многочисленные эксперименты показывают, что для невырожденных матриц А алгоритм работает нормально. В качестве компромисса можно использовать матрицу Z, столбцы которой имеют вес 1. Если т ^ 256, то для хранения такой матрицы достаточно М байтов вместо тМ/8, а умножение на неё матрицы выполняется, очевидно, более эффективно, чем в общем случае.
Более важным является снижение трудоёмкости умножения матрицы на разреженную матрицу В. При изложении алгоритма в п. 1 было сказано, что В состоит из первых М столбцов матрицы А. Тем не менее при неравномерном распределении веса столбцов имеет смысл составлять матрицу В из столбцов меньшего веса; иными словами, перед выполнением алгоритма выполнить перестановку столбцов матрицы А (перенумерацию неизвестных), чтобы снизить вес матрицы В и тем самым снизить трудоёмкость умножения на эту матрицу.
Вес строк матрицы А также может быть распределён неравномерно. В результате работы алгоритма Видемана — Копперсмита получается не одно, а сразу п решений системы. Если перед выполнением алгоритма отбросить одно уравнение системы, т. е. одну строку матрицы А, и найти п решений для такой системы, то можно получить не менее п — 1 решений исходной системы. Действительно, пусть и^,... , ^^-!, ид — решения урезанной системы, не удовлетворяющие отброшенному уравнению, 1 ^ к ^ п. Тогда к—1 векторов и! ®ик,... , ик-! ®ик будут удовлетворять и отброшенному уравнению, и, очевидно, всем остальным. Следовательно, если требуется получить небольшое число решений (например, одно), наиболее тяжёлые строки матрицы А (т. е. строки большого веса) можно изначально исключить; при исключении £ строк из полученных п решений можно сформировать не менее п — £ решений исходной системы уравнений.
4. Особенности применения алгоритма Видемана — Копперсмита в рамках метода решета числового поля факторизации целых чисел
Как уже было отмечено, нахождение решений большой разреженной системы однородных линейных уравнений над полем СЕ(2) является одним из этапов метода решета числового поля, применяемого для факторизации больших целых чисел. При этом матрица решаемой системы не является случайной и имеет ряд ярко выраженных особенностей, что необходимо учитывать при выполнении алгоритма Видемана — Копперсмита.
Не вдаваясь в подробности, можно сказать, что каждая строка матрицы соответствует очередному простому числу, а элементы строки представляют чётность показателя степени, с которой данное простое число входит в разложение некоторых выбранных целых чисел. Для случайного натурального числа п и простого р вероятность того, что рк делит п, а рк+! не делит п, примерно равна (р — 1)/рк+!, к = 0,1,... (понятно, что при к > 1с^рп эта вероятность равна нулю). Соответственно вероятность того, что элемент строки, соответствующей простому числу р, равен 1, примерно равна £ Р =-----------. Поэтому первые строки матрицы имеют достаточно большой вес
д=! Р Р + 1
(для первой строки, соответствующей р = 2, вес равен примерно N/3), а затем с ростом номера строки вес строки быстро снижается. Используя идею п. 3.4, можно отбросить несколько первых строк (т. е. уравнений), выполнить алгоритм, а затем из полученных решений урезанной системы составить несколько решений исходной системы.
Ещё один важный момент касается предварительной обработки матрицы решаемой разреженной СОЛУ. Например, если вес строки равен 1, то можно сразу выполнить исключение соответствующей переменной из решаемой системы. Аналогично, если некоторая переменная встречается только в одном уравнении (т. е. вес столбца равен 1), то соответствующие строку и столбец можно удалить перед выполнением алгоритма поиска решений. Применяются и другие операции (см., например, [10]), приводящие к снижению размеров M и N системы, но увеличивающие средний вес строки ш матрицы системы. Такая предварительная обработка матрицы называется фильтрацией и выполняется до того момента, пока эти операции приводят к снижению числа M2ш, определяющего трудоёмкость первого и третьего этапов алгоритма Видемана — Копперсмита (при реализации алгоритма на кластере необходимо также учитывать общее время передачи данных по сети, которое пропорционально M2). При факторизации 768-битного числа [4] исходная матрица 35 288 334 017 х 47 762 243 404 была за 15 дней ужата до размеров 192 795 550 х 192 796 550 со средним весом строки ш = 144.
При выполнении фильтрации решение о её прекращении принимается, исходя из оптимальных соотношений между размерами матрицы и её весом, приводящих к минимальной трудоёмкости алгоритма нахождения решений СОЛУ. Для одной и той же матрицы данные соотношения могут быть различными при использовании разных вычислителей для выполнения алгоритма Видемана — Копперсмита. Поэтому при решении сильно разреженной СОЛУ фильтрация является неотъемлемым этапом и должна согласовываться с выбранным вычислителем и параметрами алгоритма Видемана — Копперсмита.
Заключение
Рассмотрены вопросы, связанные с реализацией первого и третьего этапов алгоритма Видемана — Копперсмита. Наиболее трудоёмкой операцией этих этапов является умножение блока векторов на разреженную матрицу. Эта операция может быть эффективно распараллелена путем независимого умножения на различные строки разреженной матрицы, однако при реализации на вычислителях кластерного типа требуется дополнительный обмен между узлами для получения результата умножения.
Уделено внимание и другим операциям первого и третьего этапов алгоритма Ви-демана — Копперсмита, поскольку их неэффективная реализация может существенно увеличить время работы алгоритма. В рамках исследования особенностей реализации алгоритма на кластерах рассмотрены вопросы проверки и сохранения промежуточных результатов вычислений для защиты от сбоев в работе вычислителя. Приведены также некоторые особенности, возникающие при использовании алгоритма Видемана — Копперсмита в рамках факторизации больших целых чисел методом решета числового поля.
ЛИТЕРАТУРА
1. Coppersmith D. Fast evaluation of logarithms in fields of characteristic two // IEEE Trans.
Inform. Theory. 1984. V.IT-30(4). P. 587-594.
2. Montgomery P. L. A block Lanczos algorithm for finding dependencies over GF(2) //
EUROCRYPT’95. LNCS. 1995. V.921. P. 106-120.
3. Coppersmith D. Solving linear equations over GF(2) via block Wiedemann algorithm // Math.
Comp. 1994. V. 62(205). P. 333-350.
4. Kleinjung T., Aoki K., Franke J., et al. Factorization of a 768-bit RSA modulus // CRYPTO-
2010. LNCS. 2010. V.6223. P. 333-350.
5. Wiedemann D. H. Solving sparse linear equations over finite fields // IEEE Trans. Inform. Theory. 1986. V. IT-32(1). P. 54-62.
6. Рыжов А. С. О реализации алгоритма Копперсмита для двоичных матричных последовательностей на вычислителях кластерного типа // Прикладная дискретная математика. 2013. №3. С. 112-122.
7. Thome E. Subquadratic computation of vector generating polynomials and improvement of the block Wiedemann algorithm // J. Symb. Comput. 2002. No. 33. P. 757-775.
8. Ахо А., Хопкрофт Д., Ульман Дж. Построение и анализ вычислительных алгоритмов. М.: Мир, 1979. 536 с.
9. www.open-mpi.org/doc/v1.6 — Open MPI v1.6.4 documentation. 2013.
10. Cavallar S. Strategies for filtering in the number field sieve // Proc. ANTS IV. LNCS. 2000. V. 1838. P. 209-231.
