Программирование, 2019, № 5, стр. 56-66

I. ВВЕДЕНИЕ

Матрицы с элементами, принадлежащими самым разным кольцам и полям, используются во всех областях математики и в приложениях. В этой статье речь идет о матрицах, элементы которых принадлежат кольцу линейных разностных скалярных операторов над некоторым разностным полем $\mathbb{K}$ с автоморфизмом (сдвигом) $\sigma $. Предполагается, что поле $\mathbb{K}$ имеет характеристику 0. Изложение разворачивается вокруг задач проверки обратимости данной матрицы обсуждаемого типа и нахождения обратной матрицы, сколь скоро она существует. Мы даем беглый обзор известных алгоритмов, основанных на регуляризации, т.е. получении в невырожденном виде сопоставляемых исходной операторной матрице ведущих и трейлинговых матриц, элементы которых принадлежат $\mathbb{K}$. В некоторых алгоритмах регуляризируются не ведущая и трейлинговая, а фронтальная и тыльная матрицы (их элементы также принадлежат $\mathbb{K}$; в [1] использован общий термин “выявляющие матрицы”).

Таким образом, для матриц обсуждаемого вида уже предложено некоторое число алгоритмов обращения, и мы приводим перечень этих алгоритмов в пп. III, IV. Для них проведено исследование сложности, где сложность рассматривается в двух вариантах:

1) арифметическая сложность, т.е. сложность в худшем случае по числу арифметических операций в поле $\mathbb{K}$,

2) сдвиговая сложность, – тоже в худшем случае, но по числу применений $\sigma $ и ${{\sigma }^{{ - 1}}}$.

Некоторые алгоритмы этого множества уже были реализованы в виде Maple-процедур. Мы предоставляем недостающие реализации тех алгоритмов, которые имеют сравнительно небольшие сложности и проводим экспериментальное сравнение этих алгоритмов по фактическим затратам времени. Интересно, что алгоритмы, имеющие меньшую сложность, не всегда оказываются самыми быстрыми, и это, разумеется, не является ни в какой мере парадоксом, так как сложность, будь то арифметическая или сдвиговая, или даже сложность по общему числу тех и других операций в худшем случае, во-первых, не характеризует всех затрат на выполнение алгоритма, и, во-вторых, сравнение затрат в худшем случае не охватывает всех мыслимых случаев, а с “худшим случаем” мы можем сталкиваться весьма редко, и охватываемые набором тестов случаи могут не попадать в разряд худших. Нелишне еще заметить, что когда мы сравниваем сложности, исходя из их асимптотических оценок, то не всегда получаем адекватное представление о значениях той или иной сложности для данных не слишком большого (не “астрономического”) размера.

Для выбора наиболее удобного, например, – наиболее быстрого, алгоритма из некоторого множества, сложностные характеристики алгоритмов этого множества играют, естественно, немаловажную роль. Но если сложности разных алгоритмов являются близкими величинами, для которых мы располагаем лишь асимптотическими оценками, то эти оценки дадут нам возможность лишь предварительного отбора. Последующее экспериментальное сравнение (тестирование) оказывается очень желательным, особенно когда тестовые задания имеют размеры “типичные” для тех ситуаций, на которые ориентирована реализация алгоритмов.

Итак, в статье мы даем как некоторые сложностные оценки для алгоритмов рассматриваемого множества, так и результаты экспериментального сравнения этих алгоритмов.

Здесь надо сказать, что в [2] мы указывали ссылку на страницу в интернете, содержащую реализацию алгоритма с наименьшими в асимптотическом смысле арифметической и сдвиговой сложностью. Отмечалось, что в реализации было сделано некоторое отступление от самого алгоритма в части регуляризации ведущей и трейлинговой матриц с помощью некоторого специального варианта алгоритма EG ([3, 4]). Однако выбор подходящего варианта этой регуляризации оставлен в [2], фактически, “за кадром”. В настоящей статье сравнению алгоритмов на уровне реализации уделено особое внимание.

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

Завершим это введение в статью несколькими примечаниями и соглашениями.

Обычно в случаях операторных матриц вместо “обратимая матрица” употребляется термин унимодулярная матрица. Мы также будем пользоваться этим термином.

