The Multiplication Method with Scaling the Result for High-Precision Residue Positional Interval Logarithmic Computations

Full Text

Введение

Рост вычислительных мощностей современных компьютеров делает возможным решение прикладных задач сверхбольшой размерности с огромным количеством вычислительных операций. Неконтролируемые ошибки округления, методологически присущие стандарту вещественных чисел IEEE 754, не позволяют решить проблему точности, достоверности и воспроизводимости вычислений при решении задач данного класса [1–5]. В частности, для решения в современных постановках задач в области экспериментальной вычислительной математики [1; 2], математической физики [4], биохимии [3], астрофизики и получения достоверных результатов требуются операции с числами длиной от 100 до 1 000 десятичных цифр (с использованием специально разработанных программных библиотек высокоточных вычислений). В связи с этим актуальными направлениями исследований являются теория и способы практической реализации вычислительной математики многократной точности (высокоточная, или длинная, арифметика), оперирующей с числами произвольной длины в сверхбольших числовых диапазонах.

Для решения задач в сверхбольших числовых диапазонах в настоящее время применяются такие специализированные программные пакеты высокоточных вычислений, как ARPREC, MPFUN90, DDFUN, FMLIB, FMZM90, QD, GMP, MPFR++, NTL, PARI/GP, CLN, HPAlib, Predicates, GARPREC, GQD, MatLab, Matematica, Maple [6]. Перечисленные программные решения базируются на специально разработанных многоразрядных форматах (128-, 256-, 512-битная (и более) арифметика) в базисах классических позиционных систем счисления и правил вычислений стандарта IEEE 754. Эти решения, благодаря наличию высокоуровневых программных интерфейсов и широкого спектра реализованных библиотек математических функций, являются наиболее популярными.

Недостатком современных пакетов высокоточных вычислений является резкое снижение скорости вычислений при обработке длинных многоразрядных операндов. При выходе длины операндов за диапазон представления данных в стандарте IEEE 754 скорость вычислений снижается в десятки и тысячи раз [2; 7] из-за необходимости алгоритмической обработки цепочек межзнаковых переносов длинных операндов. В итоге время решения задачи становится неприемлемым для практической деятельности.

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

Можно выделить два направления работ, направленных на повышение скорости вычислений при выполнении расчетов в сверхбольших числовых диапазонах. Первое направление ориентировано на модернизацию и создание новых технологий гибридных вычислений и обработки данных: численная и нечисленная обработка данных реализуется в гибридных системах кодирования с использованием «гибридных» наборов операций (сигнатур) и правил их выполнения. Математические методы и их алгоритмические решения для гибридных технологий вычислений ориентируются на программную реализацию на вычислительных платформах универсального назначения и, как правило, опираются на правила выполнения операций стандарта IEEE 754. Примером успешного использования этого подхода является библиотека высокоточной модулярно-позиционной арифметики [8], где использованы системы счисления в остаточных классах (СОК), вычисления в интервальной арифметике и позиционная система счисления стандарта IEEE 754.

Вторым направлением работ для повышения скорости вычислений в сверхвысоких числовых диапазонах является разработка специализированных средств аппаратной поддержки операций над сверхбольшими операндами, разрядность которых многократно превышает разрядную сетку современных индустриальных процессоров. Популярной технологической базой для создания таких спецпроцессоров «длинной» арифметики являются программируемые логические интегральные схемы (FPGA) и системы на кристалле [9–14]. Применение таких спецпроцессоров позволяет сократить время расчетов по сравнению с программными решениями в несколько десятков раз, но недостатки, присущие позиционной длинной арифметике, сохраняются [13]. Эти методологические недостатки позиционной арифметики приводят к необходимости построения на аппаратном уровне исполнительных устройств высокой сложности, что в конечном итоге делает невозможным создание применимых технических решений. Данная проблема частично решается введением специализированных вычислительных конвейеров; однако, как показано в работе китайских ученых [12], подобный подход также ведет к резкому увеличению аппаратных затрат, поэтому на практике число ступеней сокращается до четырех сегментов.

