Программирование, 2020, № 2, стр. 67-72

1. ВВЕДЕНИЕ

Обозначения, используемые в работе

$X$ – множество переменных $X = \{ {{x}_{1}}, \ldots ,{{x}_{n}}\} $;

$\mathbb{R} = \mathbb{K}[X]$ – кольцо многочленов над кольцом $\mathbb{K}$;

$f,\;g,\;h,\;q,\;r$ – многочлены из $\mathbb{R}$;

$a,\;b,\;c$ – элементы из $\mathbb{K}$;

$F,\;G,\;H$ – конечные подмножества из $\mathbb{R}$;

$\left\langle F \right\rangle $ – идеал из $\mathbb{R}$ порожденный F;

${{\mathbb{Z}}_{{ \geqslant 0}}}$ – множество неотрицательных чисел;

$\mathbb{M} = \{ x_{1}^{{{{d}_{1}}}} \cdots x_{n}^{{{{d}_{n}}}}{\text{|}}{{d}_{i}} \in {{\mathbb{Z}}_{{ \geqslant 0}}}\} $ – множество мономов из $\mathbb{R}$;

$u,\;v,\;w,\;s,\;t$ – множество мономов или термов;

$U,\;V,\;W$ – конечные подмножества из $\mathbb{M}$;

$de{{g}_{i}}(u)$ – степень переменной ${{x}_{i}}$ в $u$;

$deg(u)$ – полная степень монома $u$;

$ \succ $ – допустимое мономиальное упорядочение с порядком переменных ${{x}_{1}} \succ \ldots \succ {{x}_{n}}$;

${\text{lt}}(f)$ – лидирующий терм относительно упорядочения $ \succ $;

${\text{lm}}(f)$ – лидирующий моном f;

${\text{lm}}(F) = \{ {\text{lm}}(f){\text{|}}f \in F\} $ – множество лидирующих мономов из F;

${\text{lcm}}(F)$ – наименьшее общее кратное множества мономов из ${\text{lm}}(F)$;

$u{\text{|}}v$ означает, что моном $u$ делит моном ${v}$.

Инволютивные базисы и базисы Грёбнера

Базис Грёбнера (и инволютивный базис) являются известным математическим инструментом, широко применяемым в задачах алгебраической геометрии, биоинформатики, вычислительной физики, робототехники и криптографии: везде, где необходимо решать системы нелинейных алгебраических, разностных и дифференциальных уравнений.

Определение и алгоритм построения базисов Грёбнера были предложены Б. Бухбергером [1].

Определение 1. [1, 2] Пусть заданно допустимое мономиальное упорядочение. Конечное подмножество $G = \{ {{g}_{1}}, \ldots ,{{g}_{m}}\} $ элементов идеала I называется его базисом Грёбнера, если

В данной работе для вычисления базиса Грёбнера применяется инволютивный подход, который описан в работах [3, 4].

Если, некоторым самосогласованным способом запретить деление по некоторым переменным, называемыми немультипликативными и разрешить по другим, называемыми мультипликативными, то мы получим сужение операции деления – инволютивное деление.

Определение 2. [3, 4] Мы будем говорить, что на множестве мономов $\mathbb{M}$ определено инволютивное деление L, если для любого конечного не пустого $U \subset \mathbb{M}$ и для любого $u \in U$ задано ${{M}_{L}}(u,U) \subseteq X$, генерирующее подмоноид $L(u,U) \subset M$, состоящий из мономов с переменными из ${{M}_{L}}(u,U)$, удовлетворяющий следующим условиям:

(a) Из $u,v \in U$ и $uL(u,U) \cap vL(v,U) \ne \emptyset $ следует $u \in vL(v,U)$ или $v \in uL(u,U)$;

(b) Из $v \in U$ и $v \in uL(u,U)$ следует $L(v,U)\, \subseteq \,L(u$, U);

(c) Из $u \in V$ и $V \subseteq U$ следует $L(u,U) \subseteq L(u,V)$ для всех $u \in V$.

Элементы ${{M}_{L}}(u,U)$ называются L-мультипликативными для $u$, элементы NM_L(u, U) = = $X{\backslash }{{M}_{L}}(u$, U) называются L-немультипликативными. Если $w \in uL(u,U)$, то $u$ называется ($L - $)инволютивным делителем $w$ и обозначается $u{{{\text{|}}}_{L}}w$. В свою очередь, моном $w$ называется ($L - $)кратным $u$.

Определение 3. [3, 4] Множество мономов называется $L$ – инволютивным, если

Определение 4. [3] Авторедуцированное множество F называется частично инволютивным до монома $v$ в допустимом мономиальном упорядочении $ \succ $ если выполнено следующее условие

Таким образом, используя конструктивное и нётерово деление (например деление Жане), пополняя частичный базис ненулевыми инволютивными нормальными формами немультипликативных продолжений элементов частично инволютивного множества F, мы в конце концов получим его инволютивный базис.