Алгоритмы проверки унимодулярности и построения обратной матрицы для дифференциального случая, когда $\mathbb{K}$ является дифференциальным полем характеристики 0 с дифференцированием $\delta = '$ и когда элементы матриц суть скалярные линейные дифференциальные операторы над $\mathbb{K}$, рассматривались в [5]. Для данной операторной матрицы $L$ как дифференциальные, так и обсуждаемые ниже разностные алгоритмы вычисляют размерность пространства решений ${{V}_{L}}$ соответствующей системы уравнений $L(y) = 0$. Делается это в предположении, что компоненты решений принадлежат связанному с L адекватному (см. п. II) расширению поля $\mathbb{K}$. Как было установлено ранее (см. [6]), операторная матрица L полного ранга (строки L независимы над кольцом скалярных линейных операторов) является унимодулярной, если и только если $dim{{V}_{L}} = 0$, т.е. само пространство ${{V}_{L}}$ является нулевым.

Будем в дальнейшем придерживаться следующих стандартных обозначений. Кольцо $n \times n$-матриц ($n$ – положительное целое) с элементами в некотором кольце или поле $R$ обозначается через ${\text{Ma}}{{{\text{t}}}_{n}}(R)$. Если M – некоторая $n \times n$-матрица, то обозначение ${{M}_{{i,*}}}$, $1 \leqslant i \leqslant n$, используется для $1 \times n$-матрицы, равной i-й строке матрицы M. Диагональная $n \times n$-матрица с элементами ${{r}_{1}},\; \ldots ,\;{{r}_{n}}$ на диагонали обозначается через ${\text{diag}}({{r}_{1}},\; \ldots ,\;{{r}_{n}})$, единичная $n \times n$-матрица – через ${{I}_{n}}$.

II. ПРЕДВАРИТЕЛЬНЫЕ СВЕДЕНИЯ

Напомним, что разностное кольцо – это коммутативное кольцо $\mathbb{K}$ с единицей и с автоморфизмом $\sigma $ (который мы часто будем называть сдвигом). Если при этом $\mathbb{K}$ является полем, то оно называется соответственно разностным полем. В дальнейшем разностные поля без оговорок предполагаются полями характеристики 0.

Кольцом констант разностного кольца $\mathbb{K}$ называется ${\text{Const}}(\mathbb{K}) = \{ c \in K{\text{|}}\sigma c = c\} $. В случае, когда $\mathbb{K}$ является разностным полем, ${\text{Const}}(\mathbb{K})$ будет подполем поля $\mathbb{K}$ (полем констант поля $\mathbb{K}$).

Пусть $\mathbb{K}$ – разностное поле с автоморфизмом $\sigma $, кольцо $\Lambda $ – разностное расширение поля $\mathbb{K}$ (соответствующий автоморфизм кольца $\Lambda $ совпадает на $\mathbb{K}$ с $\sigma $, для этого автоморфизма кольца $\Lambda $ используется то же самое обозначение $\sigma $).

Определение 1. Кольцо $\Lambda $, являющееся разностным расширением поля $\mathbb{K}$, представляет собой адекватное разностное расширение поля $\mathbb{K}$, если ${\text{Const}}(\Lambda )$ является полем и для произвольной системы

с невырожденной матрицей $A \in {\text{Ma}}{{{\text{t}}}_{n}}(\mathbb{K})$ размерность линейного над полем ${\text{Const}}(\Lambda )$ пространства принадлежащих ${{\Lambda }^{n}}$ решений равна n.

Существование некоторого адекватного разностного расширения $\Lambda $ для произвольного разностного поля доказывается достаточно элементарно – см. [7, п. 5.1]. Для произвольного адекватного расширения равенство ${\text{Const}}(\Lambda ) = {\text{Const}}(\mathbb{K})$ не гарантируется, в общем случае ${\text{Const}}(\mathbb{K})$ является для ${\text{Const}}(\Lambda )$ собственным подполем.

Примечание 1. Так называемый q-разностный случай ([8, 9]) охватывается общим разностным случаем.

Скалярный разностный оператор является элементом кольца $\mathbb{K}[\sigma ,{{\sigma }^{{ - 1}}}]$.

Определение 2. Для ненулевого скалярного оператора $f = \sum {{{a}_{i}}{{\sigma }^{i}}} $ мы определяем его ведущий и трейлинговый порядки:

