CC BY
136
ISSN 2079-3316 ПРОГРАММНЫЕ СИСТЕМЫ: ТЕОРИЯ И ПРИЛОЖЕНИЯ №1(28), 2016, с. 61-97 УДК 004.222.3+519.6
К. С. Исупов, В. С. Князьков
Арифметика многократной точности на основе систем остаточных классов
Аннотация. В работе рассматриваются алгоритмы высокоточной арифметики, основанные на использовании многомодульных систем остаточных классов (СОК) для представления мантисс чисел с плавающей точкой произвольной разрядности; порядок хранится в двоичной системе счисления. Такое представление обеспечивает большой динамический диапазон и допускает эффективное распараллеливание арифметических операций над цифрами многоразрядных мантисс по модулям СОК, что хорошо согласуется с особенностями архитектуры современных параллельных вычислительных систем. Дополнительно, в числовой формат включена атрибутивная информация, которая обеспечивает быструю оценку величины мантиссы, масштабированной относительно произведения модулей, и способствует ускоренному выполнению целого ряда проблемных для СОК немодульных процедур, таких как сравнение, контроль переполнения диапазона, округление и пр. Представлены результаты экспериментальной оценки точности и быстродействия алгоритмов, а также эффективности использования SIMD.
Ключевые слова и фразы: компьютерная арифметика, высокоточные вычисления, параллельные алгоритмы, система остаточных классов, SIMD.
Введение
С развитием параллельных аппаратных платформ и ростом масштабов вычислений одним из важнейших вопросов, требующих решения, становится обеспечение приемлемой точности арифметики с плавающей точкой. В ближайшее время суперкомпьютеры, вероятно, достигнут экзафлопсной производительности, и стоит ожидать, что получение результатов необходимой точности на таких компьютерах в машинной арифметике IEEE 754 будет непростой задачей [1]. В связи с этим в настоящее время активно развиваются методы и программные средства, обеспечивающие достоверность и воспроизводимость результатов вычислений на параллельных системах [2,3].
Исследование выполнено при финансовой поддержке РФФИ в рамках научных проектов № 14-07-31075 мол_а и № 16-37-60003 мол_а_дк.
© К. С. Исупов, В. С. Князьков, 2016 © Вятский государственный университет, 2016 © Программные системы: теория и приложения, 2016
Они основываются на алгоритмах обработки машинных чисел, позволяющих компенсировать погрешность округления. Однако возникает все больше различных задач, особенно в компьютерном моделировании, корректное решение которых может быть получено только посредством использования арифметики многократной точности, позволяющей оперировать числами произвольной разрядности [4].
Уже сейчас арифметика многократной точности востребована во многих областях, которые включают решение плохо обусловленных систем, продолжительное и / или крупномасштабное моделирование, исследование мелкомасштабных явлений, вычисление рядов, численное интегрирование и пр. [5,6]. Конкретно, арифметика многократной точности используется при изучении динамических систем (например, исследование фрактальных свойств аттрактора Лоренца [7]), для вычисления моментов произвольных порядков в хроматографии [8], при моделировании климата [9], для вычисления амплитуд рассеяния фундаментальных частиц [10], в задачах прямой кинематики [11] и строительной механики [12], в вычислительной геометрии и во многих других. С ростом производительности суперкомпьютеров и масштабов вычислений число приложений высокоточной арифметики существенно возрастает, и имеются основания полагать, что в будущем эта тенденция будет лишь усиливаться.
Большой размер задач, необходимость их быстрого решения и особенности современных архитектур (такие, как вложенный параллелизм) определяют новые требования к высокоточному программному обеспечению. К ним относятся, в первую очередь, высокая скорость, потоковая безопасность и, что важно, возможность распараллеливания многоразрядных операций. Однако классические методы длинной арифметики [13], которые лежат в основе большинства известных пакетов (GMP, MPFR, ARPREC, QD, NTL и пр.), приводят к довольно медленным и неэффективным реализациям. Фундаментальная причина этого — возникновение цепочек переносов, из-за которых алгоритмы обработки многоразрядных мантисс становятся вычислительно сложными и не распараллеливаются. В результате высокоточные расчеты сопровождаются большими затратами времени и неэффективным использованием вычислительных ресурсов.
В качестве перспективного подхода к решению обозначенной проблемы стоит отметить разработанные Н.Н. Непейводой позиционные системы счисления с перекрытием [14]. Представление чисел большой
разрядности в одной из таких систем счисления позволяет исключить необходимость распространения цепочек переносов и допускает эффективное распараллеливание вычислений.
В основе нашего подхода, рассматриваемого далее, лежит идея эффективной параллельной обработки многоразрядных мантисс за счет их представления в виде разложений по модулям в системах остаточных классов (СОК) [15-18], обладающих высоким уровнем естественного параллелизма арифметических операций.
Оставшаяся часть статьи построена следующим образом. В разделе 1 описано используемое представление чисел с плавающей точкой произвольной длины (MF-формат). В разделе 2 рассмотрен метод интервально-позиционных характеристик, позволяющий решить основные проблемы СОК, связанные с высокой сложностью немодульных вычислений. Раздел 3 содержит кодировки числовых и нечисловых величин, представимых в MF-формате. В разделе 4 определены границы числового диапазона, рассматриваются вопросы однозначного определения его переполнения. Раздел 5 содержит описание базовых алгоритмов арифметики многократной точности. В разделе 6 содержатся некоторые результаты оценки эффективности предложенного подхода. В разделе 7 определены направления дальнейших исследований, и, наконец, раздел 8 содержит краткие выводы.
1. Представление данных
Используется следующее представление чисел произвольной длины. Число с плавающей точкой представляется знаком s € {0,1}, мантиссой в системе остаточных классов X, несмещенным порядком (экспонентой) е и интервально-позиционной характеристикой (ИПХ) мантиссы I(Х/М). Порядок е представлен целым машинным числом со знаком. Мантисса X представлена кортежем остатков (х\,х2,... ,хп) по модулям СОК т2,..., тп. Каждый г-й остаток, в свою очередь, представляется неотрицательным целым машинным числом, удовлетворяющим сравнению Xi = X mod т^. Мантисса X интерпретируется как целое число без знака, изменяющееся в интервале от 0 до М — 1, причем М = П"=1 mi — динамический диапазон СОК. ИПХ мантиссы представляется двумя направленно округленными числами машинной точности — нижней границей Х/М и верхней границей Х/М — и локализует значение мантиссы X, масштабированное относительно произведения модулей М, т.е. Х/М < Х/М < Х/М. Описанное
числовое представление называется модулярно-позиционным форматом с плавающей точкой (сокращ. MF-формат от Modular-positional Floating-point, или МП-формат) [19,20]. Для компактной записи чисел в таком формате используется нотация х ^ [s,X,e,I(Х/М)}, которая выражает следующее числовое значение:
(1) х =(-1Г х
г= 1
xi|M-1|r
х 2е,
м
где Mi = M/mi = ТО1ТО2 • • • mj_imj+i • • • mn и |M— 1|mi — мультипликативная инверсия Mi по модулю то-j, представляющая собой такое целое неотрицательное число, что IM-1I„HMi = 1 mod то-j. Можно доказать существование мультипликативной инверсии для каждого Mi, при условии, что все модули TOj, i = 1, 2,..., п, попарно взаимно простые, т.е. gcd(то^, тоj) = 1 для всех i = j. Сумма по модулю М в соотношении (1) выражает десятичную величину мантиссы X и выводится из китайской теоремы об остатках [17].
MF-формат обладает информационной избыточностью: ИПХ является полем формата, но, как следует из (1), не участвует в образовании значения числа. Ее предназначение здесь — дать быструю и достоверную оценку величины мантиссы при выполнении немодульных операций, являющихся крайне проблемными в СОК (см. раздел 2). Включение ИПХ в числовое представление целесообразно, по крайней мере, по двум причинам:
• во-первых, это позволяет организовать эффективную схему вычислений с предварительной оценкой операндов (см. раздел 5.1);
• во-вторых, это дает возможность использовать интервальную арифметику при выполнении модульных операций (умножение, сложение, вычитание и пр.), позволяя рассчитать ИПХ результатной мантиссы за несколько тактов.
2. Интервальная оценка в СОК
Несомненным преимуществом систем остаточных классов, по сравнению с классическими системами счисления, является их непозиционная (невзвешенная) природа, благодаря которой для таких операций, как сложение, вычитание и умножение, большой динамический диапазон разбивается на поддиапазоны пониженной разрядности путем разделения потока вычислений на несколько независимых параллельных каналов (рис. 1).
Пусть даны числа X = {х\, Х2,..., хп) и У = у2, . . . , уп) (беззнаковые) в СОК с модулями т\, То2,..., то„. Тогда их сложение и умножение будут определяться в соответствии с алгоритмами 1 и 2.
_Алгоритм 1. Сложение чисел в СОК_
1: procedure RnsSum(x, Y) 2: for i ^ 1 to n do
3: Zi = \xi + Vi\mi > Все Zi вычисляются параллельно
4: end for
5: return Z = (z\,z2,..., zn) 6: end procedure
_Алгоритм 2. Умножение чисел в СОК_
1: procedure RnsProd(x, Y) 2: for i ^ 1 to n do
3: Zi = \xi x yi\mi > Все Zi вычисляются параллельно
4: end for
5: return Z = (z\,z2,..., zn) 6: end procedure
Представленные алгоритмы полностью параллельны и эффективны при реализации с использованием SIMD-расширений современных процессоров, легковесных GPU-потоков или аппаратных ускорителей, таких как Intel Xeon Phi. По аналогии выполняются и другие модульные операции — вычитание, деление без остатка и пр. Очевидно, что когда мы имеем дело с классической длинной арифметикой, добиться такой простоты и естественного параллелизма не представляется возможным (для сравнения можно рассмотреть алгоритмы
IntegerAddition и KaratsubaMultiply из [13]). Это обстоятельство делает весьма привлекательными методы многоразрядной арифметики на основе СОК. Однако они имеют свои недостатки, особенно в отношении таких вопросов, как выполнение немодульных операций (сравнение, определение знака, контроль переполнения и пр.). Эти операции требуют в той или иной форме оценки величины числа, представленного кортежем вычетов. Но такая оценка затруднена тем, что отдельные вычеты не позволяют судить о целочисленном значении, выражаемом ими. Пусть, например, в СОК с модулями 3, 5, 7 записаны числа X = (1,0, 6) и Y = (2,4, 0), которые нужно сравнить. Классический способ сделать это базируется на китайской теореме об остатках и состоит в вычислении двоичных представлений чисел с последующим их анализом. Так, после преобразования, X = 55, Y = 14, откуда X > Y. Такой способ, однако, является трудоемким, поскольку предполагает выполнение сложных операций умножения и сложения больших чисел и редукции по модулю M.
К числу других методов оценки величины в СОК относятся преобразование к смешанной системе (Mixed-Radix Conversion) [15,17,21], использование монотонных функций [22,23], ранг [16], контроль четности [24], ядро (Core Function) [18] и другие. Однако часто они оказываются непрактичными, так как требуют большого объема вычислений и / или хранения внушительных подстановочных таблиц.
Для решения обозначенной проблемы мы используем подход, который базируется на интервальной оценке относительных величин в СОК. В его основе лежит замечание, что для определения результата той или иной немодульной операции в большинстве случаев не требуется знать точные значения операндов, достаточно лишь иметь информацию о диапазонах их изменения. Причем, чтобы исключить необходимость работы с большими числами в диапазоне [0, M), где M — произведение всех модулей, имеет смысл оперировать масштабированными относительно M значениями [25], изменяющимися в диапазоне [0,1). В частности, возвращаясь к рассмотренному выше примеру, для сравнения чисел X и Y не обязательно преобразовывать их в двоичную систему. Достаточно знать, к примеру, что Х/М находится в интервале от 0.5 до 0.7, а Y/M — в интервале от 0.1 до 0.3. Тогда можно с уверенностью сказать, что X > Y, так как Х/М > Y/M. Подобную информацию дает ИПХ I(Х/М) = [Х/М, Х/М], которая для заданного X представляет собой интервальную функцию от остатков Х\,Х2,... ,хп, в общем случае удовлетворяющую соотношениям
Х/М < Х/М < Х/М. Границы ИПХ хранятся в виде чисел с плавающей точкой машинной точности и вычисляются с использованием направленных округлений: Х/М — с округлением вниз (toward -то), а Х/М — с округлением вверх (toward +то).
В дальнейшем будем обозначать символами , V, V, V операции сложения, вычитания, умножения и деления соответственно, выполняемые с округлением вниз. Аналогично, символы Д, Д, Д, Д будут использоваться для обозначения операций, выполняемых с округлением вверх. Кроме этого будем считать, что если перед некоторым выражением стоит символ V или Д, то все операции в рамках этого выражения выполняются с соответствующим округлением. Вопросы выполнения направленных округлений рассмотрены в работах [26,27].
Пусть в СОК с модулями т\, То2,..., топ дано X = {х\, Х2,..., хп). Тогда простейший метод вычислить границы его ИПХ, непосредственно следующий из китайской теоремы об остатках, определяется следующей парой соотношений, в которых | • |i — дробная часть:
Классическое последовательное вычисление выражений (2) и (3) выполняется за время О(п), при условии, что каждая операция производится в двоичной арифметике машинной точности за постоянное время. При распараллеливании эта оценка сокращается до О (log п). Кроме этого, в [28] предложен алгоритм, позволяющий вычислить ИПХ практически с любой априорно заданной точностью на всем диапазоне СОК1 без использования многоразрядных типов данных.
Если динамический диапазон СОК сравнительно небольшой (порядка 1000 бит), то для представления границ ИПХ достаточно стандартного формата IEEE двойной точности. Если же требуется существенно больший диапазон (М ^ 21000), то следует использовать представление с расширенной экспонентой (extended-range floatingpoint). Такое представление может быть достаточно просто построено
том числе и для малых чисел, для которых относительная ошибка ИПХ при вычислении с использованием выражений (2) и (3) значительно возрастает.
(2)
1
(3)
1
посредством объединения машинного целого г с обычным машинным числом с плавающей точкой f, и рассмотрения этой пары как числа
f х В\
где В — константа, обозначающая основание системы счисления [29]. Алгоритмы обработки таких представлений легки в реализации и не приводят к значительным вычислительным издержкам, поскольку не предполагают увеличения точности (разрядности мантиссы f ).
Важным достоинством ИПХ является возможность использования следующих интервальных соотношений:
(4) I(Х/М) +1(Y/M) = [ Х/М V Y/M, Х/М △ Y/M ] ,
(5) I(Х/М) - I(Y/M) = [ Х/М V Y/M, Х/М △ Y/M_] ,
(6) I(Х/М) х I(Y/M) = [Х/М V YTM V УМ, Х/М △ Y/M △ 1/М ],
(7) I(Х/М)/I(Y/M) = [Х/М V УЖ VV YM, Х/М △ УМ △ Y/Ж],
где I(1/М) = [ 1/М, 1/М] — интервальная аппроксимация константы 1/М. Предполагается, что Y/M, Y/M = 0 в формуле (7).
Благодаря свойствам направленных округлений [26, раздел 2.2], если Х/М G I(Х/М),Y/M G I(Y/M), 1/М G I(УМ), результатный интервал гарантированно будет содержать математически точное значение, даже если операции над границами неточны. Иными словами, при использовании выражений (4)—(7) гарантируется, что
(8) (X о Y)/М G I(Х/М) о I(Y/M)
где о G {+, —, х, /}. Свойство (8) известно также как свойство интервального включения (inclusion property) [27], с поправкой на относительность выражаемых величин. Помимо этого, выражения (4)—(7) примечательны тем, что не ограничивают возможные значения результатной ИПХ диапазоном [0,1). Это позволяет, в частности, легко и непосредственно отслеживать переполнение или определять знак результата арифметических операций, несмотря на цикличность модулярного диапазона. Другие функции для оценки величины в СОК не позволяют сделать этого (по крайней мере, в такой простой форме), поскольку они вычисляются непосредственно по значениям вычетов.
ИПХ играет важную роль в MF-формате, позволяя эффективно выполнить целый ряд немодульных процедур. Например, для того чтобы определить, превышает ли мантисса X некоторый предопределенный предел L, достаточно сравнить I(Х/М) с заранее вычисленным интервалом I(L/M): если Х/М > L/M, то X > L, и наоборот, если
ТАБЛИЦА 1. Допустимые кодировки в МР-формате
Класс Знак, s Порядок,е Знак и мантисса, ±Х
Числа {0, 1} 6min 6 6max ± (®1,Ж2, ...,хп)
{0, 1} 6max + 1 ±(0, 0,.. ., 0)
Not-a-Number {0, 1} 6max + 1 ± (Ж1,Ж2, ... ,хп), 3xi = 0
Х/М < L/M, то X < L. Если I(Х/М) и I(L/M) пересекаются, что однозначно определяется сравнением противоположных границ интервалов, то имеется возможность использовать другие методы оценки величины X. Более детальные правила применения ИПХ даны в [19].
3. Числовые кодировки
В MF-формате могут быть представлены конечные числа (включая знаковые нули), бесконечности со знаком и специальные нечисловые данные (Not-a-Number). Их кодировки приведены в таблице 1.
Как и в стандарте IEEE 754 [30], бесконечности интерпретируются в аффинном смысле, т.е. —то < х < +то, где х — любое конечное число. Основное предназначение бесконечностей состоит в том, что они позволяют получить хотя бы близкий к правильному результат вычисления в случае переполнения. Мантиссы нулей и бесконечностей представлены одинаковым кодом (0,0,..., 0), поэтому для однозначной интерпретации необходим анализ порядка е: нулю соответствует значение е = 0, а бесконечности — значение е = emax + 1. На уровне алгоритмов числовой обработки в MF-формате реализована стандартная арифметика бесконечностей. Возникновение NaNs2 является свидетельством некорректной операции или математической неопределенности результата, например, при попытке вычисления то х 0.
4. Диапазон, переполнение и потеря значимости
Множество чисел в MF-формате однозначно определяется следующими целочисленными параметрами:
• mi, Ш2,..., тп — попарно взаимно простые модули СОК;
• emin — минимальное значение порядка;
• emax — максимальное значение порядка.
2В отличие от стандарта IEEE 754 явного разделения на signaling NaNs и quiet NaNs здесь не предусмотрено.
Обозначим символами ш и П наименьшее и наибольшее по абсолютному значению конечные ненулевые числа. Поскольку мантисса X может принимать целочисленные значения из интервала [0, М — 1], где М — произведение всех TOj, то из (1) непосредственно следует, что ш =1 х 26min, а П = (М — 1) х 26max. Полный диапазон представимых чисел определяется интервалом [—0, П].
Ситуация, когда ненулевой результат х операции по абсолютному значению получился строго больше П, интерпретируется как арифметическое переполнение (overflow). В этом случае, в зависимости от режима округления, х заменяется бесконечностью или числом П. Ситуация, когда х по абсолютному значению строго меньше ш, интерпретируется как потеря значимости (underflow). В этом случае х заменяется нулем или числом ш, в зависимости от режима округления.
Следует отметить, что MF-формат ориентирован, прежде всего, на программную реализацию, вследствие чего предполагается, что для хранения порядка е используется отдельная переменная — машинное целое число со знаком. Это обеспечит очень большое расстояние между ш и П, достаточное для того, чтобы свести к минимуму вероятность переполнения или потери значимости. Например, если порядок представлен 32-битным числом и М « 21024, то ш « ю-646456993, П « 10646457300. Соотношение между этими предельными величинами составляет приблизительно 101292914293. Сложно себе представить обстоятельства, при которых результат вычислений, пусть даже промежуточный, выйдет за эти пределы. Для сравнения, отношение расстояния от Земли до наиболее удаленного из известных квазара ULAS J1120+0641 к радиусу протона составляет лишь 1.43 х 1041.
Вместе с тем, учитывая, что в общем случае для порядка е может быть отведено менее 32 бит, полностью исключить возможность переполнения или потери значимости было бы опрометчиво. Поэтому необходимы алгоритмы, позволяющие однозначно определить факт возникновения одной из этих ситуаций. Они должны учитывать, что рассматриваемое числовое представление, строго говоря, не является нормализованным (по крайней мере, в терминах IEEE 754), поскольку мантисса почти всегда будет больше 2 (базы экспоненты), а контроль единицы в старшем бите ее двоичного представления является излишне затратным. Это допускает существование для одного и того же конечного значения более чем одной кодировки3. Например, если
^Исключением является нулевое значение, для которого предусмотрен один уникальный способ записи.
СОК задана модулями 3, 5, 7, то число 0.25 может быть записано в одной из следующих форм: (1,1,1} х 2-2, (2, 2, 2} х 2-3, (1,4, 4} х 2-4, (2, 3,1} х 2-5 и т.д. Впрочем, указанная особенность не приводит к каким-либо затруднениям или ограничениям при реализации алгоритмов арифметических вычислений. Однако при выборе ет1п и етах необходимо иметь ввиду, что е может принимать значения в диапазоне от ет;п — Ш до етах + Ш, где Ш = |_^2(М — 1)]. При этом число, которому соответствует этот порядок, может лежать в пределах конечного диапазона, то есть по абсолютному значению быть не меньше ш (по крайней мере, при X = М — 1) и не превышать П (при X = 1). Как следствие, анализ на предмет переполнения или потери значимости должен выполняться на основании совместного анализа порядка е и мантиссы X. В частности, в отношении переполнения справедливы следующие положения, инвариантные к форме представления числа (кодировке) и выводимые из простых арифметических соображений:
• если е < етах, то переполнения нет при любом значении X;
• если е > етах + Ш, то переполнение имеет место при любом X;
• если е = етах + г, где г — целое число в интервале от 1 до Ш включительно, то факт переполнения имеет место только тогда, когда X > |_(М — 1)/2^.
Алгоритм проверки переполнения представлен ниже. Он возвращает 1, если переполнение возникло, и 0 — в противном случае. Алгоритм определения потери значимости CнECкUNDERFLOW задается очень схожим образом.
Алгоритм 3. Определение переполнения
procedure CheckOwerflow(x, e) if e < emax then
return 0 end if
if e > emax + W then
return 1 end if
A ^ MixedRadix(x)
> Переполнения нет
> Переполнение > Преобразование из СОК в MRS
% ^ ^ 6max
В ^ TH > Выбор MR-кода |_(М — 1)/2i\
return Compare^, В) > 0 ? 1 : 0 end procedure
В строках 8-11 мантисса X сравнивается с числом |_(М — 1)/2*J посредством системы со смешанными основаниями (Mixed-Radix System, MRS). Мантисса преобразуется из СОК в MRS алгоритмически, а MR-представление числа |_(М — 1)/2*J выбирается из подстановочной таблицы T. Далее выполняется обычное попарное сравнение MR-цифр, начиная со старшей. Конечно здесь можно было использовать ИПХ-метод, что во многих случаях позволило бы ускорить процедуру сравнения. Однако быстродействие этой процедуры не имеет принципиального значения, поскольку в подавляющем большинстве случаев алгоритм CheokOwerflow завершится раньше.
5. Основные алгоритмы
Некоторые из рассматриваемых далее алгоритмов ранее были представлены в работах [19,20]. Здесь они приведены с оптимизациями и усовершенствованиями, не меняющими сути.
Чтобы не перегружать псевдокод малозначимыми операциями копирования всюду далее предполагается, что операнды передаются по значению, а не по указателю, поэтому их локальные модификации не распространяются за пределы рассматриваемой подпрограммы.
5.1. Округление
Рассматриваемое числовое представление позволяет реализовать ряд эффективных решений при обработке многоразрядных чисел, не свойственных позиционной длинной арифметике. Одно из таких решений — новая схема округления, представленная на рис. 2.
Рис. 2. Схема округления
Согласно этой схеме решение о необходимости округления принимается на основании анализа мантисс операндов непосредственно перед арифметической операцией, а не после ее выполнения. Такая схема предпочтительна с точки зрения точности и снижения количества итераций масштабирования мантиссы, так как позволяет избежать ненужных округлений, если известно, что мантисса результата последующей операции будет находиться в допустимом диапазоне [0, М — 1]. Наибольший эффект от использования такой схемы достигается в аддитивных операциях, которые не приводят к существенному росту мантиссы результата. В качестве недостатка стоит отметить, что использование такой схемы в некоторой степени затрудняет априорный анализ ошибок округления. Эффективная реализация представленной схемы возможна благодаря наличию малоразрядной ИПХ в представлении числа, которая позволяет дать быструю оценку величины мантиссы непосредственно перед выполнением операции.
Алгоритм 4 выполняет округление. Он принимает на вход число х ^ [я,Х,е,1 (Х/М)} и количество разрядов округления г, которое определяется в вызывающей подпрограмме. Для ускорения масштабирования мантиссы степенью двойки (строка 7) разработан ускоренный итерационный метод [20], рассмотрение которого выходит за рамки данной работы. Отметим лишь, что он не требует преобразования из СОК в двоичную систему счисления и допускает эффективное распараллеливание вычислений.
Алгоритм 4. Округление
1 procedure Round(x,t)
2 if г < 0 then
3 return x
4 end if
5 sr ^ s
6 er ^ e + r
7 R ^ X/2r > С округлением до целого
8 if CheokOverflow(R, er) = 1 then > Переполнение
9 return TrapOwerflow(x)
10 end if
11 Вычисление I(R/M) > Алгоритм из [28]
12 return {sr, R, er,1 (R/M)}
13 end procedure
ТАБЛИЦА 2. Спецификация операции умножения модулей чисел |ж| и \у\
х +0 number +ТО NaN
+0 +0 +0 NaN NaN
number +0 М х 1у1 +ТО NaN
+ТО NaN +ТО +ТО NaN
NaN NaN NaN NaN NaN
Метод масштабирования из [20], вычисляет число \_Х/2ТПоэтому при его использовании представленный Алгоритм 4 работает в режиме «round toward zero», т.е. позволяет округлять до меньшего по модулю числа. Вместе с тем, реализация режима «round to nearest» (округление до ближайшего) не вызывает затруднений. Для этого достаточно модифицировать алгоритм следующим образом: после масштабирования мантиссы следует вычислить остаток в СОК Т = (t\,t2,... ,tn) = X — 2r х R и сравнить его с числом С = 2Г-1 (модулярный код для С может быть вычислен заранее и записан в подстановочную таблицу). Если окажется, что Т > С, то записать R ^ R + 1. Более того, при использовании ИПХ можно вообще не вычислять непосредственно остаток Т. Достаточно рассчитать I(Т/М) = I(Х/М) — 2Г х I(RJM), используя интервальное соотношение (5), и сравнить полученный интервал с заранее вычисленным I(СУМ) = \CSM_, С/М]. При этом если на аппаратном уровне поддерживается операция Fused Multiply-Add (FMA), то для вычисления I(Т/М) потребуется выполнить лишь две инструкции.
Если мантисса на шаге 7 округляется до ближайшего или большего по модулю числа, то возможно переполнение. В этом случае управление передается обработчику TrapOverflöw. При округлении до меньшего по модулю переполнения, очевидно, не возникнет, так как |_^/2rJ х 2e+r < X х 2е < П. Также не возникнет и потери значимости, вне зависимости от используемого режима округления.
5.2. Умножение
Таблица 2 содержит спецификацию операции |ж| х |у|. В ней «number» обозначает ненулевое конечное число. Эта таблица соответствует стандарту IEEE 754 [26, раздел 8.4]. Знак результата вычисляется сложением по модулю два знаков сомножителей.
Алгоритм 5 выполняет умножение. Он принимает сомножители х ^ isœ, X, ех, I(Х/М)}, у ^ {sy, Y, еу, I(Y/M)} и возвращает произведение z ^ {sz, Z, ez, I(Z/M)}. При умножении экстремально больших чисел может возникнуть переполнение. И напротив, если сомножители очень малы, то возможна потеря значимости. В любой из этих ситуаций управление передается соответствующему обработчику.
Благодаря использованию эффективной процедуры перемножения мантисс (RnsProd), умножение чисел многократной точности в MF-формате, как показывают тесты, выполняется в среднем на порядок быстрее, по сравнению с традиционными алгоритмами.
Алгоритм 5. Умножение
procedure Multiply(x, у)
Фильтрация вырожденных случаев I (Z/M) i I (Х/М) х I (Y/M) if Z/M > then
x i CheokAndRound(x) > Округление
у i CheokAndRound^) I (Z/M) i I (Х/М) х I (Y/M) end if
Sz i Sx ф Sy ez i cx + Cy
Z i RnsProd(x, Y) > Параллельное умножение
if CheokOverflow(z, ez) = 1 then > Переполнение
return TrapOwerflow(z) end if
if CheokUnderflow(z, ez) = 1 then > Потеря значимости
return TrapUnderflow(z) end if
return z ^ jsz, Z, ez, I (Z/M)} end procedure
В строке 2 производится анализ сомножителей на предмет соответствия одной из кодировок таблицы 1. Если хотя бы один из операндов равен нулю, бесконечности или МаМ, то, в соответствии с указанной таблицей, формируется вырожденный результат операции и алгоритм завершается. При этом может быть установлен статусный флаг, свидетельствующий о возникшем исключении.
Перед умножением производится проверка необходимости округления и округление операндов (строки 3-8). Поясним суть этой процедуры. Вначале определяется возможность выхода мантиссы произведения Z за диапазон представления посредством вычисления ИПХ в соответствии с формулой (6) и сравнения ее верхней границы с нижней границей заранее вычисленного интервала, локализующего точку (М — 1)/М. Если Z/M < M—1, то Z гарантированно будет лежать в пределах диапазона (действительно, выполнение этого неравенства означает, с учетом свойства интервального включения, что точное значение Z/M не превысит (М — 1)/М, а следовательно, Z не превысит М — 1). В противном случае Z может выйти за пределы диапазона. Для исключения такой возможности мантиссы сомножителей округляются до М' = [%/М — 1J. Количество разрядов округления г для каждого сомножителя определяется соотношением (на примере х)
(9) X/2r < М', откуда получается точное выражение для г:
(10) г = [log2 (Х/М')! = riog2 X — log2 м'1 .
Если г < 0, то округление не требуется. Вычислить (10) можно несколькими способами. Простейший способ сводится к преобразованию мантиссы X из СОК в двоичную систему счисления. Однако его высокая вычислительная сложность сводит на нет смысл его использования. Более эффективный способ получается, если обе части соотношения (9) поделить на М:
г = [log2(X/M) — log2(M '/М )1 .
Значение log2(X/M) аппроксимируется посредством логарифмирования ИПХ I(Х/М); при вычислении логарифмов ИПХ используются соответствующие направленные округления:
log2 Х/М < log2 (Х/М) < log2 Х/М.
Алгоритм 6 выполняет округление числа х ^ [s,X,e,I(Х/М)} перед умножением. Заранее вычисляется I (log 2 а) = [log2 a, log2 а\, где а = М'/М. Процедура Check возвращает пару чисел (г, d), где г — количество разрядов округления, а d — максимально возможное количество избыточных разрядов округления, показывающее на сколько в худшем случае вычисленное г может отличаться от точного значения, определяемого в соответствии с соотношением (10). Число d может быть использовано для контроля точности вычислений. Если
ИПХ мантиссы вычислена с достаточной точностью, то d = 0. Если же d превышает некоторый предел, то ИПХ требует уточнения с последующим повторным выполнением процедуры Check.
Алгоритм 6. Проверка и округление при умножении 1: procedure CheokAndRound(x)
2: (г, d) i Cheok(x) > г — число разрядов округления
3: Анализ d и оценка точности 4: return Round(x,t) 5: end procedure
6: procedure Cheok(x) 7: b i log2 X/M 8: С i log2 X/M 9: r i \c A log2 a] 10: ri i \b ▽ log2 a] 11: if r > 0 then 12: d i Г — ri
13: return (r, d)
14: else
15: return (0, 0)
16: end if
17: end procedure
5.3. Выравнивание порядков
Если порядки операндов различны, то перед сложением или вычитанием должно производиться их выравнивание. которое целесообразно осуществлять до меньшего из порядков, что подразумевает умножение мантиссы операнда с большим порядком на 2|Ле|, где Де — разность порядков. Подобный способ выравнивания используется в десятичной арифметике с плавающей точкой [31]. Однако если разность Де по абсолютному значению велика, то при умножении может произойти переполнение диапазона. Поэтому здесь используется компенсационная схема выравнивания, при которой происходит одновременное умножение мантиссы операнда с большим порядком и масштабирование мантиссы другого операнда коэффициентом, определяемым остаточной разностью порядков, которая не может быть учтена при умножении. Данная схема представлена на рис. 3.
х2'r
X -х2 X
-О-•-
div by 2r
Y <- Y
-•-O-
_ -(Ae - r
—ey — ex
О Before alignment # After alignment
Рис. 3. Схема выравнивания порядков при Де > 0
с
e
y
x
Пусть х и у - числа с мантиссами X и Y, порядками ех и еу, интервально-позиционными характеристиками I(Х/М) и I(Y/M) соответственно, и пусть Де = ех — еу > 0. Определим наименьшее целое г, такое, что при умножении X на 2Ле-г не произойдет переполнения. Тогда 2Г будет определять коэффициент масштабирования Y. Очевидно, что г находится по аналогии с числом разрядов округления в Алгоритме 6 (поэтому используется одно и то же обозначение). Так, должно выполняться соотношение
X M — 1 2 ~ 2Де '
где M - произведение модулей СОК. После преобразований, принимая во внимание, что г должно быть целым числом, получим
г > Г— log2(M — 1) + Де + log2 M + log2(X/M)].
Приняв допущение, что M достаточно большое для того, чтобы разницей между log2 (M — 1) и log2 M можно было пренебречь, получаем следующее точное выражение для г:
(11) г = [log2(X/M) + Де].
Наконец, умножив X на 2Ле-т и разделив Y на 2Г, получим операнды с одинаковым порядком ех — (Де — г) = еу + г. Вместо вычисления log2(X/M) будем вычислять log2 Х/М, где Х/М - верхняя граница интервально-позиционной характеристики I(Х/М).
Алгоритм 7 выполняет выравнивание порядков. Предполагается, что на вход поступают конечные числа, а фильтрация нулей и специальных величин выполняется в вызывающей программе.
Алгоритм 7. Выравнивание порядков
1: procedure Alignment(x, y)
2:
3: if Ae = 0 then
4: return x, y
5 end if
6: if Ae < 0 then
7: t ^ X, X ^ y, y ^ t > х и у меняются местами
8: Ae <--Ae
9: end if
10: if X/M < 2-A e then
11: r ^ 0 > Округление числа у не требуется
12: else
13: r ^ [log2 X/M A Ae]
14: c ^ log2 Y/M
15: if c < r A log2 M then > Сравниваем Y и 2Г
16: y ^ (-1)«« x 0
17: return x, y
18: end if
19: end if
20: X ^ Rnsprod(x, 2ae-r)
21: I (X/M) ^ I (X/M) x 2a e-r
22: ^ - (Ae - r)
23: y ^ Round(^, r) > Округление у на г бит
24: return x,y
25: end procedure
После выполнения алгоритма будут получены ненулевые числа с одинаковыми порядками, либо одно из них будет нулем (то, чей порядок был меньшим), а другое не изменится. При округлении у (если оно необходимо) переполнения не возникнет, так как увеличенный порядок еу, будет, по крайней мере, не больше исходного ех. Для ускорения умножения мантиссы (строка 20) нужно записать в подстановочную таблицу модулярные коды всех натуральных степеней двойки до 2^°82мJ. Размер таблицы равен |_1°ё2 М\ х п.
Так же как и в Алгоритме 6, при выравнивании достаточно просто может быть реализован контроль точности вычислений. Для этого нужно вычислить Г1 = |~log2 Х/М ▽ Де] и сравнить это значение с
Таблица 3. Выбор алгоритма сложения / вычитания чисел
Операция Знаки операндов Выполняемый алгоритм
ADD Совпадают Addition
ADD Различны Subtration
SUB Совпадают Subtration
SUB Различны Addition
ТАБЛИЦА 4. Спецификация операции сложения модулей чисел |ж| (по строкам) и |у| (по столбцам)
+ +0 number NaN
NaN +TO NaN
+0 +0 \y\ NaN
number M M + \у\ NaN
+TO +TO +TO +TO NaN
NaN NaN NaN NaN NaN
числом г, вычисляемым на шаге 13. Если разность г — ri оказывается неприемлемо большой, следует уточнить I(Х/М). Аналогичная проверка может быть реализована на шаге 14 при принятии решения об обнулении операнда с меньшим порядком.
5.4. Сложение
Знак числа представляется отдельно от мантиссы, поэтому достаточно реализовать алгоритмы сложения (Addition) и вычитания (Subtration) чисел одинаковых знаков. Выбор необходимого алгоритма осуществляется исходя из типа выполняемой операции (ADD или SUB) и знаков операндов в соответствии с таблицей 3 [32, с. 419].
Таблица 4 содержит спецификацию операции \ж\ + \у\. Она соответствует стандарту IEEE 754 [26, раздел 8.3]. В этой таблице «number» обозначает ненулевое конечное число. Знак результата сложения определяется общим знаком слагаемых.
Алгоритм 8 выполняет сложение чисел одинаковых знаков. На вход подаются х ^ {sx,X,ex,I(Х/М)} и у ^ jsy,Y,ey,1 (Y/M)}, в
результате возвращается сумма z ^ [sz,Z,ez,1 (Z/M)}. Алгоритм вычитания Sübtration очень схож с алгоритмом сложения. При этом знак разности определяется на основании ИПХ, вычисляемой по формуле (5). Заметим, что при вычитании переполнения не происходит, однако возможна потеря значимости.
Алгоритм 8. Сложение
1
2
3
4
5
6
7
8 9
10 11 12
13
14
15
16
17
18
19
20
procedure Addition(x, у)
Фильтрация вырожденных случаев
Alignment(x, у)
if х = 0 then return у
else if у = 0 then return x
end if
I(Z/M) i I(X/M) + I(Y/M) if Z/M > M]-1 then x i Round(x, 1) у i Round(^, 1) I (Z/M) i I (X/M) + I (Y/M) end if sz i sx Cz i Cx
Z i RnsSum(x, Y) if CheokOyerflow(z, ez) = 1 then
return TrapOwerflow(z) end if
return z ^ jsz, Z, ez, I (Z/M)} end procedure
> Округление на 1 разряд
> Sx Sy
> Cx
> Параллельное сложение
6. Оценка эффективности
Рассмотренные алгоритмы реализованы в высокоточной программной библиотеке языка С МЕ-ЫЬгагу [19]. Для их оценки были выполнены следующие эксперименты:
(1) Исследование точности вычислительных операций.
(2) Вычисление полинома Румпа.
(3) Вычисление рекуррентного соотношения Мюллера.
(4) Исследование эффективности 81МБ-распараллеливания.
(5) Исследование производительности операций СЕМУ и СЕММ.
Задачи 2 и 3 являются классическими примерами расчетов, которые даже при незначительном объеме вычислений сопровождаются возникновением катастрофических погрешностей. Такие задачи часто используются для верификации высокоточных программных средств.
Во всех экспериментах система остаточных классов для представления мантисс была задана 32 15-битными модулями, обеспечивающими диапазон М « 2479. Вычисления выполнялись на машине следующей конфигурации: Intel Core i7-4702MQ / 4 Cores / 6 MB Cache / 4 GB RAM / Intel C++ Compiler XE 14.0.
6.1. Точность вычислительных операций
Мантиссы чисел в MF-формате могут принимать значения в диапазоне [0, М — 1]. В общем случае для предотвращения переполнения производится их округление до [л/М — 1J. Таким образом, можно дать грубую априорную оценку точности р вычислений (в битах):
(12) р = [log2[VM-1JJ.
Например, если М = 2479, то р = 239. Однако используемая схема округления (см. пункт 5.1) позволяет производить округление лишь при необходимости, для предотвращения переполнения, что обеспечивает в ряде случаев получение результата удвоенной точности. Например, если анализ интервально-позиционных характеристик показал, что мантиссы исходных операндов не превышают [л/М — 1J, то при их умножении округления не будет, следовательно, будет получен точный результат (относительная ошибка равна нулю, если исходные операнды являются точными). Вероятность получения результата удвоенной точности возрастает при выполнении аддитивных операций (сложение, вычитание), когда порядки операндов примерно равны.
Вычислявшиеся в ходе эксперимента выражения и результаты эксперимента представлены в таблице 5. Исходными данными (xi,yi,Zi) для выражений были псевдослучайные 239-битные числа, равномерно распределенные в диапазоне (0,1), полученные с помощью алгоритма «Вихрь Мерсенна». Число слагаемых в каждой вычисляемой сумме и сомножителей в каждом вычисляемом произведении равнялось 106. Каждое выражение вычислялось 100 раз. Для оценки результатов использовалась библиотека MPFR [33], реализующая корректное с точки зрения стандарта IEEE 754 округление. Использовался режим округления «round toward zero». Оценивалась максимальная апостериорная относительная ошибка: е = max |(ж — х')/ж|, где х — точный результат, х' — вычисленный результат.
Таблица 5. Результаты экспериментальной оценки точности арифметики с представлением чисел в МТ-формате
№ Выражение
Относительная ошибка е
MF-Library (239 бит) MPFR (239 бит)
1 Y,xi
2 £ 1/xi
3 хгУг
4 ^гУг
5 Y,(1/Xi — 1/Vi)2
6 £ 1/(г + Xi)2
7 Y\xi
8 П Хг + У г
9 U(Xi — Уг )2
2
0
2.203 х 10-74 2.803 х 10-139 0
1.400 х 10-72 1.827 х 10-72 6.483 х 10-67 6.688 х 10-67 1.487 х 10-66
3.866 х 10-67 4.342 х 10-67 3.866 х 10-67 1.133 х 10-72 5.919 х 10-67 1.022 х 10-66 8.172 х 10-67 1.036 х 10-66 1.633 х 10-66
Таким образом, в общем случае полученные результаты согласуются с априорной оценкой (12). Более точные результаты в задачах 1, 3 и 4 получены благодаря отмеченной выше особенности округления.
6.2. Полином Румпа
Вычислялся следующий полином (S.M. Rump, 1988 [34]):
(13) f (а, Ь) = 333.75b6 + а2(11а2Ь2 — Ь6 — 12164 — 2) + 5.5Ъ8 + а/2Ъ,
при а = 77617.0 и Ь = 33096.0. В качестве эталонного использовалось решение, полученное с использованием 4096-битной арифметики МРЕИ. Результаты представлены в таблице 6. Здесь, как и в предыдущем эксперименте, более высокая точность МЕ-ЫЬгагу по сравнению с МРЕИ (256 бит) объясняется используемой схемой округления.
6.3. Рекуррентное соотношение Мюллера
Другой пример расчетов, сопровождающихся потерей точности, который приводит профессор КаЬап в работе [35] — это вычисление следующего, на первый взгляд тривиального, соотношения:
х0 = 4.00, Х1 = 4.25, Xi = f (Xi-1,Xi-2).
К. С. Исупов, В. С. Князьков Таблица 6. Результаты вычисления полинома (13)
Тип данных (точность) Результат Верных цифр
Float (24 бит) -6.33825 х 1029 0
Double (53 бит) -1.18059 х 1021 0
MPFR (64 бит) -5.76461 х 1017 0
MPFR (128 бит) 0.8273960599... 39
MPFR (256 бит) 0.8273960599... 77
MF-Library (239 бит) 0.8273960599... 140
MPFR (512 бит) 0.8273960599... 155
Параметры соотношения (14) подобраны так, что при точных вычислениях lim хп = 5.0. Однако из-за влияния ошибок округления
(а конкретно, из-за невозможности точного представления = 76/17 в виде числа с плавающей точкой) последовательность |жп} отходит от верного ответа к точке х* = 100.0. Номер итерации, с которой начинается такой переход, прямо пропорционален точности вычислений. Результаты эксперимента представлены на рис. 4. Графики показывают, что точность MF-Library (239 бит) сравнима с точностью MPFR (256 бит).
6.4. Эффективность SIMD-распараллеливания
SIMD-параллелизм — важная составляющая современных вычислений. В разной степени большинство из современных процессорных архитектур поддерживают SIMD-расширения или SIMD-инструкции. Производители процессоров возлагают большие надежды на этот принцип вычислений, систематически наращивая разрядность векторных регистров. В свою очередь, появляются все более совершенные команды SIMD-обработки, которые позволяют одновременно обрабатывать все большее количество пар операндов, а современные компиляторы позволяют производить автоматическую векторизацию при условии подходящих для этого структур данных. Эффективность алгоритмов с точки зрения SIMD на центральных процессорах может косвенно свидетельствовать о целесообразности реализации этих алгоритмов на других аппаратных платформах, таких как графические ускорители. Эти обстоятельства нельзя не принимать во внимание при разработке программного обеспечения и, в частности, программных библиотек
120 ■ 100 ■ 806040200-20-
100 ■ 90807060. 50403020100-
Float (24 bit)
1
0 10 20 30 40 50 60 70 п
MPFR (128 bit)
180 160 140 120 100 = 80 60 40 20 0
Double (53 bit)
0 10 20 30 40 50 60 70
100-, MPFR (256 bit)
10 20 30 40 50 60 70
10 20 30 40 50 60 70
5.04.94.8: 4.74.64.54.4-
MPFR (512 bit)
100-| MF-Llbrary (239 bit) 90-
10 20 30 40 50 60 70
10 20 30 40 50 60 70
Рис. 4. Результаты вычисления соотношения (14)
многократной точности. Однако известные библиотеки, ориентированные на центральные процессоры, как правило, не задействуют 81МБ-параллелизм, либо задействуют его лишь в незначительной степени, что обусловлено вполне очевидными причинами, связанными, прежде всего, с образованием множественных зависимостей по данным при обработке мантисс большой разрядности. При вычислениях в МЕ-формате отдельные цифры мантисс, представленных в СОК, взаимно независимы, что позволяет организовать их эффективную векторную обработку.
1.82
ЩЗ.60
0.13 I 0.15
□
12.08
□
6.16
□
13.14
7.95
□
16.28
—i—i—i—i—|—|—г-
10 15
MFLOPS
Serial
^SIMD
25.22
20
25
Рис. 5. Производительность высокоточных арифметических операций с представлением чисел в МР-формате
Для оценки эффективности SIMD-распараллеливания запуск MF-Library выполнялся в двух конфигурациях: при установленном запрете на векторизацию (прописыванием директив #pragma novector) и при использовании средств автоматической векторизации компилятора Intel C+—К В последнем случае векторизовались циклы вычисления модулярных мантисс и ИПХ. Исследованы операции сложения (add), вычитания (sub), умножения (mult), деления (div) и умножения с накоплением (mac, х ^ х + yz). Результаты представлены на рис. 5.
При векторизации всех исследованных операций, кроме деления, получено двукратное ускорение, что является довольно высоким показателем для данной тестовой конфигурации. При векторизации деления ускорение отсутствует. Это связано с особенностями используемого алгоритма (см. [20]), приводящими к необходимости затратного преобразования модулярных мантисс в двоичную систему счисления с применением многоразрядной позиционной арифметики. Лучшие показатели могут быть достигнуты при использовании более эффективного алгоритма деления. Можно ожидать большего ускорения, если вместо автоматической векторизации, выполняемой средствами компилятора, вручную использовать SIMD-расширения команд, а также при выполнении на системе, содержащей более широкие векторные регистры и более совершенные SIMD-команды.
О 100 200 300 400 500 600 700 800 900 1000
П
Рис. 6. Время последовательного и параллельного (4 ядра) выполнения операции GEMM с точностью 239 бит
п
Рис. 7. Время последовательного и параллельного (4 ядра) выполнения операции ОБЫУ с точностью 239 бит
6.5. Быстродействие BLAS
Быстродействие MF-Library исследовано на операцях обобщенного матричного умножения (GEMM, BLAS Level-3) и общего векторного умножения (GEMV, BLAS Level-2) с точностью 239 бит. Рассматривались классические последовательные, а также параллельные алгоритмы, что позволило дополнительно оценить потоковую безопасность. При выполнении GEMM (С ^ аАВ + ) матрицы А, В, С были плотными, их порядок п изменялся с шагом 50 в интервале от 100 до 1000. Для верификации результата вычислялась норма ||С||i. Данные эксперимента представлены на рис. 6. Среднее ускорение MF-Library по сравнению с MPFR v. 3.0.1 на операции GEMM составило 1.9 раза и 2.3 раза соответственно при последовательных и параллельных вычислениях. При распараллеливании скорость MF-Library возросла в 5.3 раза (сверхлинейное ускорение объясняется поддержкой Hyper-Threading). При выполнении GEMV (у ^ аАх + fy) векторы х,у и матрица А также были плотно заполнены. Размер векторов варьировался в диапазоне от 500 до 1500 с шагом 100. Корректность операции оценивалась по норме ||y||i. Результаты экспериментов представлены на рис. 7. В данном случае по сравнению с MPFR получено ускорение 2.4 раза и 2.6 раза соответственно при последовательных и параллельных вычислениях. При распараллеливании скорость библиотеки MF-Library увеличилась в среднем в 4 раза.
Другие результаты экспериментального исследования эффективности рассматриваемого в данной работе подхода в сравнении с аналогами представлены в работе [19].
7. Направления дальнейших исследований
Рассмотренная в настоящей работе техника вычислений многократной точности появилась сравнительно недавно, что обуславливает весьма широкий перечень требующих проработки вопросов, как исключительно алгоритмических, так и по части разработки программного обеспечения. К числу наиболее приоритетных мы относим следующие.
(1) Развитие методов вычисления и анализа интервально-позиционных характеристик для быстрого выполнения немодульных операций в СОК. Эта задача включает в себя: • определение отличительных, инвариантных к виду модулей и их количеству, признаков чисел в СОК, лежащих вблизи границ диапазона представления, и разработка на их основе универсального способа классификации чисел по вычетам;
• реализацию алгоритмов обработки чисел с расширенным порядком (extended-range floating-point), необходимых для корректной работы с интервально-позиционными характеристиками в больших динамических диапазонах СОК (Р > 21000);
• разработку параллельных алгоритмов уточненного вычисления интервально-позиционной характеристики;
• разработку и обоснование новых ускоренных методов модулярного масштабирования, не требующих хранения больших подстановочных таблиц и не накладывающих специфических требований на масштабирующий коэффициент.
(2) Получение строгих оценок для верхних границ ошибок округления, возникающих при выполнении основных арифметических операций. Эта задача играет важную роль в вопросах воспроизводимости и достоверности результатов вычислений.
(3) Оптимизация алгоритмов с целью более эффективного использования SIMD-параллелизма.
(4) Адаптация форматов, структур данных и алгоритмов для эффективной реализации на графических процессорах и ускорителях типа Intel Xeon Phi. Разработка соответствующих подпрограмм.
(5) Разработка эффективного алгоритма деления. Текущая реализация приводит к необходимости преобразования мантисс из СОК в позиционную систему счисления, вследствие чего является крайне медленной. Имеет смысл, вероятно, отказаться от прямых методов деления в пользу быстрых численных методов, таких как метод Ньютона-Рафсона или метод Гольдшмидта [36]. Их основное преимущество, помимо достаточно быстрой сходимости, состоит в том, что они не требуют трудоемкого деления мантисс в СОК, а сводятся к выполнению последовательности сложений и умножений, которые эффективно распараллеливаются.
(6) Расширение функциональности посредством разработки и программной реализации методов высокоточного вычисления алгебраических и трансцендентных (логарифмических и показательных, тригонометрических, гиперболических и соответствующих им обратных) функций, наиболее употребительных плотных и разреженных BLAS-операций, а также ряда численных методов высокой точности, в частности, для решения систем уравнений.
(7) Комплексная отладка и оптимизация кода, тестирование.
8. Заключение
Камнем преткновения на пути создания эффективной высокоточной арифметики является обработка многоразрядных мантисс, алгоритмически сложная и приводящая к большим временным затратам. Для решения этой проблемы мы использовали системы остаточных классов (СОК), которые позволяют представить мантиссу произвольной длины в виде набора взаимно-независимых малоразрядных вычетов с возможностью их эффективной параллельной обработки.
Традиционно СОК применялись при решении целочисленных задач в таких областях, как цифровая обработка сигналов, криптография, обнаружение и исправление ошибок в кодах, обработка изображений. Использование СОК в арифметике с плавающей точкой потребовало решения ряда специфических проблем, не свойственных вычислениям с целыми числами, таких как округление, выравнивание порядков, контроль над переполнением числового диапазона. Для этих целей в числовое представление введена атрибутивная информация — интервально-позиционная характеристика, позволяющая быстро оценить величину мантиссы по без полного ее преобразования в двоичную систему. В результате получен рассмотренный в этой работе модулярно-позиционный формат с плавающей точкой — MF-формат.
Алгоритмы высокоточной арифметики в MF-формате сочетаются с особенностями современных параллельных вычислительных архитектур, позволяя эффективно использовать SIMD, многопоточный, и (при достаточно большом числе модулей СОК) многоядерный параллелизм. Кроме этого, они поддерживают арифметику бесконечностей и обработку основных исключительных ситуаций, определенных стандартом IEEE 754 применительно к двоичной арифметике с плавающей точкой. Также в этой работе показана возможность контроля за точностью вычислений, имеющая важное значение при выполнении априорного анализа ошибок округления вычислительных алгоритмов.
Определены приоритетные задачи дальнейших исследований. Ожидается, что результатом решения этих задач будет новое быстродействующее кроссплатформенное программное обеспечение параллельной арифметики многократной точности для гибридных вычислительных систем, эффективно задействующее доступные ресурсы параллелизма как центральных, так и графических процессоров. Кроме этого, в перспективе мы планируем рассмотреть возможности аппаратной реализации разработанных алгоритмов для FPGA.
Список литературы
[1] S. Collange, D. Defour, S. Graillat, R. Iakymchuk. "Reproducible and Accurate Matrix Multiplication for High-Performance Computing", 16th GAMM-IMACS International Symposium on Scientific Computing, Computer Arithmetic and Validated Numerics, SCAN'2014 (Wurzburg, Germany, 2014), pp. 42-43. t 61
[2] S. Collange, D. Defour, S. Graillat, R. Iakymchuk. "Numerical Reproducibility for the Parallel Reduction on Multi- and Many-Core Architectures", Parallel Computing, 49 (2015), pp. 83-97. t 61
[3] J. Demmel, H. D. Nguyen. "Parallel Reproducible Summation", IEEE Transactions on Computers, 64:7 (2015), pp. 2060-2070. t 61
[4] A. Pescaru, E. Oanta, T. Axinte, A. D. Dascalescu. "Extended Precision Data Types for the Development of the Original Computer Aided Engineering Applications", IOP Conference Series: Materials Science and Engineering, 95:1 (2015), URL: http://stacks.iop.org/1757-899X/95/ i=1/a=012125t 62
[5] D. H. Bailey, R. Barrio, J. M. Borwein. "High-Precision Computation: Mathematical Physics and Dynamics", Applied Mathematics and Computation, 218:20 (2012), pp. 10106-10121. t 62
[6] D. H. Bailey, J. M. Borwein. "High-Precision Arithmetic in Mathematical Physics", Mathematics, 3 (2015), pp. 337-367. t 62
[7] D. Viswanath, S. Sahutoglu. "Complex Singularities and the Lorenz Attractor", SIAM Review, 52 (2010), pp. 294-314. t 62
[8] S. Leweke, E. von Lieres. "Fast Arbitrary Order Moments and Arbitrary Precision Solution of the General Rate Model of Column Liquid Chromatography With Linear Isotherm", Computers & Chemical Engineering, 84 (2016), pp. 350-362. t 62
[9] Y. He, C. Ding. "Using Accurate Arithmetics to Improve Numerical Reproducibility and Stability in Parallel Applications", The Journal of Supercomputing, 18 (2001), pp. 259-277. t 62
[10] C. F. Berger, Z. Bern, L.J. Dixon, C. F. Febres, D. Forde, H. Ita, D. A. Kosower, D. Maitre. "Automated Implementation of On-shell Methods for One-loop Amplitudes", Physical Review D, 78 (2008). t 62
[11] J.-P. Merlet. "On the Real-Time Calculation of the Forward Kinematics of Suspended Cable-Driven Parallel Robots", 14th World Congress in Mechanism and Machine Science (Taipei, Taiwan, 2015). t 62
[12] В. Л. Якушев, В. Н. Симбиркин, А. В. Филимонов, П. А. Новиков, И. Н. Коньшин, Г. Б. Сушко, С. А. Харченко. «Решение плохообу-словленных симметричных СЛАУ для задач строительной механики параллельными итерационными методами», Вестник ННГУ, 2012, №4(1), с. 238-246. t 62
[13 [14
[15 [16 [17 [18 [19
[20
[21 [22
[23
[24
[25
[26
R. Brent, P. Zimmermann. Modern Computer Arithmetic, Cambridge University Press, New York, NY, USA, 2010, 236 p. t 62,66 Н. Н. Непейвода, И. Н. Григоревский, Е. П. Лилитко. «О представлении действительных чисел», Программные системы: теория и приложения, 5:4(22) (2014), с. 105-121, URL: http://psta.psiras.ru/read/psta2014_ 4_105-121.pdf t 62
N. S. Szabo, R. I. Tanaka. Residue Arithmetic and its Application to Computer Technology, McGraw-Hill, New York, NY, USA, 1967, 236 p.t63'66 И. Я. Акушский, Д. И. Юдицкий. Машинная арифметика в остаточных классах, Сов. Радио, М., 1968, 440 с. t 63,66
B. Parhami. Computer Arithmetic: Algorithms and Hardware Designs, Oxford University Press, Oxford, 2000, 490 p. t 63,64,66
A. Omondi, B. Premkumar. Residue Number Systems: Theory and Implementation, Imperial College Press, London, UK, 2007, 312 p. t 63,66 K. Isupov, V. Knyazkov. "A Modular-Positional Computation Technique for Multiple-Precision Floating-Point Arithmetic", Parallel Computing Technologies, Lecture Notes in Computer Science, vol. 9251, Springer International Publishing, 2015, pp. 47-61. t 64,69,72,81,88 К. С. Исупов, А. Н. Мальцев. «Способ представления чисел с плавающей точкой большой разрядности, ориентированный на параллельную обработку», Вычислительные методы и программирование, 15:4 (2014), с. 631-643. t 64,72,73,74,86
M. Akkal, P. Siy. "A New Mixed Radix Conversion Algorithm MRC-II", Journal of Systems Architecture, 53:9 (2007), pp. 577-586. t 66 G. Dimauro, S. Impedovo, G. Pirlo. "A New Technique for Fast Number Comparison in the Residue Number System", IEEE Transactions on Computers, 42:5 (1993), pp. 608-612. t 66
G. Pirlo, S. Impedovo. "A New Class of Monotone Functions of the Residue Number System", International Journal of Mathematical Models and Methods in Applied Sciences, 7:9 (2013), pp. 802-809. t 66 J.-H. Yang, C.-C. Chang, C.-Y. Chen. "A High-Speed Division Algorithm in Residue Number System Using Parity-Checking Technique", International Journal of Computer Mathematics, 81:6 (2004), pp. 775-780. t 66
C. Y. Hung, B. Parhami. "An Approximate Sign Detection Method for Residue Numbers and Its Application to RNS Division", Computers & Mathematics with Applications, 27:4 (1994), pp. 23-35. t 66
J.-M. Muller, N. Brisebarre, F. de Dinechin, C.-P. Jeannerod, V. Lefevre, G. Melquiond, N. Revol, D. Stehle, S. Torres. Handbook of Floating-Point Arithmetic, Birkhauser, Boston, 2010, 572 p. t 67,68,74,80 U. Kulisch. Computer Arithmetic and Validity — Theory, Implementation and Applications, de Gruyter, Berlin, 2008, 409 p. t 67,68
[28] К. С.Исупов. «Об одном алгоритме сравнения чисел в системе остаточных классов», Вестник АГТУ. Серия «Управление, вычислительная техника, и информатика», 2014, №3, с. 40-49. t 67,73
[29] J. R. Hauser. "Handling Floating-Point Exceptions in Numeric Programs", ACM Transactions on Programming Languages and Systems, 18:2 (1996), pp. 139-174. t 68
[30] IEEE Standard for Floating-Point Arithmetic. Introduced 2008-08-29, Institute of Electrical and Electronics Engineers, New York, NY, USA, 2008, 70 p. t 69
[31] M. F. Cowlishaw. "Decimal Floating-Point: Algorism for Computers", 16th IEEE Symposium on Computer Arithmetic, ARITH-16'03, IEEE Computer Society, Washington, DC, USA, 2003, pp. 140-147. t 77
[32] M. D. Ercegovac, T. Lang, Digital Arithmetic, The Morgan Kaufmann Series in Computer Architecture and Design, Morgan Kaufmann, San Francisco, 2004, 709 p. t 80
[33] L. Fousse, G. Hanrot, V. Lefevre, P. Pelissier, P. Zimmermann. "MPFR: A Multiple-Precision Binary Floating-Point Library With Correct Rounding", ACM Transactions on Mathematical Software, 33:2 (2007). t 82
[34] S. M. Rump. "Algorithms for Verified Inclusions — Theory and Practice", Reliability in Computing: The Role of Interval Methods in Scientific Computing, ed. R. E. Moore, Academic Press Professional, San Diego, CA, USA, 1988, pp. 109-126. t 83
[35] W. Kahan. How Futile are Mindless Assessments of Roundoff in Floating-Point Computation, https://www.cs.berkeley.edu/~wkahan/ Mindless.pdf, 2006. t 83
[36] I. Kong, E. E. Swartzlander. "A Goldschmidt Division Method With Faster Than Quadratic Convergenc", IEEE Transactions on Very Large Scale Integration (VLSI) Systems, 19:4 (2011), pp. 696-700. t 89
Рекомендовал к публикации Программный комитет
Третьего национального суперкомпьютерного форума НСКФ-2015
Об авторах:
Константин Сергеевич Исупов
Ведущий научный сотрудник кафедры электронных вычислительных машин Вятского государственного университета, кандидат технических наук. Область научных интересов: программное обеспечение высокопроизводительных систем, параллельные алгоритмы, высокоскоростная компьютерная арифметика, вычисления многократной точности.
e-mail: [email protected]
Владимир Сергеевич Князьков
Директор научно-образовательного центра супервычислительных технологий и систем, профессор кафедры электронных вычислительных машин Вятского государственного университета, доктор технических наук. Область научных интересов: параллельные вычислительные архитектуры, высокопроизводительные вычисления, реконфигурируемые вычислительные системы.
e-mail: [email protected]
Пример ссылки на эту публикацию:
К. С. Исупов, В. С. Князьков. «Арифметика многократной точности на основе систем остаточных классов», Программные системы: теория и приложения, 2016, 7:1(28), с. 61-97.
URL: http://psta.psiras.ru/read/psta2016_1_61-97.pdf
Konstantin Isupov, Vladimir Knyazkov. Parallel multiple-precision arithmetic based on residue number system,.
Abstract. This paper deals with algorithms of multiple-precision arithmetic, based on the use of multi module residue number systems for representing of arbitrary length significands of floating-point numbers; the exponent is represented in the binary number system. Such number representation provides a large dynamic range and allows for effective parallelization of arithmetic operations on the digits of multiple-precision significands across RNS modules. This agrees well with the architectural features of modern parallel computing systems. Additionally, the attributive information which provides a fast estimation for the relative value of significand and allows you to increase the speed of executing complex non-modular operations in RNS, such as comparison, overflow control, rounding, etc., is included into the number format. Results of an experimental study on precision, performance and SIMD efficiency of multiple-precision algorithms are presented. (In Russian).
Key words and phrases: computer arithmetic, high-precision computations, parallel algorithms, residue number system, SIMD.
References
[1] S. Collange, D. Defour, S. Graillat, R. Iakymchuk. "Reproducible and Accurate Matrix Multiplication for High-Performance Computing", 16th GAMM-IMACS International Symposium on Scientific Computing, Computer Arithmetic and Validated Numerics, SCAN'2014 (Wurzburg, Germany, 2014), pp. 42-43.
[2] S. Collange, D. Defour, S. Graillat, R. Iakymchuk. "Numerical Reproducibility for the Parallel Reduction on Multi- and Many-Core Architectures", Parallel Computing, 49 (2015), pp. 83-97.
[3] J. Demmel, H. D. Nguyen. "Parallel Reproducible Summation", IEEE Transactions on Computers, 64:7 (2015), pp. 2060-2070.
[4] A. Pescaru, E. Oanta, T. Axinte, A. D. Dascalescu. "Extended Precision Data Types for the Development of the Original Computer Aided Engineering Applications", IOP Conference Series: Materials Science and Engineering, 95:1 (2015), URL: http://stacks.iop.org/1757-899X/95/i=1/a=012125
[5] D. H. Bailey, R. Barrio, J. M. Borwein. "High-Precision Computation: Mathematical Physics and Dynamics", Applied Mathematics and Computation, 218:20 (2012), pp. 10106-10121.
[6] D. H. Bailey, J. M. Borwein. "High-Precision Arithmetic in Mathematical Physics", Mathematics, 3 (2015), pp. 337-367.
[7] D. Viswanath, S. Sahutoglu. "Complex Singularities and the Lorenz Attractor", SIAM Review, 52 (2010), pp. 294-314.
[8] S. Leweke, E. von Lieres. "Fast Arbitrary Order Moments and Arbitrary Precision Solution of the General Rate Model of Column Liquid Chromatography With Linear Isotherm", Computers & Chemical Engineering, 84 ( 2016), pp. 350-362.
[9] Y. He, C. Ding. "Using Accurate Arithmetics to Improve Numerical Reproducibility and Stability in Parallel Applications", The Journal of Supercomputing, 18 (2001), pp. 259-277.
© K. Isupov, V. Knyazkov, 2016
© Vyatka State University, 2016
© Program systems: Theory and Applications, 2016
[10] C. F. Berger, Z. Bern, L. J. Dixon, C. F. Febres, D. Forde, H. Ita, D. A. Kosower, D. Maitre. "Automated Implementation of On-shell Methods for One-loop Amplitudes", Physical Review D, 78 (2008).
[11] J.-P. Merlet. "On the Real-Time Calculation of the Forward Kinematics of Suspended Cable-Driven Parallel Robots", 14th World Congress in Mechanism and Machine Science (Taipei, Taiwan, 2015).
[12] V. L. Yakushev, V. N. Simbirkin, A. V. Filimonov, P. A. Novikov, I. N. Kon'shin, G. B. Sushko, S. A. Kharchenko. "Solution of ill-conditioned symmetric slae for structural mechanics problems by parallel iterative methods", Vestnik NNGU, 2012, no.4(1), pp. 238-246 (in Russian).
[13] R. Brent, P. Zimmermann. Modern Computer Arithmetic, Cambridge University Press, New York, NY, USA, 2010, 236 p.
[14] N. N. Nepeyvoda, I. N. Grigorevskiy, Ye. P. Lilitko. "New representation of real numbers", Programmnye Sistemy: Teoriya i Prilozheniya, 5:4 (2014), pp. 105-121 (in Russian), URL: http://psta.psiras.ru/read/psta2014_4_105-121.pdf
[15] N. S. Szabo, R. I. Tanaka. Residue Arithmetic and its Application to Computer Technology, McGraw-Hill, New York, NY, USA, 1967, 236 p.
[16] I. Ya. Akushskiy, D. I. Yuditskiy. Machine Arithmetic in Residue Classes, Sov. Radio, M., 1968 (in Russian), 440 p.
[17] B. Parhami. Computer Arithmetic: Algorithms and Hardware Designs, Oxford University Press, Oxford, 2000, 490 p.
[18] A. Omondi, B. Premkumar. Residue Number Systems: Theory and Implementation, Imperial College Press, London, UK, 2007, 312 p.
[19] K. Isupov, V. Knyazkov. "A Modular-Positional Computation Technique for Multiple-Precision Floating-Point Arithmetic", Parallel Computing Technologies, Lecture Notes in Computer Science, vol. 9251, Springer International Publishing, 2015, pp. 47-61.
[20] K. S. Isupov, A. N. Mal'tsev. "A parallel-processing-oriented method for the representation of multi-digit floating-point numbers", Vychislitel'nyye metody i programmirovaniye, 15:4 (2014), pp. 631-643 (in Russian).
[21] M. Akkal, P. Siy. "A New Mixed Radix Conversion Algorithm MRC-II", Journal of Systems Architecture, 53:9 (2007), pp. 577-586.
[22] G. Dimauro, S. Impedovo, G. Pirlo. "A New Technique for Fast Number Comparison in the Residue Number System", IEEE Transactions on Computers, 42:5 (1993), pp. 608-612.
[23] G. Pirlo, S. Impedovo. "A New Class of Monotone Functions of the Residue Number System", International, Journal of Mathematical Models and Methods in Applied Sciences, 7:9 (2013), pp. 802-809.
[24] J.-H. Yang, C.-C. Chang, C.-Y. Chen. "A High-Speed Division Algorithm in Residue Number System Using Parity-Checking Technique", International Journal of Computer Mathematics, 81:6 (2004), pp. 775-780.
[25] C. Y. Hung, B. Parhami. "An Approximate Sign Detection Method for Residue Numbers and Its Application to RNS Division", Computers & Mathematics with Applications, 27:4 (1994), pp. 23-35.
[26] J.-M. Muller, N. Brisebarre, F. de Dinechin, C.-P. Jeannerod, V. Lefèvre, G. Melquiond, N. Revol, D. Stehle, S. Torres. Handbook of Floating-Point Arithmetic, Birkhauser, Boston, 2010, 572 p.
[27] U. Kulisch. Computer Arithmetic and Validity — Theory, Implementation and Applications, de Gruyter, Berlin, 2008, 409 p.
[28] K. S. Isupov. "On an algorithm for number comparison in the residue number system", Vestnik AGTU. Seriya "Upravleniye, vychislitel'naya tekhnika i informatika", 2014, no.3, pp. 40—49 (in Russian).
[29] J. R. Hauser. "Handling Floating-Point Exceptions in Numeric Programs", ACM Transactions on Programming Languages and Systems, 18:2 (1996), pp. 139—174.
[30] IEEE Standard for Floating-Point Arithmetic. Introduced 2008-08-29, Institute of Electrical and Electronics Engineers, New York, NY, USA, 2008, 70 p.
[31] M. F. Cowlishaw. "Decimal Floating-Point: Algorism for Computers", 16th IEEE Symposium on Computer Arithmetic, ARITH-16'03, IEEE Computer Society, Washington, DC, USA, 2003, pp. 140-147.
[32] M. D. Ercegovac, T. Lang, Digital Arithmetic, The Morgan Kaufmann Series in Computer Architecture and Design, Morgan Kaufmann, San Francisco, 2004, 709 p.
[33] L. Fousse, G. Hanrot, V. Lefevre, P. Pelissier, P. Zimmermann. "MPFR: A Multiple-Precision Binary Floating-Point Library With Correct Rounding", ACM Transactions on Mathematical Software, 33:2 (2007).
[34] S. M. Rump. "Algorithms for Verified Inclusions — Theory and Practice", Reliability in Computing: The Role of Interval Methods in Scientific Computing, ed. R. E. Moore, Academic Press Professional, San Diego, CA, USA, 1988, pp. 109-126.
[35] W. Kahan. How Futile are Mindless Assessments of Roundoff in Floating-Point Computation, https://www. cs.berkeley.edu/~wkahan/Mindless.pdf, 2006.
[36] I. Kong, E. E. Swartzlander. "A Goldschmidt Division Method With Faster Than Quadratic Convergenc", IEEE Transactions on Very Large Scale Integration (VLSI) Systems, 19:4 (2011), pp. 696-700.
Sample citation of this publication:
Konstantin Isupov, Vladimir Knyazkov. "Parallel multiple-precision arithmetic based on residue number system", Program systems: theory and applications, 2016, 7:1(28), pp. 61-97. (In Russian).
URL: http://psta.psiras.ru/read/psta2016_1_61-97.pdf