В связи с этим при создании средств аппаратной поддержки длинной арифметики актуален подход, ориентированный на создание вычислительных платформ, поддерживающих на аппаратном уровне технологии гибридных вычислений, что позволяет сократить аппаратные затраты по сравнению с позиционной системой счисления. Серьезный вклад в развитие этого направления внесли И. Я. Акушский, Д. И. Юдицкий, В. М. Амербаев и целый ряд не менее значимых специалистов. Наиболее широко в системах гибридных вычислений используются системы счисления в остаточных классах [8; 15; 16] и логарифмические системы счисления [17–19]. Например, системы остаточных классов успешно используются для решения задач криптографии [20; 21] и цифровой обработки сигналов [22–24].

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

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

Обзор литературы

Основной проблемой высокоточных расчетов в сверхбольших числовых диапазонах с применением вычислительных операций по правилам стандарта IEEE 754 является выполнение контроля ошибок округления, контроля переполнения диапазона и масштабирования чисел при выполнении аддитивных и мультипликативных операций. Особенно это касается длительных итерационных и автоматных вычислений с накоплением при обработке массивов данных большого объема. Накопление ошибок при некорректно организованном контроле приводит к получению недостоверных результатов. Для обеспечения требуемой точности, достоверности и воспроизводимости расчетов в настоящее время применяются вычисления с использованием длинной позиционной арифметики, реализованной в современных специализированных библиотеках высокоточных вычислений. Основным недостатком современных библиотек длинной арифметики является неприемлемое для практики увеличение времени решения прикладных задач [4; 6]. Так, в работе китайских и итальянских ученых [7] для задач оптимального управления время вычислений возрастает с 5 с при использовании стандартного типа данных двойной точности до 980 с при использовании точности в 128 бит и до 35 ч – при использовании точности в 400 бит. Аналогичные результаты представлены в работе А. Вороса [2], где время вычисления возрастает с 4 минут до 22,5 дней. В задачах криптографии проблема ускорения арифметических операций над длинными целыми числами является не менее острой, чем в задачах моделирования [13].

Задача повышения скорости вычислений при выполнении расчетов в сверхбольших числовых диапазонах частично решается за счет применения специализированных процессоров-ускорителей для поддержки вычислений с использованием длинной арифметики. Например, японские ученые [9] представили семейство процессоров на базе FPGA, реализующих длинную арифметику типа «double-double» и «quad-double». Скорость решения задачи вычисления интегралов Фейнмана [10] с использованием данных процессоров приблизительно в 80–200 раз выше, чем скорость расчета с применением программных реализаций таких вычислений. Американскими учеными [11] приведены результаты реализации на FPGA целочисленных вычислений на длинных (разрядность – 64 000 бит) операндах в сравнении с вычислениями с применением библиотеки GMP: расчеты ускорились минимум в 5 раз при операциях сложения/вычитания и в 9 раз – при операции умножения.

Для ускорения выполнения операции умножения длинных чисел (1 024–2 048 бит) китайскими учеными [12] представлен конвейерный метод на базе 64-разрядных умножителей с глубиной конвейера до четырех ступеней (увеличение числа ступеней приводит к неоправданному росту аппаратных затрат). Как показывает анализ исследований, в области аппаратных решений умножителей в основном применяются базовые алгоритмы умножения квадратичной сложности в позиционной системе счисления [13; 14], поскольку аппаратная реализация асимптотически быстрых алгоритмов затруднена [26].

Для ускорения выполнения арифметических операций (кроме операции деления) над длинными целыми числами наиболее эффективными с точки зрения аппаратных затрат являются модулярные системы счисления. Например, исследователями [23] представлено устройство эллиптической криптографии, ускоряющее выполнение операции умножения Монтгомери с использованием 40 модулярных 15-битовых каналов. В работе австралийских ученых [16] модулярный вычислитель имеет 108 модулярных каналов с разрядной сеткой в 19 бит каждый, что позволяет работать в диапазоне чисел до 2 048 бит. Другими авторами [20] представлены модулярные устройства, позволяющие работать в диапазоне 1 024–4 096 бит.

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

Алгоритмы масштабирования в СОК представлены в достаточно большом количестве исследований. Разработанные методы масштабирования либо предназначены для специальных наборов модулей [27; 28], либо для масштабирования используются специальные подстановочные таблицы [22; 29; 30]. Последний подход практически неприемлем для масштабирования модулярных чисел при использовании произвольных наборов модулей большой разрядности из-за огромного (до Тбайта) объема подстановочных таблиц.