а также порядок ${\text{ord}}f\, = \,\overline {{\text{ord}}} f\, - \,\underline {{\text{ord}}} f.$ Полагаем $\overline {{\text{ord}}} 0$ = = –∞ , $\underline {{\text{ord}}} 0 = \infty $, ${\text{ord}}0 = - \infty $.

Для конечного набора $F$ скалярных операторов (вектора, матрицы, строки матрицы и т.д.) определим ведущий порядок $\overline {{\text{ord}}} F$ как максимум ведущих порядков его элементов, трейлинговый порядок $\underline {{\text{ord}}} F$ – как минимум трейлинговых порядков его элементов и полагаем ordF = $\overline {{\text{ord}}} F\, - \,\underline {{\text{ord}}} F$.

Разностная операторная матрица – это матрица, принадлежащая ${\text{Ma}}{{{\text{t}}}_{n}}(\mathbb{K}[\sigma ,{{\sigma }^{{ - 1}}}])$. Далее мы сопоставляем такой операторной матрице некоторые матрицы, принадлежащие ${\text{Ma}}{{{\text{t}}}_{n}}(\mathbb{K})$. Мы будем коротко называть операторную матрицу просто оператором.

Оператор имеет полный ранг (или: является оператором полного ранга), если его строки линейно независимы над $\mathbb{K}[\sigma ,{{\sigma }^{{ - 1}}}]$.

Если

и $L \ne 0$, то L можно записать в развернутом виде где ${{A}_{l}},{{A}_{{l - 1}}},\; \ldots ,\;{{A}_{t}} \in {\text{Ma}}{{{\text{t}}}_{n}}(\mathbb{K})$, при этом матрицы A_l, ${{A}_{t}}$ (ведущая и трейлинговая матрицы исходного оператора) ненулевые.

Определение 3. Пусть ведущие порядки строк оператора L равны ${{\alpha }_{1}},\; \ldots ,\;{{\alpha }_{n}}$, а трейлинговые – ${{\beta }_{1}},\; \ldots ,\;{{\beta }_{n}}$. Фронтальной матрицей оператора L назовем ведущую матрицу оператора PL, где

Соответственно тыльной матрицей оператора $L$ назовем трейлинговую матрицу оператора QL, где

Мы говорим, что оператор L строго редуцирован, если невырождены как фронтальная, так и тыльная матрицы этого оператора.

Определение 4. Оператор $L \in {\text{Ma}}{{{\text{t}}}_{n}}(\mathbb{K}[\sigma ,{{\sigma }^{{ - 1}}}])$ называется унимодулярным или обратимым, если для него существует обратный ${{L}^{{ - 1}}} \in {\text{Ma}}{{{\text{t}}}_{n}}(\mathbb{K}$[σ, σ^–1]): $L{{L}^{{ - 1}}} = {{L}^{{ - 1}}}L$ = I_n. Множество унимодулярных $n \times n$-операторов будет обозначаться через ${{\Upsilon }_{n}}$. Два оператора ${{L}_{1}},{{L}_{2}} \in {\text{Ma}}{{{\text{t}}}_{n}}(\mathbb{K}[\sigma ,{{\sigma }^{{ - 1}}}])$ назовем эквивалентными, если ${{L}_{1}} = U{{L}_{2}}$ для некоторого $U \in {{\Upsilon }_{n}}$.

Далее мы обозначаем через $\Lambda $ некоторое фиксированное адекватное разностное расширение исходного разностного поля $\mathbb{K}$ с автоморфизмом σ. Обозначение ${{V}_{L}}$ привлекается для пространства принадлежащих ${{\Lambda }^{n}}$ решений системы $L(y)$ = 0. Для краткости будем иногда говорить о ${{V}_{L}}$ как о пространстве решений оператора L.

Теорема 1. ([10]) Пусть $L \in {\text{Ma}}{{{\text{t}}}_{n}}(\mathbb{K}[\sigma ,{{\sigma }^{{ - 1}}}])$ имеет полный ранг. Тогда

(i) Если оператор $L$ строго редуцирован, то dimV_L = = $\sum\nolimits_{i = 1}^n \,{\text{ord}}{{L}_{{i,*}}}.$

(ii) $L \in {{\Upsilon }_{n}} \Leftrightarrow {{V}_{L}} = 0$.

