2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему На страницу Пред.  1, 2, 3, 4
 
 Re: Быстрое вычисление функции exp(x)
Сообщение27.03.2014, 06:19 
Заслуженный участник


12/07/07
4135
Донецк, Украина
Pphantom, приведите, пожалуйста, текст используемой IFC функции exp. Если получится, я его вставлю в сообщение post841313.html#p841313.

 Профиль  
                  
 
 Re: Быстрое вычисление функции exp(x)
Сообщение27.03.2014, 17:14 
Заслуженный участник


12/07/07
4135
Донецк, Украина
Посмотрел текcт ассемблерного листинга gcc в сообщении mustitz. В сообщении 840724 «с ходу» оптимизировал не очень: встроенные константы fpu желательно было хранить в стеке fpu, а не загружать в цикле. Кроме того, если допустить использование не только инструкций 386/387, но и Pentium II (если не путаю), то комбинацию
fucom
fstsw ax
sahf

следует заменить на
fcomi.

В результате к тексту в сообщении 840724 "добавилась" новая функция MainA4
код: [ скачать ] [ спрятать ]
Используется синтаксис ASM
 procedure MainA4();
var
  Sum: Extended;
begin
 asm
  fld    Xmax
  fld    Delta
  fldl2e
  fld1
  fldz
  fld    X0
@1:
  fld    st(0)
  fmul   st, st(4)     { y := x*log2e;      }
  fld    st(0)         { i := round(y);     }
  frndint
  fsub   st(1),   st   { f := y - i;        }
  fxch   st(1)         { z := 2**f          }
  f2xm1
  fadd  st(0), st(4)
  fscale               { result := z * 2**i }
  fstp   st(1)
  faddp  st(2),  st    {Sum := Sum + Exp(X);}
  fadd   st,     st(4) {X := X + Delta}
  fcomi  st,     st(5) {X < Xmax ?}
  jbe   @1             {until X < Xmax}
  fstp   st(2)
  fstp   st(1)
  fstp   st(1)
  fstp   st(1)
  fstp   st(1)
  fstp   Sum
 end;
 WriteLn(Sum);
end;

У меня получилось, что разница во времени выполнения MainA2 и MainA4 около 3%. Выигрыш небольшой, но заметный. Преимущественно за счет хранения констант в стеке fpu.
mustitz, спасибо!

upd Компилировал для выполнения уже не в Delphi 5, а в Delphi 7.
[Встроенный ассемблер Delphi 5 поддерживает инструкции до Pentium I включительно, но не выше. В частности, fcomi он не поддерживает.]

 Профиль  
                  
 
 Re: Быстрое вычисление функции exp(x)
Сообщение27.03.2014, 20:24 
Супермодератор
Аватара пользователя


09/05/12
24060
Кронштадт
GAA в сообщении #841428 писал(а):
Pphantom, приведите, пожалуйста, текст используемой IFC функции exp. Если получится, я его вставлю в сообщение post841313.html#p841313 .

Нет, не получится. Или, во всяком случае, сходу не получилось. А в результате дизассемблирования сам черт ногу сломит...

 Профиль  
                  
 
 Re: Быстрое вычисление функции exp(x)
Сообщение29.03.2014, 17:51 
Заслуженный участник


12/07/07
4135
Донецк, Украина
Нашёл по теме в сети (возможно, будет интересно):
1. Презентация “The VDT mathematical library: a modern reimplementation of Cephes” содержит сравнение скорости и точности вычисления элементарных функций различными библиотеками (VDT, SVML & VC). SVML (Short Vector Math Library, Intel)— самая быстрая, но без открытого кода и немного менее точная по сравнению с VDT.
Upd. Статья Hauth T., Innocente V., Piparo D. Development and Evaluation of Vectorised and Multi-Core Event Reconstruction Algorithms within the CMS Software Framework// J. Phys.: Conf. Ser. 396 052065.
2. Есть исходные тексты VDT на Си. (Прямые ссылки: Files\File list\math\vdt\include, Files\File list\math\vdt\tests.) [Но без SSE и AVX, или я не нашёл.]
На ассемблере На Си, но с SSE2 [тоже на базе Cephes, но не VDT]: Julien Pommier “Simple SSE and SSE2 (and now NEON) optimized sin, cos, log and exp”.