Основной сложностью при выполнении масштабирования является операция расширения базиса. Методы расширения базиса были исследованы авторами статьи ранее [25]. В результате исследований было установлено, что наиболее быстрые методы расширения базиса выполняются с использованием приближенной китайской теоремы об остатках – вычисления так называемой позиционной характеристики модулярного числа. Учеными [8] представлена интервальная позиционная характеристика (ИПХ) числа, в которой использованы преимущества интервальной арифметики [31; 32]. Использование ИПХ позволяет учитывать в явном виде ошибки округления, а также определять достоверность вычисления данной величины. Главным недостатком ИПХ является необходимость использования операций с плавающей точкой с направленным округлением, в то время как все остальные операции в СОК выполняются над целыми числами малой разрядности.

Использование логарифмической системы счисления (ЛСС) позволяет упростить выполнение операций умножения и деления, включая масштабирование [17; 33]. ЛСС превосходят по скорости и энергоэффективности арифметику с плавающей точкой на низкой разрядности операндов: до 16 бит – на любом наборе арифметических операций [18], до 32 бит – с преобладанием операций умножения и деления [19]. Дальнейшее увеличение разрядности ЛСС приводит к экспоненциальному росту сложности выполнения операций сложения и вычитания, поэтому при больших разрядностях ЛСС значительно уступает позиционной арифметике и используется только в приложениях, не требующих высокой точности [19].

В данной статье предлагается объединить преимущества СОК, ЛСС и интервальных вычислений: для высокоточных вычислений в сверхбольших числовых диапазонах рекомендуется модулярно-позиционная интервально-логарифмическая форма представления чисел и апробация эффективности таких гибридных вычислений на примере выполнения операции умножения с масштабированием.

Материалы и методы

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

1. Мантисса вещественного числа представляется в виде целого числа в системе остаточных классов набором n остатков ‹m₁, m₂, ..., m_n› от деления позиционного значения мантиссы на каждый из n модулей {p₁, p₂, ..., p_n}

$\text{M}\stackrel{СОК}{\to }⟨m_1,m_2,\dots ,m_n⟩,$

где ${\text{m}}_i=\text{M}\text{ mod }{p}_i\equiv {\left|\text{M}\right|}_{p_i}$  – i-й остаток от деления числа M по i-му модулю p_i:

$\begin{array}{c}{\text{m}}_i={\left|\text{M}\right|}_{p_i}=\text{M}-\left[\frac{\text{M}}{p_i}\right]\cdot p_i,\\ i=\mathrm{1,}\mathrm{2,}\dots \mathrm{.,}n,\end{array}$

где $\left[\frac{M}{p_i}\right]$   – целая часть частного $\frac{M}{p_i}$  ; {p₁, p₂, ..., p_n} – набор оснований или базис СОК. При этом диапазон представления модулярных мантисс определяется произведением всех модулей СОК, то есть  $M\in \left[0;P=p_1\cdot p_2\cdot \dots \cdot p_n\right)$ . Для кодирования цифр мантиссы используются целые числа без знака, представленные в позиционной системе счисления, но операции над цифрами мантиссы выполняются по правилам модулярной арифметики. Любая модульная операция $\circ \in \left\{+\text{,}-\text{,}\times \right\}$  над двумя числами ‹x₁, x₂, ..., x_n› и ‹y₁, y₂, ..., y_n›, представленными в СОК, выполняется независимо по каждому модулю:

$\left\{z_1,z_2,\dots ,z_n\right\}=\left\{{\left|x_1\circ y_1\right|}_{p_1},{\left|x_2\circ y_2\right|}_{p_2},\dots ,{\left|x_n\circ y_n\right|}_{p_n}\right\}.$

2. Характеристика абсолютной величины мантиссы вещественного числа представляется в виде логарифмического интервала (в интервально-логарифмической системе счисления):

$M\stackrel{ИЛСС}{\to }\left\{\underset{¯}{L_M}=\underset{¯}{{\log }_{b}M};\overline{L_M}=\overline{{\log }_{b}M}\right\},$

где $\underset{¯}{{\log }_{b}M}$  , $\overline{{\log }_{b}M}$  – логарифм числа по основанию b, вычисленный с округлением к –∞ и +∞ соответственно; M – модуль числа, представленный в позиционной системе счисления. Для кодирования характеристики мантиссы вещественного числа используется двоичная позиционная система счисления, но операции над значениями характеристик чисел выполняются по правилам интервальной арифметики и логарифметики. Результат умножения двух логарифмических интервалов $\left[\underset{¯}{L_X}=\underset{¯}{{\log }_{b}X};\overline{L_X}=\overline{{\log }_{b}X}\right]$  и  $\left[\underset{¯}{L_Y}=\underset{¯}{{\log }_{b}Y};\overline{L_Y}=\overline{{\log }_{b}Y}\right]$ определяется следующим образом:

