Быстрое вычисление функции exp(x) [Free Pascal, asm]

[pi-314]

Всем привет,

сразу оговорюсь - не программист, не математик, а прикладной гидромеханик.
Задача численного решения системы уравнений мат.физики столкнула с программированием на Free Pascal (в среде Lazarus).

Код написал, численная схема работает, но задача ускорения выполнения программы осложнена медленной работой встроенной в Free Pascal функции для вычисления exp(x).

Вопрос гуглил, для Pascal решений не нашел, просьба помочь.

[Sender]

Каковы требуемые диапазон значений аргумента и точность вычислений?

[pi-314]

Sender в сообщении #831048 писал(а):

Каковы требуемые диапазон значений аргумента и точность вычислений?

диапазон значений аргумента от -1e-10 до 100; точность 1e-10

[_Ivana]

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

[pi-314]

_Ivana в сообщении #831065 писал(а):

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

тип переменной extended разве не решает эти трудности ?

[kp9r4d]

pi-314 в сообщении #831073 писал(а):

тип переменной extended разве не решает эти трудности ?

Под тип extended в FPC выделяется 10 байт. 10 байт — это 80 бит, в самой грубой оценке (не учитывая биты, которые тратятся на знак и на экспоненту) можно хранить $\log_{10}(2^{80}) \sim 24$ значащих цифр. Вам же надо $\log_{10}(e^{100}) \sim 43$ значащие цифры и ещё 10 цифр после запятой, то есть 53 значащие цифры. Не влезает явно.
Какую задачу вы решаете? Может вам больше матпакеты вроде Wolfram Mathematica подойдут?

[rockclimber]

kp9r4d в сообщении #831113 писал(а):

pi-314 в сообщении #831073 писал(а):

тип переменной extended разве не решает эти трудности ?

Под тип extended в FPC выделяется 10 байт. 10 байт — это 80 бит, в самой грубой оценке (не учитывая биты, которые тратятся на знак и на экспоненту) можно хранить $\log_{10}(2^{80}) \sim 24$ значащих цифр. Вам же надо $\log_{10}(e^{100}) \sim 43$ значащие цифры и ещё 10 цифр после запятой, то есть 53 значащие цифры. Не влезает явно.
Какую задачу вы решаете? Может вам больше матпакеты вроде Wolfram Mathematica подойдут?

53 значащие цифры для гидромеханика-практика? Вы, скорее всего, неправильно поняли запись

pi-314 в сообщении #831053 писал(а):

диапазон значений аргумента от -1e-10 до 100; точность 1e-10

которая обычно означает "инженерный формат" записи числа: "1e-10" означает $10^{-10}$, то есть автору скорее всего надо 10 - 12 значащих цифр.

pi-314 в сообщении #831041 писал(а):

Код написал, численная схема работает, но задача ускорения выполнения программы осложнена медленной работой встроенной в Free Pascal функции для вычисления exp(x).

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

[kp9r4d]

rockclimber в сообщении #831118 писал(а):

которая обычно означает "инженерный формат" записи числа: "1e-10" означает $10^{-10}$, то есть автору скорее всего надо 10 - 12 значащих цифр.

«Диапазон значений аргумента» означает, что $x$ может принимать значения до $10^{-10}$ если уж так. Хотя да, это уже совсем неадекватно, скорее всего имелся в виду как раз диапазон значений самой экспоненты. Тогда да, с экстендедом никаких проблем не возникнет.

[pi-314]

Цитата:

53 значащие цифры для гидромеханика-практика? Вы, скорее всего, неправильно поняли запись , то есть автору скорее всего надо 10 - 12 значащих цифр.

Вы правы, нужны 10-12 значащих цифр)

Цитата:

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

Профайлера под Lazarus (Win) не нашел, поэтому проверял кустарно через :

procedure log_on(t:integer); {входная процедура счетчика продолжительности выполнения участков кода}
begin
tt2[t] := GetTickCount;
end;

procedure log_off(t:integer);{выходная процедура счетчика продолжительности выполнения участков кода}
begin
tt[t] := tt[t] + GetTickCount - tt2[t];
end;
и по участкам кода расставив типа:

log_on(N);
log_off(N);

А также общее время выполнения программы через тот же GetTickCount;
Больше всего времени среди всех вычислений тратится на power(x) (заменил на intpower(x) где можно) и на exp (x).

-- 27.02.2014, 20:46 --

Пока нагуглил такое решение - быстрее стандартной Lazarus-ской exp(x), но точность в районе 2-го знака после запятой относительно оригинальной exp(x):

  _fexpA:single=$800000/ln2;  // $800000 = mantissa range
  _fexpF=ln2*ln2/2;
  _fexpB=$400000;
  _fexpE:single=_fexpF/_fexpB;
  _fexpD:single=1-_fexpF*_fexpF;

function  fast_exp(const x:double):double;
  var y:double;
 begin
 {$ASMMODE intel}
 asm
   fld dword ptr [x]
   fld _fexpA
   fmul
   fistp dword ptr [x]
   mov eax,dword ptr [x]
   add eax,$3F800000
   mov dword ptr [x],eax
   fld dword ptr [x]
   and eax, 2*_fexpB-1  // get the remainder
   sub eax, _fexpB      // make middle =0
   mov dword ptr [x],eax
   fild dword ptr [x]    // load remainder
   fld _fexpE
   fmul                 // factor=(r*fexpE)^2 + fexpD
   fmul st(0),st(0)
   fld _fexpD
   fadd
   fmul                // multiply by initial solution
   fstp y
  end;
  fast_exp:=y;
 end;

[venco]

Надеюсь вы используете не Real тип.

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

[arseniiv]

venco в сообщении #831187 писал(а):

Надеюсь вы используете не Real тип.

А во FreePascal он не синоним Double? (Или, если там насчёт этого есть разные режимы компиляции — тогда не по умолчанию ли такой?) В Delphi, помню, Real Turbo Pascal стал Real48, а Real сделали синонимом Double. С седьмой версии точно, хотя хотелось бы верить, что это началось ещё в первой.

[venco]

arseniiv в сообщении #831208 писал(а):

venco в сообщении #831187 писал(а):

Надеюсь вы используете не Real тип.

А во FreePascal он не синоним Double?

Не знаю, я с Паскалем знаком по Турбо. Так там была разница.

[pi-314]

Цитата:

Какую задачу вы решаете? Может вам больше матпакеты вроде Wolfram Mathematica подойдут?

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

ИМХО - язык программирования типа Pascal - гибче, хотя может все дело в привычке )

-- 28.02.2014, 08:50 --

venco в сообщении #831187 писал(а):

Надеюсь вы используете не Real тип.

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

1. А что с Real не так ? (повторюсь - не программист ))
2. Обычный ПК - Intel(R) Core(TM) i5-2500 CPU @ 3.30GHz, Win7 (32 bit)
3. Оказывается в набор инструкций FPU не входит операция вычисления exp(x).

Обнаружил статью по теме (http://www.rsdn.ru/?article/alg/fastpow.xml), изучаю,скорее всего, нужно расширить предложенное там решение на отрицательные значения аргумента, понятно как это сделать на Pasal, но может кто может доработать сразу ASM ? ))

Из статьи следует, что если Free Pascal в той или иной степени базировался на библиотеках Borland - проблема exp (x) - давняя болезнь Pascal/Delphi...

[Sender]

pi-314 в сообщении #831244 писал(а):

А что с Real не так ?

Что-то подсказывает, что лучше использовать Double.

pi-314 в сообщении #831244 писал(а):

скорее всего, нужно расширить предложенное там решение на отрицательные значения аргумента

На первый взгляд, должно работать и для отрицательных.

[rockclimber]

pi-314 в сообщении #831244 писал(а):

1. А что с Real не так ? (повторюсь - не программист ))

Когда-то очень давно, в далекой галактике я читал о том, что с появлением процессоров 80286 тип real стал автоматически преобразовываться процессором в extended, а потом результат вычислений преобразовывался обратно. Из-за этого он работал немного медленнее. Но это было во времена Turbo Pascal 7.0. Сейчас, да еще и во freepascal, легко может оказаться, что тип real является просто алиасом для extended.