[Добавлено 15.04.2014]
GAA в сообщении #841142 писал(а):
Я впервые залез в RTL FPC, поэтому мог не туда посмотреть.
Исходя из текста в \inc\systemh.inc “win64 doesn't support the legacy fpu”, для Win64 посмотрел не туда.

В system.inc находим
Код:
{$ifndef FPUNONE}
{$ifdef FPC_USE_LIBC}
{ Include libc versions }
{$i cgenmath.inc}
{$endif FPC_USE_LIBC}
{ Include processor specific routines }
{$I math.inc}
{ Include generic version }
{$I genmath.inc}
{$endif}

В genmath.inc находим
код: [ скачать ] [ спрятать ]
Используется синтаксис Delphi
{$ifndef FPC_SYSTEM_HAS_EXP}
{$ifdef SUPPORT_DOUBLE}
    {
     This code was translated from uclib code, the original code
     had the following copyright notice:

     *
     * ====================================================
     * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
     *
     * Developed at SunPro, a Sun Microsystems, Inc. business.
     * Permission to use, copy, modify, and distribute this
     * software is freely granted, provided that this notice
     * is preserved.
     * ====================================================
     *}


    {*
     * Returns the exponential of x.
     *
     * Method
     *   1. Argument reduction:
     *      Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
     *  Given x, find r and integer k such that
     *
     *               x = k*ln2 + r,  |r| <= 0.5*ln2.
     *
     *      Here r will be represented as r = hi-lo for better
     *  accuracy.
     *
     *   2. Approximation of exp(r) by a special rational function on
     *  the interval [0,0.34658]:
     *  Write
     *      R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
     *      We use a special Reme algorithm on [0,0.34658] to generate
     *  a polynomial of degree 5 to approximate R. The maximum error
     *  of this polynomial approximation is bounded by 2**-59. In
     *  other words,
     *      R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
     *      (where z=r*r, and the values of P1 to P5 are listed below)
     *  and
     *      |                  5          |     -59
     *      | 2.0+P1*z+...+P5*z   -  R(z) | <= 2
     *      |                             |
     *  The computation of exp(r) thus becomes
     *                     2*r
     *      exp(r) = 1 + -------
     *                    R - r
     *                         r*R1(r)
     *             = 1 + r + ----------- (for better accuracy)
     *                        2 - R1(r)
     *  where
                         2       4             10
     *     R1(r) = r - (P1*r  + P2*r  + ... + P5*r   ).
     *
     *   3. Scale back to obtain exp(x):
     *  From step 1, we have
     *     exp(x) = 2^k * exp(r)
     *
     * Special cases:
     *  exp(INF) is INF, exp(NaN) is NaN;
     *  exp(-INF) is 0, and
     *  for finite argument, only exp(0)=1 is exact.
     *
     * Accuracy:
     *  according to an error analysis, the error is always less than
     *  1 ulp (unit in the last place).
     *
     * Misc. info.
     *  For IEEE double
     *      if x >  7.09782712893383973096e+02 then exp(x) overflow
     *      if x < -7.45133219101941108420e+02 then exp(x) underflow
     *
     * Constants:
     * The hexadecimal values are the intended ones for the following
     * constants. The decimal values may be used, provided that the
     * compiler will convert from decimal to binary accurately enough
     * to produce the hexadecimal values shown.
     *
    }

    function fpc_exp_real(d: ValReal):ValReal;compilerproc;
      end;

Текст функции