$\underset{¯}{L_Z}=\underset{¯}{{\log }_{b}X\cdot Y}=\underset{¯}{{\log }_{b}X+{\log }_{b}Y}=\underset{¯}{L_X}+\underset{¯}{L_Y},$

$\overline{L_Z}=\overline{{\log }_{b}X\cdot Y}=\overline{{\log }_{b}X+{\log }_{b}Y}=\overline{L_X}+\overline{L_Y.}$

3. Масштаб (порядок) числа представляется в позиционной системе счисления в виде целого числа со знаком; операции выполняются также в позиционной системе счисления.
4. Знак числа представляется в позиционной системе счисления в виде одноразрядного числа со знаком; причем знак равен –1, если число отрицательное, 1 – если число положительное, и 0 – в случае равенства числа нулю. Дополнительный признак нуля вводится с целью представления интервальной логарифмической характеристики нулевого операнда, для которого невозможно вычисление логарифма.

Таким образом, число в гибридной модулярно-позиционной интервально-логарифмической форме представляется в следующем виде:

$\text{X}\stackrel{\text{МПИЛ-СС}}{\to }\left[{\text{m}}_1,{\text{m}}_2,\dots ,{\text{m}}_n,,,\lambda ,\sigma \right],$  где M = m₁, m₂, ..., m_n – модулярная мантисса числа; λ – масштаб (порядок) числа; $\underset{\_}{L},\overline{L}$  – границы интервальной логарифмической характеристики мантиссы числа; σ – знак числа.

При этом позиционное значение мантиссы вещественного числа X определяется как [X · b^e], где e – целое число, определяемое необходимой точностью. Например, мантисса вещественного числа, представленного в формате с плавающей точкой стандарта IEEE 754 как $X={\left(-1\right)}^s\times 1.f\times 2^{E-E_0}$ , где s – знак числа; 1.f – нормализованная мантисса; E – E₀ – порядок. При переводе в гибридную форму вычисляется так:

м$M=2^t\times 1.f$ , где t – разрядность мантиссы, определяемая конкретным типом данных.

Итак, позиционное значение данного числа определяется следующим образом:

$X=\sigma \cdot b^{\lambda }\cdot {\displaystyle \sum }_{i=1}^n{\left|m_i\cdot {\left|P_i^{-1}\right|}_{p_i}\right|}_{p_i}\cdot P_i,$  где  $P_i=\frac{P}{p_i}$ , ${\left|P_i^{-1}\right|}_{p_i}$  – мультипликативная инверсия P_i по модулю p_i, определяемая из соотношения ${\left|P_i^{-1}\cdot P_i\right|}_{p_i}\equiv 1$  ;  $i\in \left[1\right.,\left.n\right]$ ; n – количество модулей.

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

Масштабирование модулярных чисел выполняется на основании общего алгоритма масштабирования: пусть K – коэффициент масштабирования; Y – результат масштабирования числа X коэффициентом K; тогда результат масштабирования вычисляется по формуле:

$Y=\frac{X-{\left|X\right|}_K}{K},$  где |X|_K – остаток от деления числа X по модулю K.

Для случая масштабирования модулярных чисел коэффициентом, взаимно простым с основаниями СОК, используется итерационный алгоритм на основе алгоритма, предложенного сингапурскими и австралийскими учеными [27; 29].

1. Определение |X|_K, или так называемый этап расширения базиса, – получение остатка x_n+1от деления числа, представленного в СОК остатками x₁, x₂, ..., x_n по модулям p₁, p₂, ..., p_n, на число p_n+1 = K.
2. Непосредственно масштабирование по каждому модулю выполняется по формуле:

$y_i={\left|{\left|x_i-{\left|X\right|}_K\right|}_{p_i}\cdot {\left|{\text{K}}^{-1}\right|}_{p_i}\right|}_{p_i},$  где |K^–1|_pi – мультипликативная инверсия по модулю p_i коэффициента K.

Основные алгоритмы расширения базиса, анализ их вычислительной сложности были рассмотрены авторами в прошлой работе [25]. В данной статье используется быстрый метод масштабирования на основании китайской теоремы об остатках (КТО).

