Конечно, лучше выбирать язык и компилятор более подходящий для научных расчетов. Но, если такой возможности нет, то возникает желание написать заплатку на ассемблере. Чтобы оценить выигрыш от такой заплатки я написал черновики с использованием SSE 4.2. При помощи директив условной компиляции (УК) подпрограммы могут быть скомпилированы:
(a) с проверкой аргументов на попадание в диапазон (переменная ExpRangeCheck) или без оной [далее кратко: R+ или R-];
(b) с модификацией без сохранения значений регистров xmm0—xmm3 (официальное соглашение: переменная УК DelphiRSC) или с модификацией регистров xmm0—xmm5;
(с) с масштабированием денормализованных результатов при помощи добавления подразумеваемой целой 1 и сдвига (переменная УК PsrlqScalling), или при помощи добавления фиксированного значения к показателю и умножения для компенсации этого добавления [далее кратко: PS и MS];
(d) c возможностью генерации исключительных ситуаций (переменная УК RaiseException) или без оной [далее кратко: RE+ или RE-].
Несмотря на то, что в документации Embarcadero указывается: без сохранения значений можно модифицировать только xmm0—xmm3, в библиотечных функциях они сами [разработчики Embarcadero] не сохраняют значения xmm4 и xmm5, а директивы .SAVENV xmm4 и .SAVENV xmm5 не приводят к формированию в прологе и эпилоге кода. Поэтому в приводимых ниже тестах программа была скомпилирована без определенной переменной DelphiRSC. Также в реализациях на асме не контролировались исключения «неточной операции» (“Precision Error”) и «денормализованный результат».
1. В прикреплении DPExp64.txt подпрограммы основаны на реализованном в cephes на С алгоритме, базирующемся на Паде аппроксимации.
Модуль DPExp64.pas содержит три подпрограммы:
pexpsd(x: Double): Double — для вычисления экспоненты от одного аргумента (это слегка доработанная и переименованная функция dexp из
приводившегося выше в теме модуля BMath64);
pexppsd(x1, x2: Double, out d: TDD) — процедура для вычисления за один вызов значений для двух аргументов, использующая при необходимости функцию pexpsd;
pexppd(x1, x2: Double, out d: TDD) — аналогичная процедура вычисления за один вызов двух значений, но не использующая pexpsd.
Сравнение скоростей выполнения проводилось на процессоре с архитектурой Sandy Bridge. В модификации R+, RE- время выполнения двух вызовов функции pexpsd составляет приблизительно 0.51, а время выполнения pexppd — 0.38 от времени выполнения двух вызовов стандартной функции
exp модуля system XE4. Т.е. выигрыш от вычисления за один вызов двух значений — мал.
(Детали сравнения)
В таблице приведены времена выполнения в единицах времени выполнения стандартной функции exp (XE4 64bit). В титульной строке указаны типы возвращаемых результатов:
н, н — оба нормализованные;
н, д — первый нормализованный, второй денормализованный;
д, н — первый денормализованный, второй нормализованный;
д, д — оба денормализованные.
Таб. 1 — Сравнение времен выполнения для приближенного вычисления экспоненты
Из первых трех строк таблицы видно: (даже при отключенной проверке попадания в диапазон) вычисление за один вызов процедуры сразу двух значений — даёт не много.
Из сравнения PS и MS видно: для случая «оба результата нормализованные» PS не хуже MS, а если хотя бы один из результатов денормализованный, то PS существенно быстрее.
Добавление возможности генерировать исключения (RE+) заметно снижает скорость выполнения (за счет создания пролога, эпилога и замены ret на «jmp на конец процедуры»).
В прикреплении DDExp64.txt две функции и процедура для вычисления экспоненты
с погрешностью не больше одного бита в последнем знаке мантиссы:
expsd_Tang (x: Double): Double — реализация алгоритма из работы Tang P.T.P. Table-Driven Implementation of the Exponential Function in IEEE Floating-Point Arithmetic // ACM Transactions on Mathematical Software, Vol. 15, Issue 2, P. 144–157 (1989);
expsd (x: Double): Double — модификация алгоритма SunPro, 1993 г.
exppd (x1, x2: Double; d: TDD) — процедура вычисления за один вызов двух экспонент, «векторизация» функции expsd.
Скорость expsd_Tang оказалась не выше, чем у expsd. В модификации R+, RE+ время выполнения двух вызовов expsd составляет приблизительно 0.69, а одного вызова exppd — 0.46 от времени выполнения двух вызовов стандартной функции
exp модуля system XE4. При включении возможности генерации исключений, выигрыш от использования «приближенной» процедуры pexppd вместо «точной» exppd становится мизерным.
(Детали сравнения)
Таб.2. — Сравнение времен выполнения для подпрограмм вычисления экспоненты с погрешностью не больше одного бита в последнем знаке мантиссы
Из первых трёх строк Таб.2 видно: времена выполнения expsd_Tang и expsd приблизительно одинаковы. Однако алгоритм Sun лучше «векторизуется». Поэтому процедура для вычисления экспоненты «для двух аргументов» написана только для алгоритма SunPro.
Как и п.1, использование MS приводит при получении денормализованных результатов к заметному уменьшению скорости выполнения.
expsd выполняется несколько медленней оригинальной версии SunPro при малых значениях аргумента, когда масштабирование не нужно (данные в таблице не приведены).
В качестве примера использования
pexpsd и
exppd был взят слегка изменённый пример «с двумя экспонентами» из сообщения
post840657.html#p840657 pi-314 (Он не вполне осмысленный, ну да ладно.) В конце сообщения приведен проект тестирования.
Для случая R-, PS, RE- время выполнения с использованием pexpsd составляло 0.35 от времени выполнения с использованием стандартной функции exp, а время выполнения с использованием exppd — 0.30. Т.е. скорость выполнения возросла в 2.88 и 3.34 раза соответственно.
Для случая R+, PS, RE+ время составило 0.35 и 0.31, т.е. добавление проверки на попадание в диапазон и возможности генерации исключений практически не замедлило выполнение.
Итого. «Векторизация», использование SSE4.2 вместо SSE2 (Embarcadero) и не отслеживание «неточный результат» и денормализованный операнд позволило значительно ускорить выполнение для большинства значений аргументов. По тексту подпрограмм встречаются места, в которых полезно использовать «неразрушающие возможности» AVX. Т.е. использование AVX позволит немного увеличить скорость выполнения даже при вычислении двух значений за один вызов процедуры.
Пожелания и вопрос. Помимо очевидного «Буду рад всем конструктивным замечаниям по тексту подпрограмм», есть пожелания и вопрос.
Я хотел выполнить сравнение и для Visual C++, но времени не оказалось. В тексте достаточно комментариев. На мой взгляд, его очень быстро можно перенести на другой 64-битный компилятор c SSE «соглашением» о передаче параметров и возврате результата функциями. Интересно было бы посмотреть на результаты сравнения скорости написанных функций и, например, используемых Visual C++.
Кроме двух указанных мною в теме, я не нашел полезных ссылок на посвященные вычислению exp статьи с алгоритмами удобными для реализации с использованием SSE или AVX . Если такие ссылки имеются, то приведите, пожалуйста, в ветке.
64-битный компилятор Embarcadero XE4 неправильно генерирует код для перемещения учетверённого слова из xmm-регистра в регистр общего назначения и наоборот. Факт известный (см., например,
QualityCentral/Delphi-BCB/Compiler/Delphi/BASM/Report 112094). Это обоснованная особенность или откровенный баг? Если баг, то исправили ли его в последней версии? И как обстоят дела с этим в других компиляторах?
program Tst;
uses DDExp64, DPExp64, sysutils;
{$APPTYPE CONSOLE}
function stdTst: Double;
const
Delta: Double = 0.0000001;
var
x: Double;
Sum: Double;
begin
x := -10.0;
Sum := 0.0;
repeat
if X > 10.0 then Break;
Sum := Sum + exp(x)*Delta/exp(10*x);
x := x + Delta;
until False;
Result:= Sum;
end;
function pexpTst: Double;
const
Delta: Double = 0.0000001;
var
x: Double;
Sum: Double;
begin
x := -10.0;
Sum := 0.0;
repeat
if x > 10.0 then Break;
Sum := Sum + pexpsd(x)*Delta/pexpsd(10*x);
x := x + Delta;
until False;
Result:= Sum;
end;
function exppdTst: Double;
const
Delta: Double = 0.0000001;
var
x: Double;
Sum: Double;
d: DDExp64.TDD;
begin
x := -10.0;
Sum := 0.0;
repeat
if x > 10.0 then Break;
exppd(x, 10*x, d);
Sum := Sum + d.Lo*Delta/d.Hi;
x := x + Delta;
until False;
Result:= Sum;
end;
var Time1, Time2, d1, d2, d3: TDateTime;
f: text;
begin
Assign(f, 'PITst.txt');
Append(f);
Writeln(f, SS);
Time1:= Time;
StdTst;
Time2:= Time;
d1:= Time2 - Time1;
Writeln(f, ' S=', d1*24 * 60 * 60:10:3);
Time1:= Time;
pexpTst;
Time2:= Time;
d2:= Time2 - Time1;
Writeln(f, ' P=', d2*24 * 60 * 60:10:3, ' P/S =', d2/d1:4:2, ' S/P =', d1/d2:4:2);
Time1:= Time;
exppdTst;
Time2:= Time;
d3:= Time2 - Time1;
Writeln(f, ' D=', d3*24 * 60 * 60:10:3, ' D/S =', d3/d1:4:2, ' S/D =', d1/d3:4:2);
close(f);
end.
![$[6/4]_S$](https://dxdy-02.korotkov.co.uk/f/1/3/0/130ec5d84f577c65984f15c83a01557e82.png "$[6/4]_S$")
имеет вид![$$[6/4]_S =\frac Q P = \frac {2+ \frac 5 {22} x^2 + \frac 1 {396} x^4 + \frac 1 {332640} x^6} {1 + \frac 1 {33} x^2 + \frac 1 {7920} x^4}$$](https://dxdy-01.korotkov.co.uk/f/8/f/4/8f4ec1d616ff1bf3a572e862eaf0918782.png "$$[6/4]_S =\frac Q P = \frac {2+ \frac 5 {22} x^2 + \frac 1 {396} x^4 + \frac 1 {332640} x^6} {1 + \frac 1 {33} x^2 + \frac 1 {7920} x^4}$$")
В формате с плавающей точкой значения коэффициентов будут
генерировались 1000 миллионов значений 
и вычислялись: функция dexp, использующая приведенные выше в сообщении коэффициенты, и стандартная функция exp, использующая FPU. Максимальная разница в возвращаемых значениях составляла 2 двоичных знака, максимальная погрешность — приближенно 
. Такое же расхождение в количестве двоичных знаков (и максимальная погрешность) было и для значений из cephes (но при другом значении аргумента).