{$else SUPPORT_DOUBLE}

    function fpc_exp_real(d: ValReal):ValReal;compilerproc;
    {*****************************************************************}
    { Exponential Function                                            }
    {*****************************************************************}
    {                                                                 }
    { SYNOPSIS:                                                       }
    {                                                                 }
    { double x, y, exp();                                             }
    {                                                                 }
    { y = exp( x );                                                   }
    {                                                                 }
    { DESCRIPTION:                                                    }
    {                                                                 }
    { Returns e (2.71828...) raised to the x power.                   }
    {                                                                 }
    { Range reduction is accomplished by separating the argument      }
    { into an integer k and fraction f such that                      }
    {                                                                 }
    {     x    k  f                                                   }
    {    e  = 2  e.                                                   }
    {                                                                 }
    { A Pade' form of degree 2/3 is used to approximate exp(f)- 1     }
    { in the basic range [-0.5 ln 2, 0.5 ln 2].                       }
    {*****************************************************************}
    const  P : TabCoef = (
           1.26183092834458542160E-4,
           3.02996887658430129200E-2,
           1.00000000000000000000E0, 0, 0, 0, 0);
           Q : TabCoef = (
           3.00227947279887615146E-6,
           2.52453653553222894311E-3,
           2.27266044198352679519E-1,
           2.00000000000000000005E0, 0 ,0 ,0);

           C1 = 6.9335937500000000000E-1;
            C2 = 2.1219444005469058277E-4;
    var n : Integer;
        px, qx, xx : Real;
    begin
      if( d > MAXLOG) then
          float_raise(float_flag_overflow)
      else
      if( d < MINLOG ) then
      begin
        float_raise(float_flag_underflow);
      end
      else
      begin

     { Express e**x = e**g 2**n }
     {   = e**g e**( n loge(2) ) }
     {   = e**( g + n loge(2) )  }

        px := d * LOG2E;
        qx := Trunc( px + 0.5 ); { Trunc() truncates toward -infinity. }
        n  := Trunc(qx);
        d  := d - qx * C1;
        d  := d + qx * C2;

      { rational approximation for exponential  }
      { of the fractional part: }
      { e**x - 1  =  2x P(x**2)/( Q(x**2) - P(x**2) )  }
        xx := d * d;
        px := d * polevl( xx, P, 2 );
        d  :=  px/( polevl( xx, Q, 3 ) - px );
        d  := ldexp( d, 1 );
        d  :=  d + 1.0;
        d  := ldexp( d, n );
        result := d;
      end;
    end;
{$endif SUPPORT_DOUBLE}

 Профиль  
                  
 
 Re: Быстрое вычисление функции exp(x)
Сообщение07.04.2014, 18:37 
Заслуженный участник


12/07/07
4135
Донецк, Украина
Аппроксимация Паде имеет вид отношения полиномов, которые могут вычисляться параллельно. Это наводит на мысль о плодотворности использования SSE при вычислении экспоненты с удвоенной точностью. Для проверки этого предположения и эффективности алгоритма, реализованного в cephes, я попробовал перевести функцию exp на ассемблер и сравнить скорость выполнения стандартной функции exp Delphi и написанной на ассемблере (Текст модуля для w32: Bmath.pas; проект, при помощи которого сравнивались времена выполнения: BMathTst.dpr; см. в конце сообщения в прикреплённом файле BMath.txt).

Результат тестирования: в Win32 использующая SSE процедура более чем в 1.5 раза медленней стандартной функции (FPU).

(Детали)

Для сравнения были использованы три компилятора:
    * Borland Delphi 7 (B7) [понимает инструкции до SSE2 включительно];
    * CodeGear Delphi for MS Windows 2007 (C2007) [до SSE3 включительно];
    * Embarcadero XE4, 21013 (XE4) [понимает SSE4].
При вызове функции exp (как и в случае некоторых других подпрограмм модуля system) B7 передаёт вещественный параметр через вершину стека FPU; а C2007 и XE4 — через программный стек (как и для любых подпрограмм). За исключением входа и выхода, код exp одинаков в библиотеках этих компиляторов. Т.е. выполнение system.exp B7 будет самым быстрым. Поэтому время выполнения цикла с этой (system.exp) функцией в скомпилированной B7 программе было выбрано за единицу.

Тестирование выполнялось на процессорах с архитектурами Penryn и Sandy Bridge.

Результаты для Penryn (XP, тестировался только код для SSE2)
\small \begin{array}{|c | c | c | c | c |} 
\hline
\text{ C2007 станд. } & \text{ B7, C2007 | SSE2} & \text{B7 SSE2 FPU Scale} \\
\hline
1.02 & 1.53 & 1.75 \\ 
\hline
\end{array}
В столбце B7 SSE2 FPU Scale приведен результат для масштабирования результата вычисления exp при помощи FPU (компиляция без определённой переменной CPU_SCALE). Изменение типа вызова функции приводит к её незначительному замедлению (1.02).
Удивительно, но использование FPU для масштабирования приносит существенное замедление.

Результаты для Sandy Bridge (Win7, 64-битная)
\small \begin{array}{|c | c | c | c | c |} 
\hline
\text{C2007;XE4 | ст. } & \text{B7;C2007;XE4 | SSE2} & \text{XE4 SSE4} & \text{C2007 SSE2FPU Scale} & \text{XE4 SSE4 FPU Scale }\\
\hline
\text{1.2} & \text{2.5} & \text{2.1} & \text{2.4} & \text{2.0} \\ 
\hline
\end{array}

Недостатками компиляторов B7 и C2007 являются:
1) невозможность задать 16-битное выравнивание данных (что приводит к необходимости вместо непосредственного смещения использовать адресацию по базе);
2) невозможность указать компилятору: в каком регистре возвращается результат, через какой регистр передаётся вещественный параметр.

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