Согласно КТО, позиционное значение числа $X\in \left[0\right.,\left.P\right)$ , представленного в СОК остатками ‹x₁, x₂, ..., x_n› по основаниям {p₁, p₂, ..., p_n}, вычисляется по формуле:

$X={\left|{\displaystyle \sum }_{i=1}^n{\left|x_i\cdot {\left|P_i^{-1}\right|}_{p_i}\right|}_{p_i}{P}_i\right|}_P={\displaystyle \sum }_{i=1}^n{\left|x_i\cdot {\left|P_i^{-1}\right|}_{p_i}\right|}_{p_i}\cdot P_i-R\cdot P,$

где $P_i=\frac{P}{p_i}$ , |P_i^–1 |_pi – мультипликативная инверсия P_i по модулю p_i; $i\in \left[1\right.,\left.n\right]$  ; n – количество модулей; R – позиционный индекс.

Зная значение коэффициента R, можно вычислить остаток от деления по новому основанию без перевода модулярного числа в позиционное представление:

${\left|X\right|}_{p_{n+1}}={\left|{\left|{\displaystyle \sum }_{i=1}^n{\left|x_i\cdot {\left|P_i^{-1}\right|}_{p_i}\right|}_{p_i}\cdot {\left|P_i\right|}_{p_{n+1}}\right|}_{p_{n+1}}-{\left|R\cdot {\left|P\right|}_{p_{n+1}}\right|}_{p_{n+1}}\right|}_{p_{n+1}}.$

Для вычисления коэффициента R авторами разработан алгоритм с использованием целочисленных интервалов на основе приближенной интервальной оценки величины:

$\tilde{X}={\displaystyle \sum }_{i=1}^n{\left|x_i\cdot P_i^{-1}\right|}_{p_i}\cdot \frac{1}{p_i}=\frac{X+R\cdot P}{P}=R+\frac{X}{P},$

где целую часть величины $\tilde{X}$  определяет коэффициент R, а дробную – значение X/P. Процесс вычисления коэффициента R с использованием вещественных интервалов с направленным округлением и необходимые условия корректности вычислений представлены в работе К. Исупова и В. Князькова [8]; метод вычисления коэффициента R с использованием целочисленных интервалов описан в патенте [34].

Результаты исследования

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

Для вычисления произведения $Z=\left[⟨z_1,z_2,\dots ,z_n⟩,\underset{¯}{L_Z},\overline{L_Z},{\lambda }_Z,{\sigma }_Z\right]$  чисел $X=\left[⟨x_1,x_2,\dots ,x_n⟩,\underset{¯}{L_X},\overline{L_X},{\lambda }_X,{\sigma }_X\right]$  и  $Y=\left[⟨y_1,y_2,\dots ,y_n⟩,\underset{¯}{L_Y},\overline{L_Y},{\lambda }_Y,{\sigma }_Y\right]$  необходимо:

– вычислить знак произведения σ_Z = = σ_X · σ_Y путем алгебраического умножения знаков сомножителей;

– вычислить верхнюю границу интервальной логарифмической характеристики результата  $\underset{¯}{L_Z}=\underset{¯}{L_X}+\underset{¯}{L_Y}$ путем алгебраического сложения значений нижних границ $\underset{¯}{L_X}$  и $\underset{¯}{L_Y}$  ИЛХ операндов в позиционной системе счисления;

– вычислить верхнюю границу ИЛХ результата $\overline{L_Z}=\overline{L_X}+\overline{L_Y}$  путем алгебраического сложения значений нижних границ $\overline{L_X}$  и $\overline{L_Y}$  ИЛХ операндов в позиционной системе счисления;

– вычислить порядок результата λ_Z = λ_X + λ_Y путем алгебраического сложения порядков сомножителей;

– выполнить умножение модулярных мантисс путем нахождения значений $z_i={\left|x_i\cdot y_i\right|}_{p_i}=x_i\cdot y_i-\frac{x_i\cdot y_i}{p_i}\cdot p_i$  для всех  $i\in \left[1;n\right]$ ; при этом вычисления выполняются над операндами, представленными в позиционной системе счисления по правилам модулярной арифметики.

В данной статье отсутствует описание обработки исключительных ситуаций, таких как получение машинного нуля, переполнение и т. п. Более подробно метод описан в патенте [34].