III. РЕГУЛЯРИЗАЦИЯ: СЕМЕЙСТВА АЛГОРИТМОВ EG И RR

IV. ПРОВЕРКА УНИМОДУЛЯРНОСТИ, ОБРАЩЕНИЕ ОПЕРАТОРОВ

В основе имеющихся алгоритмов проверки унимодулярности оператора L и построения обратного оператора ${{L}^{{ - 1}}}$ в случае его существования, лежит вычисление $dim{{V}_{L}}$. Такое вычисление, в свою очередь, основывается на алгоритмах регуляризации (см. п. III и примечание 3).

Эквивалентные оператору L операторы

в силу примечания 3 являются строго редуцированным. По теореме 1 мы можем вычислить ${\text{dim}}{{V}_{L}}$ и установить, является ли L унимодулярным. Если $dim{{V}_{L}} = 0$, то в силу теоремы 1(i) соответствующий оператор из (5) принадлежит ${\text{Ma}}{{{\text{t}}}_{n}}(\mathbb{K})$, и ${{L}^{{ - 1}}}$ вычисляется легко.

Алгоритмы из [10] разработаны по образцу алгоритмов для дифференциального случая [5] с использованием RR, в то время как алгоритм из [2] построен на $\Delta {\text{EG}}$, это позволило снизить сдвиговую сложность – см. примечание 4. В [2] предложены алгоритмы Unimodularity Testing и Inverse Operator, сложности которых (арифметическая и сдвиговая) суть (3), (4) и соответственно

Построение обратного оператора выполняется с привлечением расширенных версий основных алгоритмов, в то время как проверка унимодулярности обходится основными версиями этих же алгоритмов.

Как показано в [10], основанные на RR алгоритмы имеют сложности $\Theta ({{n}^{4}}{{d}^{2}})$ и $\Theta ({{n}^{4}}{{d}^{3}})$ (если сохранять результаты всех сдвигов, иначе сдвиговая сложность будет $\Theta ({{n}^{4}}{{d}^{5}})$. Для алгоритмов, основанных на EG, сложности равны $\Theta ({{n}^{4}}{{d}^{2}})$, $\Theta ({{n}^{4}}{{d}^{2}})$).

V. РЕАЛИЗАЦИЯ, ЭКСПЕРИМЕНТАЛЬНОЕ СРАВНЕНИЕ

Реализация рассматриваемых алгоритмов²² выполнена в среде Maple [17] (предварительная версия реализации представлена в [2]). За основу взята описанная в [1] реализация алгоритмов EG, ΔEG, RR для дифференциального случая. Процедуры переработаны для использования в разностном случае, и дополнены реализацией расширенных версий алгоритмов. Также реализован алгоритм проверки унимодулярности оператора и построения обратного оператора, использующий в качестве базового алгоритма один из алгоритмов EG, ΔEG, RR по выбору пользователя. В этой реализации сдвиг определен как $\sigma y(x) = y(x - 1)$ и требуется специальный выбор заменяемого уравнения при выполнении шага исключения в алгоритмах EG, ΔEG для гарантии завершения работы. Существует вариант этих алгоритмов, который не требует такого специального выбора заменяемого уравнения. В этом варианте для шага сдвига используется оператор $\sigma - 1$ (использование аналогичного варианта алгоритма EG для q-разностного случая упомянуто в [18, примеч. 3]). Но при использовании алгоритмов EG, ΔEG для реализации алгоритма проверки унимодулярности оператора и поиска обратного оператора существенно использование именно версии с контролем порядка исключений (см. примечание 3).

Алгоритмы реализованы в виде процедур нового пакета EGRRext. Пакет предоставляет процедуры:

• EG – реализует алгоритм EG и его расширенную версию ExtEG;

• RR – реализует алгоритм RR и его расширенную версию ExtRR;

• TriangleEG – реализует алгоритм ΔEG и его расширенную версию ExtΔEG;

• IsUnimodular – реализует алгоритмы Unimodularity Testing и Inverse Operator.

Оператор $L = {{A}_{l}}{{\sigma }^{l}} + {{A}_{{l - 1}}}{{\sigma }^{{l - 1}}} + \cdots + {{A}_{t}}{{\sigma }^{t}}$ во входных параметрах процедур представляется в виде списка