Невозможность приказать компилятору передавать параметр через xmm регистр тоже приводит к увеличению времени выполнения на проценты (менее чем в 1.1 раза). Также не приводит к существенным задержкам возврат через стек FPU (заодно можно выполнить масштабирование средствами fpu, или масштабировать переменную в памяти средствами cpu по ходу перебрасывании в fpu).

Другими словами, подгоняя под другие компиляторы эту ассемблерную процедуру, мы бы выиграли проценты.

Замечательно то, что с переходом к Sandy Bridge, то ли FPU улучшили, то ли SSE ухудшили. Использование FPU выглядит заметно предпочтительней.

В случае с компиляцией 64-битных программ ситуация изменяется. Ряд компиляторов, в частности XE4, используют SSE для вещественной арифметики. Extended (long double) считается эквивалентным double. В подпрограмму вещественный параметр передаётся через xmm0, и через него же возвращается.

В файле BMath64.txt приведен модуль BMath64 (с функцией для 64-битного режима XE4) и проект BMathTst64 для тестирования. Функции проекта написаны на ассемблере для минимизации влияния окружения вызова функции на оценку скорости вычисления экспоненты. Тестирование выполнялось только на Sandy Bridge. Для того чтобы оценить плодотворность одновременного вычисления числителя и знаменателя в аппроксимации Паде тестировалось две модификации функции: в одной числитель и знаменатель вычислялись параллельно, в приводимой ниже таблице SSE P (от packed), а в другой — последовательно SSE S (от scalar).

Результаты для Sandy Bridge, отношения времён выполнения (Win7, 64-битная)
\small \begin{array}{| c | c | c |} 
\hline
\text{XE4/FPU} & \text{SSE P/FPU} & \text{SSE S / SSE P} \\
\hline
\text{ 1.01} & \text{0.56} & \text{ 1.06} \\ 
\hline
\end{array}
Основанная на SSE и используемая по умолчанию в XE4 функция имеет приблизительно такую же скорость, что и стандартная функция, использующая FPU. Функция, основанная на аппроксимации Паде, заметно быстрее. Параллельное вычисление почти не даёт выигрыша во времени по сравнению со скалярным вариантом: слишком много скалярных вычислений + дополнительные действия для выполнения параллельного вычисления.

Upd. 16.04.2014 Я слегка оптимизировал код, добавил комментарии и скалярный вариант функции, а также отредактировал конец этого сообщения.

Upd2 4.05.2014 Для win64 (Embarcadero XE4) убрал сохранение регистров (см. начало post851310 ниже в теме), добавил вычисление функции от больших по модулю отрицательных значений аргумента (когда результат денормализованное число), а также генерацию исключительных ситуаций (Underflow и Overflow). В случае получения денормализованного результата время выполнения написанной реализации функция составляет приблизительно 0.10 от времени выполнения использующей FPU функции, а время выполнения стандартной SSE функции XE4 — приблизительно 0.6 от использующей FPU функции.
Медленное вычисление экспоненты (в случае больших отрицательных значений аргумента) функцией использующей FPU связано с медленным выполнением инструкции fscale для таких значений аргумента. Это очевидно, но, на всякий случай, для непосредственной проверки я написал простенькую программу (fscale.dpr) под Win32. При значениях аргумента -1070 или -1025 время выполнения приблизительно в 19 раз больше чем при 1, -900 или -1025.
код: [ скачать ] [ спрятать ]
Используется синтаксис Delphi
program fscale; {Win32}
uses sysutils;

const
 n1:   Integer   = 1;
 n2:   Integer   = -900;
 n3:   Integer   = -1020;
 n4:   Integer   = -1025;
 n5:   Integer   = -1070;
 rep: LongWord = 500000000;
var
 d: Double;

procedure scale;
asm
        fldl2e
        fscale
        fstp     st(1)
end;

Procedure MainF(n: Integer); stdcall;
var i: LongWord;
 CW: Word;