Следует отметить, что поскольку мантиссы чисел, представленные в СОК, ограничены диапазоном [0; P), то при выполнении умножения двух мантисс результат может выйти за пределы диапазона представления, то есть $M_Z=M_X\cdot M_Y\ge P$  . Для того чтобы мантисса результата была представима в СОК, необходимо выполнить операцию масштабирования:

$\frac{M_X\cdot M_Y}{b^a},$  где  $a=\left[{\log }_b\frac{M_X\cdot M_Y}{P}\right]$ ; M_X, M_Y – позиционные значения мантисс; P – позиционное значение диапазона представления; [ ] означает округление к наибольшему целому.

Рассмотрим предельный случай, когда числа из диапазона [0; P) могут появиться с равной вероятностью. Вероятность того, что произведение двух мантисс выйдет за пределы диапазона представления модулярных мантисс, то есть M_X · M_Y ≥ P, равна

$p\left(M_Z\ge P\right)\approx \frac{{\left(P-1\right)}^2-\left(P-1\right)\ln \left(P-1\right)}{P^2}\approx \frac{P-\ln P}{P}\approx 1.$

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

В случае использования ИЛХ коэффициент a рассчитывается следующим образом:

$a=\overline{L_X}+\overline{L_Y}-L_P,$   где $\overline{L_X}$  , $\overline{L_Y}$  – верхние границы интервальных логарифмических характеристик чисел X и Y; $L_P={\log }_{b}P$  – константа для конкретного диапазона представления.

Таким образом, при использовании ИЛХ для вычисления коэффициента масштабирования не требуется преобразования в позиционную систему счисления и вычисления логарифма.

При умножении модулярных мантисс целесообразно выполнять масштабирование обоих операндов до непосредственного выполнения умножения; причем, если величина обоих операндов превышает значение $\sqrt{P}$ , следует распределять коэффициент масштабирования между операндами таким образом, чтобы отмасштабированные операнды не превышали величину $\sqrt{P}$  :

$b^a=b^{a_X}+b^{a_Y},$   где b^aX – масштабирующий коэффициент, применяемый к первому сомножителю; b^aY – масштабирующий коэффициент, применяемый ко второму сомножителю; a_X и a_Y – значения, определяемые соотношениями ИЛХ операндов следующим образом.

Пусть $L_1=\overline{L_X}+\overline{L_Y}-L_P$  , $L_2=\overline{L_X}-\overline{L_Y}$  ; тогда $a_X=\frac{L_1+L_2}{2},a_Y=\frac{L_1-L_2}{2}$  . Если только один из операндов превышает величину $\sqrt{P}$  , к нему необходимо применить масштабирующий коэффициент b^a.

Таким образом, если $\overline{L_Z}\ge L_P$  , необходимо выполнить масштабирование мантисс операндов, после чего выполнить умножение отмасштабированных мантисс, а также скорректировать значение порядка результата λ_z = λ_z + L₁ и значение верхней и нижней границы интервальной логарифмической характеристики результата $\underset{¯}{L_Z}=\underset{¯}{L_Z}-L_1$  , $\overline{L_Z}=\overline{L_Z}-L_1$ .

Процесс масштабирования является итерационным, поскольку за один шаг выполняется масштабирование коэффициентом, не превышающим 2^q, где q – разрядность модулей СОК.

На каждом шаге масштабирования вычисляется значение коэффициента R и остаток от деления модулярного числа на p_n+1 = 2^α, где α ≤ q:

$x_{n+1}={\left|M_X\right|}_{2^{\alpha }}={\left|{\left|{\displaystyle \sum }_{i=1}^n{\left|{\left|x_i\cdot {\left|P_i^{-1}\right|}_{p_i}\right|}_{p_i}\right|}_{2^{\alpha }}\cdot {\left|P_i\right|}_{2^{\alpha }}\right|}_{2^{\alpha }}-{\left|R\cdot {\left|P\right|}_{2^{\alpha }}\right|}_{2^{\alpha }}\right|}_{2^{\alpha }}.$

Затем выполняется масштабирование коэффициентом 2^α:

${\tilde{x}}_i={\left|{\left|x_i-{\left|M_X\right|}_{2^{\alpha }}\right|}_{p_i}\cdot {\left|{\left(2^{\alpha }\right)}^{-1}\right|}_{p_i}\right|}_{p_i}$