где A – явная матрицаразмера $n \times n(l - t + 1)$. Явная матрица A представляется в виде стандартного Maple-объекта Matrix. Элементами явной матрицы являются рациональные функции одной переменной, которые в свою очередь представляются стандартным образом. Если t = 0, то оператор может быть также представлен только одной явной матрицей A.

Процедура IsUnimodular возвращает true или false в качестве результата проверки унимодулярности данного на вход оператора. Если указан необязательный входной параметр процедуры – имя переменной, то также вычисляется соответствующий обратный оператор в случае его существования, и вычисленный обратный оператор присваивается этой переменной. Обратный оператор также представляется в виде списка из его явной матрицы и его ведущего и трейлингового порядков. Если необязательное имя переменной не задано, то процедура использует алгоритм Unimodularity Testing, иначе – алгоритм Inverse Operator (см. п. IV). При этом процедура имеет еще один необязательный параметр method, который может принимать значения EG, TEG или RR и позволяет указать какой из алгоритмов EG, ΔEG, RR, соответственно, использовать в качестве базового (если этот параметр не указан, то используется EG). В первой версии реализации, представленной в [2], в качестве базового алгоритма использовался только EG.

Пример 1. Рассмотрим применение процедуры IsUnimodular к оператору

Явная матрица оператора равна

с $l = 1$ и t = 0. Процедура вызывается дважды для каждого из методов: первый раз только для проверки унимодулярности, а во второй раз и для построения обратного оператора. Дополнительно выведено и время вычислений:

> L := Matrix([[0, –1/x, 1, 0],

          [0, –x/2, x^2/2, 1]]);

> st:=time():

  IsUnimodular(L, x, ’method’=’EG’);

  time()-st;

> st:=time():

  IsUnimodular(L, x, ’InvL_EG',

                         ’method’=’EG’);

  time()-st;

> InvL_EG;

> st:=time():

  IsUnimodular(L, x, ’method’=’TEG’);

  time()-st;

> st:=time():

  IsUnimodular(L, x, ’InvL_TEG',

             ’method’=’TEG’);

  time()-st;

> InvL_TEG

> st:=time():

 IsUnimodular(L, x, ’method’=’RR’);

 time()-st;

> st:=time():

  IsUnimodular(L, x, ’InvL_RR',

                       ’method’=’RR’);

  time()-st;

> InvL_RR;

Таким образом

Пример 2. Оператор

где 0_n – нулевая матрица размера $n \times n$, M₁, ${{M}_{2}} \in {\text{Ma}}{{{\text{t}}}_{n}}(\mathbb{K}[\sigma ,{{\sigma }^{{ - 1}}}])$ – произвольные операторы, является унимодулярным для любых ${{M}_{1}}$ and ${{M}_{2}}$. Обратный оператор равен

Выполнена серия экспериментов. Для каждого из них генерировалась пара из разреженных на 50% операторов ${{M}_{1}}$ и ${{M}_{2}}$ с элементами в виде случайных полиномов степени не больше 2. Затем для (2) находился обратный оператор. Операторы ${{M}_{1}}$ и ${{M}_{2}}$ имели одинаковое число строк, оператор ${{M}_{1}}$ имел порядок 2, а порядок оператора ${{M}_{2}}$ принимал значения $d = 3,5,7,9,11$ (тот же порядок имел и оператор M). Число строк ${{M}_{1}}$ и ${{M}_{2}}$ равнялось $n = 2,3,4$ (соответственно, число строк M равнялось $k = 4,6,8$). Обратный оператор M в каждом эксперименте вычислялся всеми тремя методами. Полученные результаты представлены в табл. 1. Использование метода, основанного на RR, в ряде экспериментов не дало возможности получить результат – вычисления были принудительно остановлены, когда был исчерпан доступный объем физической памяти используемого компьютера. Это помечено в таблице знаком $ > $ с указанием времени до принудительной остановки вычислений.

VI. АНАЛИЗ РАСХОЖДЕНИЙ