asm
  mov    ecx,   [rep]
  mov    [i],   ecx
  mov    CW, 1372h
  fldcw  CW
@loop:
  fild   dword ptr  [n]
  call   scale
  fstp   qword ptr  [d]
  dec    [i]
  jnz    @loop
end;

var Time1, Time2, d0, d1, d2, d3: TDateTime;
  f: text;

begin
  Assign(f, 'tstscale.txt');
  Append(f);
  Writeln(f);
  Writeln(f, 'rep=', rep);

  Time1:= Time;
  MainF(n1);
  Time2:= Time;
  d0:= Time2 - Time1;
  Writeln(f, 'n = ', n1:7, d0*24 * 60 * 60:10:4);

  Time1:= Time;
  MainF(n2);
  Time2:= Time;
  d1:= Time2 - Time1;
  Writeln(f, 'n = ', n2:7, d1*24 * 60 * 60:10:4);
  if d0 <> 0 then Writeln(f, 'FPU2/FPU1=', d1/d0:7:3);

  Time1:= Time;
  MainF(n3);
  Time2:= Time;
  d2:= Time2 - Time1;
  Writeln(f, 'n = ', n3:7, d2*24 * 60 * 60:10:4);
  if d0 <> 0 then Writeln(f, 'FPU3/FPU1=', d2/d0:7:3);

  Time1:= Time;
  MainF(n4);
  Time2:= Time;
  d3:= Time2 - Time1;
  Writeln(f, 'n = ', n4:7,  d3*24 * 60 * 60:10:4);
  if d2 <> 0 then Writeln(f, 'FPU4/FPU1=', d3/d2:7:3);

  Time1:= Time;
  MainF(n5);
  Time2:= Time;
  d3:= Time2 - Time1;
  Writeln(f, 'n = ', n5:7,  d3*24 * 60 * 60:10:4);
  if d2 <> 0 then Writeln(f, 'FPU4/FPU1=', d3/d2:7:3);

  close(f);
end.


Вложения:
BMath64.txt [16.08 Кб]
Скачиваний: 169
BMath.txt [29.43 Кб]
Скачиваний: 117
 Профиль  
                  
 
 Re: Быстрое вычисление функции exp(x)
Сообщение07.04.2014, 19:08 


27/02/14
22
2 GAA:
Вот это research! :appl:

Большая просьба от тех, кто в asm полный ноль - выложите, пожалуйста, Ваш самый быстрый вариант на asm для случая единичного вычисления exp(), без нахождения суммы экспонент внутри функции.
Позаимствую для решаемой задачи! :-)

 Профиль  
                  
 
 Re: Быстрое вычисление функции exp(x)
Сообщение17.04.2014, 12:39 
Заслуженный участник


12/07/07
4135
Донецк, Украина
Я добавил в сообщение post842741.html#p842741 исходный текст функций из RTL FPC [оказалось, что «всё в нём есть» :)] и отредактировал конец своего предыдущего сообщения. Попытка оптимизация не привела к заметным изменениям результата.

На мой взгляд, при сравнении скорости генерируемого различными компиляторами кода нужно учитывать следующие моменты.

1. Если нужны особо точные вычисления, то следует использовать 80-битный вещественный тип. Этот тип не поддерживается SSE. Следовательно, сравнивать имеет смысл только скорости использующего FPU кода. Delphi и, по-видимому, FPC генерируют такой код для win32, но не для win64. При этом необходимо удостовериться, что в FPU установлена точность вычислений 80 бит.

2. Если для расчётов достаточно 64-битного типа, то имеет смысл сравнивать код, который обеспечивает точность Double. В этом случае возможен как вариант с SSE, так и вариант FPU с установленной точностью 64-бита.
[Я для Win64 (Embarcadero XE4) специально сравнивал плохо сравнимые вещи: FPU с 80-битной точностью вычислений и SSE. Даже при таком «нечестном» сравнении FPU не проиграл (см. первый столбец табл. для Win64 в предыдущем сообщении). Можно предположить: в XE4 неудачный вариант (по умолчанию) функции для вычисления экспоненты (SSE). Но я сомневаюсь, что он очень уж сильно проигрывает функциям в библиотеках других компиляторов. Хотя немного проигрывать он вполне может. XE4, например, не поддерживает AVX. Да и вообще Delphi ещё со времён Borland не очень активно учитывает новые инструкции, особенно в функциях system.]