где  ${\left|{(2^{\alpha })}^{-1}\right|}_{p_i}$  – мультипликативная инверсия числа 2^α по модулю p_i – константа для конкретного значения модуля p_i.

Все вычисления производятся над целыми числами, представленными в позиционной системе счисления, по правилам модулярной арифметики.

Если a_X > q, процедура масштабирования повторяется над уже масштабированной мантиссой ${\tilde{x}}_1,{\tilde{x}}_2,\dots ,{\tilde{x}}_n$ и так далее, пока не будет выполнено полное масштабирование коэффициентом 2^αX. Аналогичным образом выполняется масштабирование второго сомножителя Y.

Подробный алгоритм деления модулярной мантиссы числа, представленного в модулярно-позиционной интервально-логарифмической форме, на число 2^α (масштабирование степенью двойки), также представлен в патенте [34].

Среднее время выполнения разработанного метода равно $T=p\left(M_Z\ge P\right)\cdot t_1+\frac{n}{k}{t}_2,$  где p(M_Z ≥ P) – вероятность того, что произведение двух мантисс выйдет за пределы диапазона представления модулярных мантисс (в предельном случае p(M_Z ≥ P) = 1); t₁ – время выполнения операции масштабирования; t₂ – время выполнения операции умножения по модулю; k – количество параллельных модулярных каналов.

Среднее время выполнения операции масштабирования определяется следующим образом:

$t_1=j\cdot \left(t_3+\frac{n}{k}{t}_4+\frac{n}{k}{t}_5\right),$ 

где t₃ – время выполнения операции расширения базиса; t₄ – время выполнения операции вычитания по модулю; t₅ – время выполнения операции умножения по модулю на константу; j – число итераций масштабирования.

Диапазон представления модулярных мантисс  $P={\displaystyle \prod }_{i=1}^n{p}_i\approx 2^{n\cdot q}$ . Минимальный коэффициент масштабирования равен 2¹, максимальный равен  $2^{\frac{n\cdot q}{2}}$ . Максимальное количество шагов масштабирования примем равным  $\frac{n}{2}$ .

Таким образом, минимальное и максимальное время выполнения операции масштабирования равны соответственно:

$t_1{}_{min}=t_3+\frac{n}{k}{t}_4+\frac{n}{k}{t}_5,$

$t_1{}_{max}=\frac{n}{2}\left(t_3+\frac{n}{k}{t}_4+\frac{n}{k}{t}_5\right).$

Время выполнения операции расширения базиса равно [34]:

$t_3=\frac{n}{k}{t}_5+\frac{n}{k}{t}_6+t_7,$

где t₅ – время выполнения операции умножения по модулю на константу; t₆ – время выполнения операции скалярного произведения вектора на вектор-константу; t₇ – время выполнения операции сложения.

Время выполнения операции умножения по модулю двух произвольных чисел приблизительно в 2 раза выше, чем время выполнения стандартного целочисленного умножения; время выполнения операции умножения по модулю на константу приблизительно равно времени выполнения стандартного целочисленного умножения [13]. Таким образом, минимальное и максимальное время выполнения разработанного метода равно:

$T_{min}=\left(\frac{n}{k}{t}_5+\frac{n}{k}{t}_6+t_7+\frac{n}{k}{t}_4+\frac{n}{k}{t}_5\right)+\frac{n}{k}{t}_2=\left(\frac{2n}{k}{t}_c+\frac{n}{k}{t}_c+t_c+\frac{n}{k}{t}_{\text{c}}\right)+\frac{2n}{k}{t}_c=\left(\frac{6n}{k}+1\right)t_c,$

$T_{max}=\frac{n}{2}\left(\frac{n}{k}{t}_5+\frac{n}{k}{t}_6+t_7+\frac{n}{k}{t}_4+\frac{n}{k}{t}_5\right)+\frac{n}{k}{t}_2=\frac{n}{2}\left(\frac{2n}{k}{t}_c+\frac{n}{k}{t}_c+t_c+\frac{n}{k}{t}_{\text{c}}\right)+\frac{2n}{k}{t}_c=\left(\frac{5n^2+4n}{2k}+1\right)t_c,$

где t₇ – длительность такта.

