Введение В связи с активным развитием технологий параллельного программирования приобретают актуальность методы вычислений, в основе которых лежат непозиционные системы счисления. Наибольшее распространение среди таких систем получили системы остаточных классов [1, 2]. Система остаточных классов (СОК, модулярная система) определяется набором взаимно простых модулей (оснований) p1, p2, …, pn. Целое число X из интервала [0, P - 1], где P - произведение всех модулей, многоразрядное в позиционной системе, представляется в СОК в виде нескольких малоразрядных остатков (модулярных разрядов): , (1) где для всех i = 1, 2, …, n. Использование СОК - перспективный способ организации высокоскоростной арифметической обработки чисел большой разрядности. Все модульные операции над модулярными числами вида (1), такие как сложение, умножение, возведение в степень и пр., выполняются отдельно и независимо по каждому модулю, следовательно, легко и быстро в связи с малой разрядностью остатков xi. Кроме этого, структуры данных, основанные на принципах модулярного представления, эффективно векторизуются, а также могут быть распределены по ядрам многоядерных процессоров, что обусловливает возможность распараллеливания вычислений на уровне арифметических операций. Это особенно актуально в связи с активным развитием высокопроизводительных вычислительных систем. Однако наибольшей проблемой на пути создания быстродействующих алгоритмов модулярной обработки остается высокая сложность немодульных операций. К таким операциям относятся, в частности, сравнение чисел по величине, вычисление знака, оценка выхода результата операции за допустимый диапазон представления, масштабирование и пр. Для выполнения этих операций, в отличие от модульных, недостаточно знания отдельных остатков xi - требуется оценка значения числа X в целом, которая затруднена в связи с непозиционной природой СОК. В связи с вышеизложенным нами рассматривается новый алгоритм сравнения чисел в СОК, основанный на использовании интервально-позиционных характеристик. Способы интерпретации позиционного значения модулярного кода Для сравнения чисел в СОК необходимо обладать информацией об их позиционной величине. Величина модулярного числа (1) определяется соотношением , (2) где B1, B2, …, Bn - ортогональные базисы СОК [1], каждый i-й из которых суть произведение чисел и . Здесь - это вес ортогонального базиса (мультипликативная инверсия от Pi по модулю pi). Величина ортогональных базисов пропорциональна произведению модулей, из-за чего прямое вычисление суммы (2) оказывается трудоемким. Вследствие этого для выполнения немодульных операций используются различные функции от остатков, позволяющие оценить величину модулярного числа без ее непосредственного вычисления. Обобщенно такие функции принято называть позиционными характеристиками (ПХ). К числу наиболее известных ПХ относятся ранг, след, функция ядра [1, 2], диагональная функция [3], четность [4], представление в системе счисления со смешанными основаниями (Mixed Radix System) [5] и пр. Однако алгоритмы вычисления этих ПХ также отличаются высокой сложностью. Применение методов табличной выборки, в свою очередь, приводит к неприемлемым затратам памяти при работе в больших числовых диапазонах. В связи с этим наиболее перспективными представляются подходы, основанные на оценке относительной величины модулярных чисел. Определение. Относительная величина E(X/P) модулярного числа X - это отношение его позиционного значения к произведению P всех модулей СОК: (3) Точное значение E(X/P) в общем случае не представимо в ЭВМ с ограниченной разрядной сеткой, поэтому на практике вместо нее вычисляют приближенную позиционную характеристику (ППХ) - округленное значение относительной величины [2, с. 194; 6-8]. Вычисление ППХ выполняется на порядок быстрее, чем точных ПХ (за линейное и логарифмическое время по числу модулей соответственно при последовательной и параллельной реализации алгоритма). Однако из-за сложности, а часто - принципиальной невозможности учета всех факторов, влияющих на величину погрешности ППХ, вызванной округлениями (в особенности в условиях неоднозначности модели вычислений с плавающей точкой, реализуемой современными процессорами и компиляторами языков программирования), при ее использовании не всегда гарантируется корректность результата. Применительно к операции сравнения это означает, что не исключается возможность возникновения ошибок вычислений, связанных с неправильным выбором ветви при условном переходе, и особенно актуально, когда базис СОК задан большим числом модулей (в этом случае погрешности вычисления ППХ контролировать наиболее сложно). Таким образом, актуальна задача разработки новых алгоритмов выполнения немодульных операций в СОК, и в частности сравнения, которые обладают вычислительной сложностью, соизмеримой со сложностью вычисления ППХ, но при этом имеют строгое обоснование и обеспечивают математическую достоверность результата. Интервально-позиционная характеристика Рассмотрим метод достоверной оценки позиционной величины модулярных чисел, в основе которого лежит интервальная аппроксимация относительной величины (3). Центральным в нем является понятие интервально-позиционной характеристики модулярного числа [9]. Определение. Интервально-позиционная характеристика (ИПХ) - это замкнутый вещественный интервал c направленно округленными границами, которые удовлетворяют условию включения: . Интервально-позиционная характеристика проецирует диапазон СОК на полуинтервал [0, 1), ассоциируя всякое модулярное число X с парой округленных позиционных чисел - границ, которые при некоторых вполне определенных условиях локализуют его относительную величину (рис. 1). Рис. 1. Интервально-позиционная характеристика модулярного числа Для вычисления границ ИПХ справедливы формулы (4) где стрелки соответствуют направленным округлениям при вычислении и суммировании слагаемых: ¯ - округление с недостатком («вниз»), - округление с избытком («вверх»). В последовательном случае для вычисления этих выражений требуется O(n) элементарных операций, в параллельном - O(log n). Для сравнения: преобразование в систему со смешанными основаниями (Mixed Radix Conversion, MRC) требует соответственно O(n2) и O(n) операций [5, 10]. Интервально-позиционная характеристика является позиционной аппроксимацией относительной величины (3), поэтому может быть использована для выполнения таких немодульных операций в СОК, как сравнение, определение знака числа, контроль четности, определение переполнения суммы или произведения чисел, масштабирование и пр. Остановимся подробнее на операции сравнения, поскольку она является основной в классе немодульных процедур, а оставшиеся из перечисленных могут быть реализованы на ее основе. Алгоритм модулярного сравнения Пусть в СОК заданы числа X и Y. Тогда для их сравнения в нетривиальном случае (X ¹ Y) необходимо вычислить ИПХ I(X/P) и I(Y/P), проверить формальные условия корректности результата и, если они выполняются, определить результат, сравнив ИПХ по правилам интервальной арифметики. Нарушение формальных условий корректности результата свидетельствует о недостаточной точности вычисления ИПХ. Рассмотрим эти условия. Абсолютная ошибка ИПХ отражается в ее диаметре: (5) За счет направленных округлений погрешности, возникающие при вычислении ИПХ, приводят лишь к увеличению диаметра, не оказывая в общем случае влияния на свойство включения E(X/P) Î I(X/P). Однако, поскольку область значений границ ИПХ ограничена полуинтервалом [0, 1), в ряде случаев, из-за недостаточной точности их вычисления, указанное свойство может нарушаться. Это происходит тогда, когда число X очень мало́ по отношению к P либо, наоборот, находится в непосредственной близости с точкой P - 1. В первом случае неправильно вычисляется нижняя граница ИПХ, во втором - верхняя. В любом случае diam I(X/P) < 0, т. е. нижняя граница больше верхней. Такая ИПХ называется неправильной по Каухеру или просто неправильной. Первое формальное условие корректного сравнения - правильность ИПХ операндов. Если это условие выполняется, то окончательный вывод о достоверности результата формулируется на основании проверки второго формального условия. Применительно к операции сравнения этим условием является отсутствие пересечения (коллизии) ИПХ, которое в терминах интервального исчисления определяется интервалом Если диаметр (5) этого интервала меньше нуля, то ИПХ не пересекаются в стандартном теоретико-множественном смысле, т. е. не содержат общих точек. Конъюнкция первого и второго условий образует признак (достаточное условие) корректности результата операции сравнения. Выполнение достаточного условия свидетельствует об однозначной интерпретации результата при текущей точности вычисления ИПХ. Напротив, нарушение достаточного условия свидетельствует о слишком малой точности ИПХ, по крайней мере для одного из операндов. В этом случае появляется возможность повысить их точность или, если это неосуществимо в рамках машинной арифметики (без применения многоразрядных форматов), использовать другой метод сравнения, например преобразование в систему со смешанными основаниями с последующим сопоставлением полиадических кодов (MRC-метод). С учетом сказанного алгоритм модулярного сравнения формулируется следующим образом. Алгоритм (сравниваются числа и , представленные в СОК с модулями p1, p2, …, pn). Шаг 1. Выполнить попарное сравнение остатков для исключения тривиального случая: если xi = yi для всех i = 1, 2, …, n, то X = Y и алгоритм завершается. Шаг 2. Вычислить ИПХ I(X/P) и I(Y/P). Шаг 3. Проверить первое формальное условие корректного сравнения: если и , то условие выполняется. В этом случае перейти к шагу 4, иначе - к шагу 6. Шаг 4. Если , то и E(X/P) > E(Y/P), поскольку E(X/P) Î I(X/P) и E(Y/P) Î I(Y/P). В этом случае завершить алгоритм с результатом X > Y. Иначе перейти к шагу 5. Шаг 5. Если , то завершить алгоритм с результатом X < Y. Невыполнение этого неравенства на данном шаге говорит о нарушении второго формального условия корректного сравнения из-за слишком малой точности вычисления ИПХ. При этом следует перейти к шагу 6. Шаг 6. Если повышение точности вычисления ИПХ неосуществимо в рамках машинной арифметики, то преобразовать числа X и Y из СОК в систему со смешанными основаниями и сравнить цифры полиадических кодов, начиная со старшей (использовать MRC-метод). Пример. В СОК с модулями {7, 9, 11, 13} и весами ортогональных базисов {6, 5, 9, 10} даны числа и . Требуется выполнить их сравнение. 1. Числа, очевидно, не равны, поэтому вычисляем их ИПХ в соответствии с формулами (4) с округлением до двух цифр: 2. Интервалы I(X/P) = [0,01, 0,04] и I(Y/P) = [0,02, 0,06] правильные (первое условие выполнено), но пересекаются, т. к. больше нуля. Результат сравнения не может быть однозначно определен при текущей точности вычисления ИПХ. 3. Повысим точность на один разряд: 4. В данном случае , поэтому второе условие корректного сравнения также выполнено. Поскольку , то делаем вывод, что X < Y. Проверка преобразованием в десятичную систему: X = 270, Y = 310. Оценка трудоемкости алгоритма Трудоемкость представленного алгоритма сравнения зависит от сочетания входных данных: при неудачной их конфигурации и недостаточной точности ИПХ могут нарушаться условия достоверности результата. В этом случае будем полагать, что для окончательного его определения используется MRC-метод. В соответствии с этим в приведенных ниже оценках используется метрика ρMRC - вероятность использования MRC-метода (перехода к шагу 6), которая является композицией из вероятностей нарушения первого и второго условий. Будем считать также, что указанные события возникают с равной вероятностью, которая, как показывает анализ графа возможных состояний алгоритма, равна . Определим прежде всего трудоемкость MRC-метода. Для вычисления всех цифр смешанного представления с использованием последовательного алгоритма [5] требуется n(n - 1)/2 арифметических вычитаний и столько же умножений по модулю. Кроме этого, n операций потребуется в худшем случае для сравнения цифр смешанного представления. Следовательно, общая трудоемкость MRC-метода сравнения чисел в СОК (6) С учетом оценки (6) трудоемкость представленного алгоритма сравнения на основе ИПХ, в котором MRC-преобразование используется с вероятностью ρMRC, составляет (7) где n операций необходимы для проверки полного совпадения чисел и еще около 10n операций - для вычисления двух ИПХ [9]. Оценка (7) не учитывает операции, необходимые для переключения режимов округления процессора, а также операции позиционного сравнения границ ИПХ, которые при большом числе модулей n не вносят сколь-либо существенного вклада в общую трудоемкость. Таким образом, предлагаемый алгоритм сравнения чисел в СОК обеспечивает достоверность результата и, в зависимости от сочетания входных данных, позволяет до 0,18n раз (при ρMRC = 0) повысить скорость вычислений по отношению к MRC-алгоритму. На рис. 2 представлен график оценочного ускорения алгоритма для различных сочетаний n и ρMRC. Точки ρ1, ρ2 и ρ3 соответствуют значениям вероятности ρMRC, при которых предлагаемый алгоритм по крайней мере не медленнее MRC-аналога при n = 16, n = 32 и n = 64 соответственно. Рис. 2. Зависимость оценочного ускорения предлагаемого алгоритма модулярного сравнения относительно MRC-алгоритма Сходные оценки могут быть получены и для параллельных алгоритмов. В этом случае вычисление ИПХ требует O(log n) элементарных операций [9], а преобразование к системе счисления со смешанными основаниями выполняется за линейное время при использовании техники, представленной в [10]. К вопросу о точности вычисления интервально-позиционной характеристики Формулы (4), при присущей им низкой вычислительной сложности, обладают недостатками, связанными с невысокой точностью при вычислении ИПХ для малых модулярных чисел. Рассмотрим их. В качестве меры точности ИПХ используем ее относительную погрешность: (8) где и определяются отношениями абсолютных погрешностей границ к точной относительной величине (3). При большой погрешности (8) оценка относительной величины становится неинформативной, а условия корректности результата могут нарушаться, что приводит к росту времени вычислений. Пусть ε - предел допустимой относительной ошибки ИПХ, а ее границы вычисляются по формулам (4) в k-значной арифметике с плавающей точкой (ПТ), т. е. разрядность их мантисс равна k. Показано [11], что , только если число X удовлетворяет неравенству (9) где ψ = n2-k/e, n - число модулей СОК, а P - их произведение. Если X < ψP, то формулы (4) не обеспечивают заданную точность. Уменьшение ψ за счет увеличения разрядности k влечет необходимость использования длинной арифметики при большом P, что приводит к снижению быстродействия. Таким образом, при вычислении ИПХ выдвигаются два требования - точность и скорость. Рассмотрим алгоритм, позволяющий достичь компромисса в их удовлетворении. Условимся прежде всего, что границы ИПХ представлены числами с ПТ: В этом случае их деление на натуральные степени двойки выполняется быстро и безошибочно (без модификации мантиссы, при условии, что результат лежит в пределах нормализованного диапазона). Это позволяет, когда X не отвечает условиям (9), вычислить смещенную ИПХ - характеристику числа , для которого условия (9) выполняются, а затем без увеличения погрешности перейти к искомой ИПХ масштабированием смещенных границ коэффициентом 2v. Однако при таком подходе на смещение накладываются ограничения: для всех (10) которые, как легко убедиться, не являются взаимно согласованными, если P ³ 1/2ψ2. Это противоречие устраняется декомпозицией критической области [1, ψP]. При этом вместо скалярного определяется векторное смещение: где для всех j = 1, 2, …, g, причем . Такой способ задания последовательности показателей {vj} требует, чтобы число ψ было строго меньше 0,5. Элементы вектора v разбивают отрезок [1, ψP] на g интервалов и обеспечивают выполнение условий (10) в каждом из них по отдельности. Для ускорения расчетов предварительно вычисляется подстановочная таблица - матрица смещенных весов ортогональных базисов, строки которой ассоциированы с элементами вектора v, а столбцы - с модулями СОК: С учетом сказанного алгоритм ISaC (Iterative Shift and Correction) высокоточного вычисления ИПХ формулируется следующим образом. Алгоритм (для числа в k-значной арифметике с ПТ вычисляется I(X/P) с погрешностью , не превышающей допустимого предела e). Шаг 1. По второй формуле в (4) вычислить значение верхней границы . Шаг 2. Если , то, при допущении о пренебрежительной малости погрешности верхней границы относительно числа ψ, условие (9) выполняется, поэтому уточненное вычисление ИПХ не требуется. В этом случае перейти к шагу 3, иначе - к шагу 4. Шаг 3. По первой формуле в (4) вычислить и перейти к завершающему шагу 10. Шаг 4. Принять j = 1 и перейти к следующему шагу. Шаг 5. Вычислить смещенную верхнюю границу по формуле (11) где Mj,i - i-й элемент j-й строки матрицы М. Шаг 6. Если , то заданная точность вычислений достигнута. В этом случае выйти из тела цикла, перейдя к шагу 8. Иначе перейти к шагу 7. Шаг 7. Инкрементировать индекс смещения j и вернуться к шагу 5. Шаг 8. По найденному в итерационном блоке индексу j выбрать j-ю строку матрицы М и вычислить смещенную нижнюю границу по аналогии с формулой (11), но с округлением «вниз» каждой дроби и их суммы. Шаг 9. По индексу смещения j выбрать из вектора v элемент и масштабировать им смещенные границы, получив границы искомой ИПХ: Шаг 10. Закончить выполнение алгоритма, искомым результатом являются вычисленные значения границ и . Схема алгоритма представлена на рис. 3. Конечность алгоритма обеспечивается способом задания вектора v, который обеспечивает полное перекрытие критической области, т. е. гарантирует, что , где . Максимальное количество итераций уточняющего цикла равно размеру g вектора v. Технический предел алгоритма (кроме естественного условия ψ < 0,5) - ограничение минимального показателя степени (порядка) emin в формате с ПТ, используемом для представления границ ИПХ: необходимо, чтобы разность e - vj была не меньше emin, где e - порядок в экспоненциальном представлении границы. В частности, emin = −1022 для формата двойной точности спецификации IEEE-754 и emin = −16382 для формата двойной расширенной точности. Рис. 3. Высокоточный алгоритм вычисления интервально-позиционной характеристики При использовании алгоритма отпадает необходимость в применении многоразрядной арифметики, а вычисление сумм в формулах смещенных границ может быть выполнено с распараллеливанием по данным за логарифмическое время. Выполненные экспериментальные исследования, отраженные в [11], показывают, что при работе в 32-модульной СОК быстродействие алгоритма ISaC при одинаковой точности вычисления ИПХ более чем в 11 раз выше классического многоразрядного аналога (реализующего формулы (4) с использованием 490-битной арифметики библиотеки MPFR). Произведенная оптимизация программной реализации позволила улучшить данные показатели по меньшей мере на 18 %. Заключение 1. Разработанный алгоритм сравнения чисел в СОК прост и не предполагает использование многоразрядной арифметики. Его быстродействие, в зависимости от сочетания входных данных, выше до 0,18n раз по сравнению с алгоритмом на основе преобразования модулярного кода в систему счисления со смешанными основаниями. При распараллеливании вычислений в неисключительном случае (при выполнении всех условий корректности результата) алгоритм характеризуется логарифмическим временем работы по числу модулей СОК. По аналогии с представленным алгоритмом сравнения формулируются также и другие алгоритмы выполнения немодульных операций, в частности вычисления знака числа, представленного в симметричной СОК, контроль четности, оценка переполнения диапазона при сложении и умножении. 2. При использовании ИПХ для оценки относительной величины числа в модулярном представлении осуществляется естественный учет погрешностей округления, не требующий рассмотрения специфики используемой модели машинных вычислений. Это позволяет простым образом контролировать корректность результата немодульной операции, вне зависимости от количества модулей СОК и их разрядности. 3. Разработан новый алгоритм вычисления ИПХ, который позволяет получить высокоточную информацию о величине числа в модулярном представлении без использования трудоемкого преобразования в позиционную систему счисления. Алгоритм основан на быстром итерационном смещении вычисляемых границ с их последующей безошибочной корректировкой масштабированием степенью двойки и не требует работы с длинными числами, что, при равной точности, обусловливает его высокое быстродействие по сравнению с классическим аналогом на основе многоразрядной арифметики. 4. Полученные теоретические и практические результаты могут быть использованы при реализации высокоточных вычислений в модулярной арифметике, когда базис СОК задан большим числом модулей, а применение для выполнения немодульных операций MRC-, CRT-методов или табличной выборки сопряжено с большими временными или емкостными затратами соответственно.