3. Если обеспечивать точность Double нет необходимости, то, очевидно в этом случае есть больше возможностей для манёвра и возможно заметное ускорение вычислений. В этом случае при сравнении хорошо бы указывать погрешность вычислений. Иначе возможно сравнение «несравнимого».

 Профиль  
                  
 
 Re: Быстрое вычисление функции exp(x)
Сообщение18.04.2014, 08:34 


27/02/14
22
2 GAA: большое спасибо за исходники!

 Профиль  
                  
 
 Re: Быстрое вычисление функции exp(x)
Сообщение18.04.2014, 14:35 
Заслуженный участник


12/07/07
4135
Донецк, Украина
С сохранением и восстановлением xmm регистров для Embarcadero XE4 в функциях на ассемблере я был не уверен. Сегодня нашел
In line with the x64 Application Binary Interface (ABI), the contents of the following registers must be preserved and restored within inline assembly functions: R12, R13, R14, R15, RDI, RSI, RBX, RBP, RSP, XMM4, XMM5, XMM6, XMM7, XMM8, XMM8, XMM9, XMM10, XMM11, XMM12, XMM13, XMM14, and XMM15.
Т.е., если вдруг возникнет желание применять выложенную выше функцию для Win64, то надо проверить по используемой версии компилятора: скорее всего сохранение и восстановление xmm-регистров можно удалить. (На скорость вычислений это влияет крайне мало. Результаты тестов почти не изменятся.)

Там еще что-нибудь, возможно, надо доделать. Например, масштабирование при помощи SSE (а не CPU), генерацию исключительных ситуаций (это зависит от компилятора, и вообще, привязка к конкретному компилятору дело тонкое). В идеале скомпилировать на внешнем ассемблере и затем прилинковывать. Но тут уже начинается борьба за крохи. И тестировать надо, увы, старательно…

Добавлено 25.04.14

Как получены в cephes или Free Pascal выражения для коэффициентов многочленов P и Q, я не знаю, но можно было воспользоваться таким стандартным подходом.
Из определения гиперболического котангенса вытекает
$\exp (x) = \frac{S+x}{S-x}= 1+ \frac{2x}{S-x},$
где $S = x \cth (x/2)$
Форма Паде $[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}$$В формате с плавающей точкой значения коэффициентов будут
Q[0] = 2.0
Q[1] = 2.272727272727272727273e-1
Q[2] = 2.525252525252525252525e-3
Q[3] = 3.006253006253006253006e-6
P[0] = 1.0
P[1] = 3.03030303030303030303e-2
P[2] = 1.26262626262626262626e-4
P[3] = 0.0

Значения этих коэффициентов несколько отличаются от используемых в cephes или Free Pascal значений. Для проверки насколько они могут быть плохи, я выполнил следующее простенькое тестирование. Равномерно на $[-\ln(2)/2, \ln(2)/2)$ генерировались 1000 миллионов значений $x$ и вычислялись: функция dexp, использующая приведенные выше в сообщении коэффициенты, и стандартная функция exp, использующая FPU. Максимальная разница в возвращаемых значениях составляла 2 двоичных знака, максимальная погрешность — приближенно $3\cdot 10^{-16}$. Такое же расхождение в количестве двоичных знаков (и максимальная погрешность) было и для значений из cephes (но при другом значении аргумента).

 Профиль  
                  
 
 Re: Быстрое вычисление функции exp(x)
Сообщение04.10.2014, 22:28 
Заслуженный участник


12/07/07
4135
Донецк, Украина
Конечно, лучше выбирать язык и компилятор более подходящий для научных расчетов. Но, если такой возможности нет, то возникает желание написать заплатку на ассемблере. Чтобы оценить выигрыш от такой заплатки я написал черновики с использованием 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 — Сравнение времен выполнения для приближенного вычисления экспоненты
\small \begin{tabular}{|c | c | c | l | c |c |c |c |} 
\hline
 &  &  &  & н, н & н, д & д, н & д, д \\
\hline
       &       &    & pexpsd  & 0.48 & 0.20 & <- & 0.13\\
       &       & PS & pexppsd & 0.42 & 0.18 & 0.17 & 0.11 \\
       & R-  &    & pexppd  & 0.36 & 0.16 & 0.16 & 0.10 \\
\cline{3-8}
       &      &     & pexpsd   & 0.49 & 0.85&  <- & 0.98\\
       &      & MS & pexppsd & 0.44 & 0.80 & <- & 0.50 \\