Приведем упрощенный алгоритм 1 для вычисления инволютивного базиса, не использующий критерии для пропуска редукций p к нулю. Полный алгоритм приведен в [3].

Связь между инволютивным базисом и базисом Грёбнера определяется следующей теоремой:

Теорема 1. [3] Инволютивный базис, авторедуцированный в смысле обычного деления, является редуцированным базисом Грёбнера.

2. МОТИВАЦИЯ

Получение базисов Грёбнера (и инволютивных базисов) является тяжелой вычислительной задачей как по расходу оперативной памяти, так и по времени, необходимому для завершения алгоритма. Наибольшую трудность представляют вычисления нормальной формы полинома, которые занимают более 90% времени работы. Также известно, что сложность можно снизить сначала вычислив базисы с коэффициентами из кольца модулярных чисел (или с числами Фарея [5]) и потом восстановив базис с коэффициентами из $\mathbb{Q}$. Предлагаемое представление полиномов делает вычисления более подходящими современной архитектуре процессоров и адаптированными к способам работы с оперативной памятью. Также данное представление позволяет перенести часть операций с полиномами на GPU.

3. ПРЕДСТАВЛЕНИЕ ПОЛИНОМОВ

При вычислениях с модулярными коэффициентами в известных представлениях данных для полиномов от многих переменных (списках термов, кучах [6], геобакетах [7] и т.п.) большую часть места в оперативной памяти занимают именно мономы. В задачах из области криптографии [8, 9] и задачах SAT-выполнимости могут присутствовать тысячи переменных, что ведет к “разбуханию” этих структур данных до таких масштабов, что 16Гб RAM не хватает даже на несколько секунд вычислений. В работе [10] автор уже предлагал один из методов снижения потребления памяти, развитие идей этой работы и привело к разработке нового представления полиномов, основанного на следующих идеях:

• нам не надо представлять все возможные комбинации мультииндексов, в задачах вычисления базисов Грёбнера множество возможных мономов ограничено широко известной фигурой staircase [11];

• в процессе вычислений в множестве полиномов Q, которые необходимо рассмотреть в ходе алгоритма, образуется большое количество полиномов, содержащих много одинаковых мономов.

Подсчитав количество различных мономов, встречающихся в ходе вычислений базисов на примерах из [12] (см. таблицу 1), можно утверждать, что оно очень мало по сравнению с теоретической оценкой, которую можно найти в работе [13]. Наибольшее количество мономов, встреченных автором в примере noon9: 131 860, а теоретическая оценка составляет $\left( {\begin{array}{*{20}{c}} {9 + 17} \\ 9 \end{array}} \right)$ = 3 124 550.

Представим все полиномы, которые находятся в частичном базисе и в множестве $Q$, в виде единой таблицы ${{T}_{K}}$, столбцам которой сопоставим мономы, в данный момент входящие в базис, а строкам – сами полиномы. В ячейках таблицы находятся коэффициенты при мономах.

Для обеспечения произвольного доступа к коэффициентам полинома по мономам и для упрощения работы с данными к таблице добавим еще несколько структур:

• связанный список мономов ${{L}_{ \prec }}$, отсортированный по убыванию по упорядочению $ \prec $, каждому элементу сопоставим номер его столбца в таблице ${{T}_{K}}$;

• хэш-таблицу ${{H}_{{col}}}$, по ключу-моному выдающему номер его столбца в таблице ${{T}_{K}}$;

• хэш-таблицу ${{H}_{{mon}}}$, по номеру столбца таблицы ${{T}_{K}}$ выдающую моном, ему соответствующий.

4. ДЕТАЛИ РЕАЛИЗАЦИИ

Для задач вычисления базисов Грёбнера и инволютивных базисов представление полиномов должно обеспечивать высокую эффективность следующих операций, перечисленных в порядке снижения важности:

• редукция полиномов (как по лидирующим, так и по хвостовым термам), включающая поиск подобных термов;

• умножение полинома на переменную (продолжение);

• нахождении старшего монома в полиноме

Умножение полинома на переменную фактически является переносом всех его коэффициентов в другие столбцы таблицы, соответствующие мономам, которые получаются при умножении каждого монома полинома на переменную. Это является достаточно сложной операцией (см. алгоритм 2), которая потенциально приводит к необходимости расширения таблицы при добавлении новых столбцов, соответствующих мономам, которых еще не было в полиномах частичного базиса. С течением времени все необходимые для вычислений мономы будут найдены и изменения таблицы прекратятся. Очевидно, что этот же алгоритм применяется для умножения полинома на моном.

Редукция полинома f по $g$ выполняется по алгоритму 3.

После некоторого количества редукций обязательно надо проверить плотность таблицы и осуществить компрессию убрав нулевые столбцы.