Сравним разработанный метод с конвейерным методом умножения длинных чисел с плавающей точкой (обозначим его как Lei et al [12]), а также со стандартными методами умножения (обозначим их как Schulte ‒ Swarzlander [14] и Ishii [13]). В качестве контроля рассмотрим асимптотически быстрый метод, используемый для организации некоторых целочисленных двоичных умножителей (алгоритм Карацубы [26]). В табл. 1 представлены оценки времени выполнения разработанного метода и аналогов; n – количество слов базовой длины, необходимых для представления длинного числа; k – количество параллельных модулярных каналов. В табл. 2 выполнено сравнение времени выполнения разработанного метода и аналогов для чисел разрядности 1 024 и 2 048 бит (n = 16 64-разрядных слов и n = 32 64-разрядных слова соответственно). На рис. 1; 2 представлены расчеты времени выполнения разработанного метода и аналогов для разрядности сомножителей от 256 до 2 048 бит с использованием 64-разрядных (рис. 1) и 16-разрядных слов (рис. 2) соответственно.

 

 

 

Рис. 1. Сравнение быстродействия разработанного метода с аналогами
(с использованием 64-разрядных умножителей)

 

Fig. 1. The comparison of the developed method speed-in-action with analogues
using 64 bit multipliers

 

 

 

 

 

Рис. 2. Сравнение быстродействия разработанного метода с аналогами
(с использованием 16-разрядных умножителей)

 

Fig. 2. The comparison of the developed method speed-in-action with analogues
using 16 bit multipliers

 

 

 

Таблица  1 Сравнение временной сложности разработанного метода и аналогов

Table  1 Time complexity comparison with analogues

 

Метод умножения /Multiplication method	Время выполнения, тактов /Execution time, clock cycles
Lei at al.	$\frac{n^2}{4}+2n+8$
Schulte ‒ Swarzlander	n2 + n + 12
Ishii	n2 + n + 7
Предложенный метод с масштабированием /Proposed method with scaling	pn+1 = 2α
Предложенный метод с масштабированием параллельный (k = n) /Proposed method with scaling parallel (k = n)	$\frac{5n+6}{2}$

 

 

Таблица  2 Сравнение быстродействия разработанного метода и аналогов (в тактах)

Table  2 The speed-in-action comparison (clock cycles) with analogues

 

Метод умножения / Multiplication method	Время выполнения, тактов / Execution time, clock cycles		Ускорение / Speed-up
	Аналоги / Analogues	Предложенный метод (k = n) / Proposed method (k = n)
1 024 бит / 1 024 bit
Lei at al.	104	43	2,4
Schulte ‒ Swarzlander	284	43	6,6
Ishii	279	43	6,4
2 048 бит / 2 048 bit
Lei at al.	328	83	4,0
Schulte ‒ Swarzlander	1 068	83	12,9
Ishii	1 063	83	12,8

 

 

Обсуждение и заключение

Разработаны новые быстрые методы умножения модулярных чисел с плавающей точкой; проведена оценка быстродействия разработанных методов; выполнено сравнение с работами других авторов. Предложенные методы в 2,4–4,0 раза быстрее конвейерного метода умножения и в 6,4–12,9 раз быстрее классических методов умножения.

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

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

Авторами рассмотрен алгоритм масштабирования коэффициентом, равным 2^a; при этом интервальная характеристика представлена в виде логарифмов по основанию 2. В то же время разработанный способ представления информации и выполнения операции умножения может быть использован и при других значениях основания логарифма. Так, разработанные методы могут быть применены для операций над числами вида M · 10^E. В таком случае коэффициент масштабирования будет равен 10^a. Данное преимущество может быть дополнительно использовано в задачах, критичных к ошибкам округления при вводе десятичной информации.

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

 

About the authors

Anastasia S. Korzhavina

Vyatka State University

Author for correspondence.
Email: as_korzhavina@vyatsu.ru
ORCID iD: 0000-0001-8270-2097
ResearcherId: S-1877-2018

Senior Lecturer, Chair of Electronic Computing Machines

Russian Federation, 36 Moskovskaya St., Kirov 610000

Vladimir S. Knyazkov

Penza State University; Vyatka State University

Email: kniazkov@list.ru
ORCID iD: 0000-0003-3820-6541
ResearcherId: Т-4089-2018

Chief Researcher, Research Institute of Applied and Fundamental Research, Professor, Chair of Electronic Computing Machines, D.Sc. (Engineering)

Russian Federation, 40 Krasnaya St., Penza 440026; 36 Moskovskaya St., Kirov 610000

References

Supplementary files