Из табл. 1 видно, что результаты в целом соответствуют теоретическим оценкам сложности. При этом относительный рост времени вычислений с ростом числа строк и порядка оператора происходит заметно быстрее, чем это можно было ожидать, исходя из этих теоретических оценок. Чтобы понять причину этих расхождений нами был реализован режим работы процедуры IsUnimodular, в котором дополнительно выводятся характеристические параметры промежуточных результатов вычислений. Эти параметры выводятся после каждого появления нулевой строки в выявляющей матрице (после выполнения редукции в алгоритмах EG и RR или после выполнения одной или нескольких замен строк в алгоритме ΔΕΓ). В качестве таких параметров используются параметры строки оператора, в которой появилась нулевая строка в выявляющей матрице, а именно объем данных в представлении этой строки: используется стандартная Maple-функция length, а также максимум сумм степеней числителя и знаменателя коэффициентов в этой строке оператора. Эти дополнительные данные показали, что в ходе промежуточных вычислений коэффициенты оператора значительно “распухают”. Несмотря на то, что в эксперименте суммы степеней числителей и знаменателей коэффициентов в исходном операторе M не превышает 2, а в обратном к нему не превышает 4, в промежуточных вычислениях эта величина достигала сотен и даже тысяч. Столь же значительно рос и объем данных.

В табл. 2 представлены сводные данные по тем же экспериментам, которые представлены в 1. Для каждого эксперимента представлено пять параметров:

• максимальный среди всех промежуточных строк объем используемых данных,

• максимальная среди всех промежуточных строк величина максимума суммы степеней числителя и знаменателя,

• число промежуточных строк,

• средний по всем промежуточным строкам объем используемых данных,

• средняя по всем промежуточным строкам величина максимума суммы степеней числителя и знаменателя.

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

Представленный в табл. 1 значительный рост времени вычислений в экспериментах с ростом числа строк и порядка оператора согласуется с “распуханием” коэффициентов в промежуточных вычислениях, представленных в табл. 2 характеристическими параметрами. В результате такого “распухания” фактические вычислительные затраты на выполнение одной арифметической операции, как и операции сдвига, начинают существенно расти, и этот рост оказывается более существенным в сравнении с ростом числа самих таких операций, отраженных в сложности алгоритмов. Для разных методов, это “распухание” происходит с разной скоростью. В наших экспериментах с этой точки зрения особенно плохо проявил себя алгоритм RR. Это согласуется с результатами сравнения аналогичных алгоритмов в дифференциальном случае, представленными в [1]. Отметим также, что в условиях, когда “распухание” в разных алгоритмах имеет схожие характеристики, то есть затраты на выполнения одной операции в каждом алгоритме сопоставимы, то на первый план снова выходит число этих операций, и относительное время вычислений больше соответствует теоретическим оценкам сложности.

И отметим еще одну особенность алгоритмов EG и RR. Для выполнения редукции в этих алгоритмах используется поиск линейных зависимостей в выявляющих матрицах. В реализации для этого используется процедура NullSpace пакета LinearAlgebra системы Maple. Эта процедура может найти более чем одну линейную зависимость, что дает возможность различного продолжения вычислений, в зависимости от выбора для исключений той или иной из найденных линейных зависимостей. В нашей реализации выбиралась первая из найденных линейных зависимостей после их сортировки с помощью процедуры ComplexitySort пакета SolveTools. В табл. 3, 4 представлены результаты вычислений для части тех же тестовых систем, что и в табл. 3, 4, табл. 1, 2, но с использованием выбора не первой, а последней линейной зависимости после сортировки. Представленные результаты показывают, что порядок исключений также существенно влияет на “распухание” коэффициентов в промежуточных вычислениях, а значит и на общее время вычислений. Этот фактор также не принимается во внимание в оценках сложности, однако с практической точки зрения является весьма существенным.

Исходя из проведенных экспериментов, а также из ранее проведенных экспериментов для дифференциального случая ([1]), с практической точки зрения разумно использовать версию алгоритмов Unimodularity Testing и Inverse Operator, основанных на алгоритме EG. Именно поэтому такой вариант реализации был представлен в [2] и используется по умолчанию (если не указан параметер method) в реализации в этой работе. Эти эксперименты также показывают, что с практической точки зрения полезно иметь реализации различных алгоритмов. Для конкретного примера более эффективным может оказаться алгоритм, уступающий другим в сложностном отношении или времени счета в других примерах. Поэтому в нашей реализации оставлена возможность использования и других алгоритмов по выбору пользователя.

БЛАГОДАРНОСТИ

Работа выполнена при частичной поддержке РФФИ, грант № 19-01-00032-a.

Авторы благодарны А.А. Рябенко за замечания по первому варианту статьи и полезные советы.