Одним из преимуществ предлагаемого представления является ориентированность на архитектуру современных процессоров и учет кэширования данных памяти. Процессоры намного быстрее обрабатывают данные из кэша, чем из оперативной памяти, разница может достигать нескольких десятков раз. Поэтому очень важно, чтобы процессор мог предсказать, какие данные и в каком порядке будут обработаны, чтобы успеть доставить их в кэш. Если выполнять редукцию по алгоритму 3, то процессор уверен, какие данные следует получить, что выражается в значительном снижении кэш-промахов. Данное утверждение можно проверить, используя инструмент valgrind-tool=chegrind [14]. Для примера reimer6 списочное представление имеет 9.6% кэш-промахов, а табличное – всего 5.4%, аналогичная картина наблюдается и для остальных примеров.

Также код, более подходящий для современных процессоров, должен содержать как можно меньше условий, которые могут приводить к сбросу конвейера выполнения инструкций и вследствие этого приводящих к уменьшению скорости работы суперскалярных блоков выполнения CPU. Поэтому редукции с нулевыми коэффициентами не обходятся, а выполняются так же, как и значащие редукции, поскольку условный их обход приводил бы к понижению производительности.

Нахождение старшего монома в полиноме. После проведения операции редукции может понадобиться обновить указатель на старший моном полинома. Нахождение его производится простой итерацией по списку ${{L}_{ \prec }}$ до первого ненулевого коэффициента в соответствующем столбце строки полинома в таблице ${{T}_{K}}$.

5. ВОЗМОЖНОСТЬ ВЫЧИСЛЕНИЯ БАЗИСОВ НА GPU

Современные GPU являются мощнейшими процессорами для вычислений общего типа, обладающими некоторыми особенностями, которые надо учитывать при адаптации алгоритмов для них. Во-первых, нужно учитывать архитектуру SIMD: на GPU запускается так называемое ядро вычислений (kernel), которое во множество потоков обрабатывает одинаковым образом входные данные, организованные в batch и warp разделы (используется терминология CUDA от фирмы NVidia). Во-вторых, наиболее эффективно происходят вычисления с большими объемами данных, поскольку GPU имеют собственную локальную память, в которую копируются исходные данные и из которой они потом должны быть скопированы обратно в память, принадлежащую CPU. Из этого возникает третье условие эффективной адаптации: передач CPU $ \leftrightarrow $ GPU должно быть как можно меньше.

На данный момент автором реализован простейший вариант адаптации: перенос части алгоритма 3 на GPU. В процессе вычислений на CPU конструируются два полинома в табличном представлении, которые передаются на GPU, редуцируются ядром

_ _global_ _ void gpuReduceEle

        (mtype* aWhat, mtype* aByWhat,

        int aN,

        const mtype aA, const mtype aB,

        const mtype aM,

        const int aK, const mtype aMod) {

    const int i = blockIdx.x * blockDim.x

            + threadIdx.x;

    if (i < aN) {

       mtype c = aWhat[i] * aA

                    + aByWhat[i] * aB;

        gpuBarretRed(c, aM, aK, aMod);

        aWhat[i] = c;

        aByWhat[i] = 0;

    }

}

и передаются обратно. Такая реализация является скорее концептом, поторый позволяет удостовериться, что вычисления инволютивных базисов и базисов Грёбнера с использованием GPU возможно. В таблице 2 показаны результаты профилирования запуска вычислений примера reimer6 на видеокарте начального уровня NVidia GT710, способной запускать одно ядро вычислений использующее до 1024 потоков обработки данных с использованием CUDA 10.0 (то, что сумма процентов больше 100 является особенностью работы профилировщика от NVidia). Видно, что собственно сами вычисления в ядре занимают лишь малую часть, основные же затраты идут на копирование данных и накладные расходы на запуск вычислений. Времена вычислений на GPU приведены в данной работе не будут, поскольку целью было показать именно принципиальную возможность вычислений используя предложенное представление данных, так как довольно очевидно, что ни одно из классических представлений полиномов не позволит получить даже такого результата.

6. РЕЗУЛЬТАТЫ ТЕСТОВЫХ ПРОГОНОВ

В таблице 3 представлены результаты сравнения реализаций алгоритма вычисления инволютивных базисов с использованием табличного и классического списочного представлений полиномов. В программах различались только представления данных, остальная реализация (ввод-вывод, работа с инволютивными делениями и т.д.) одинакова. Вычисления велись на компьютере с CPU Intel Core i5-3470@3.20GHz, RAM 8Gb под управлением Ubuntu 16.04LTS. В качестве хэш-таблицы использовалась spp::sparse_hash_map от Google [15], примеры взяты из [12]. Видно, что табличное представление сравнимо со списочным, в большинстве случаев медленнее, но в двух из представленных примеров существенно быстрее.

7. ЗАКЛЮЧЕНИЕ

Несмотря на то, что в большинстве примеров при вычислениях на CPU текущая реализация табличного представления оказалось медленнее списочного, автор считает данное представление полиномов весьма перспективным как с точки зрения ускорения расчетов на CPU, так и с точки зрения переноса вычислений на GPU.