RE-  &      &     & pexppd   & 0.38 & 0.77 & <- & 0.48 \\
\cline{2-8}
       &      &     & pexpsd    & 0.51 & 0.22 & <- & 0.14 \\
       &      &  PS & pexppsd  & 0.42 & 0.18 & 0.18 & 0.11 \\
       & R+&     & pexppd    & 0.38 & 0.17 & 0.16 & 0.11\\
\cline{3-8}
       &      &     & pexppd   & 0.51 &  0.85 & <-   & 0.99 \\
       &      & MS & pexppd   & 0.47 & 0.81 & 0.80 & 0.50\\
       &      &      & pexppd   & 0.43 & 0.79 & 0.78 & 0.49\\
\hline
RE+ &R+ & PS   & pexppd   & 0.43 & 0.18 & 0.18 & 0.11 \\
        &      & MS & pexppd   & 0.47 & 0.82 &  0.81 & 0.51\\
\hline
\end{tabular}
Из первых трех строк таблицы видно: (даже при отключенной проверке попадания в диапазон) вычисление за один вызов процедуры сразу двух значений — даёт не много.
Из сравнения PS и MS видно: для случая «оба результата нормализованные» PS не хуже MS, а если хотя бы один из результатов денормализованный, то PS существенно быстрее.
Добавление возможности генерировать исключения (RE+) заметно снижает скорость выполнения (за счет создания пролога, эпилога и замены ret на «jmp на конец процедуры»).

2. В прикреплении 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. — Сравнение времен выполнения для подпрограмм вычисления экспоненты с погрешностью не больше одного бита в последнем знаке мантиссы
\small \begin{tabular}{|c | c | c | l | c |c |c |c |} 
\hline
 &  &  &  & н, н & н, д & д, н & д, д \\
\hline
       &       & & expsd\_Tang  & 0.63 & 0.26 & <- & 0.16\\
       &       & PS & expsd &  0.62 & 0.25 & <-  &  0.17\\
       & R-  &    & exppd  & 0.40 & 0.19 &0.18 & 0.12 \\
\cline{3-8}
       &      & & expsd\_Tang & 0.63  & 0.87 &  <- & 0.99\\
       &      & MS & expsd    & 0.61 &  0.90 & <- & 1.04\\
RE-  &      &     & exppd     & 0.44 & 0.85 & <- & 0.53 \\
\cline{2-8}
       &      & & expsd\_Tang & 0.66 & 0.28 & <- & 0.18 \\
       & R+&  PS & expsd     & 0.66 & 0.27 & <- & 0.18 \\
       &      &     & exppd       & 0.45 & 0.20 & 0.20 & 0.13\\
\hline
        &  & & expsd\_Tang & 0.69 & 0.28 & <- & 0.19 \\
RE+  &R+ & PS & expsd      & 0.69 & 0.28 &  <- & 0.19\\
        &      &      & exppd      & 0.46 & 0.20 &  0.20 & 0.13\\
\hline
\end{tabular}
Из первых трёх строк Таб.2 видно: времена выполнения expsd_Tang и expsd приблизительно одинаковы. Однако алгоритм Sun лучше «векторизуется». Поэтому процедура для вычисления экспоненты «для двух аргументов» написана только для алгоритма SunPro.
Как и п.1, использование MS приводит при получении денормализованных результатов к заметному уменьшению скорости выполнения.
expsd выполняется несколько медленней оригинальной версии SunPro при малых значениях аргумента, когда масштабирование не нужно (данные в таблице не приведены).

3. В качестве примера использования 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). Это обоснованная особенность или откровенный баг? Если баг, то исправили ли его в последней версии? И как обстоят дела с этим в других компиляторах?

код: [ скачать ] [ спрятать ]
Используется синтаксис Delphi
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.


Вложения:
Комментарий к файлу: Модуль DPExp64.pas и проект для тестирования скоростей выполнения (Embarcadero Delphi XE4 64-bit)
DPExp64.txt [53.15 Кб]
Скачиваний: 140
Комментарий к файлу: Модуль DDExp64.pas и проект для тестирования скоростей выполнения (Embarcadero Delphi XE4 64-bit)
DDExp64.txt [45.64 Кб]
Скачиваний: 120
 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 55 ]  На страницу Пред.  1, 2, 3, 4

Модераторы: maxal, Toucan, PAV, Karan, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group