From: Robert Lee Subject: Fast math Date: 14 Dec 1999 00:00:00 GMT Message-ID: <38567DAD.8A1743E5@nwu.edu> Content-Transfer-Encoding: 7bit Organization: Another Netscape Collabra Server User X-Accept-Language: en Content-Type: text/plain; charset=us-ascii Mime-Version: 1.0 Newsgroups: borland.public.delphi.basm I was going to wait a just put this stuff on my web site, but it is so painfully quite here that I post it and an attempt to generate something that at least resembles assembler discussion: I ran across a fast exponentation routine (Neural Computation, May 1999) that uses the storage format of floating numbers to generate a quick but low precision exp function. Not new in and of itself, but they figured out a way to improve the accuracy by letting the integer math for the exponent spill over into the mantissa and thus get a form of interpolation. Pretty nifty. I've converted it to basm (quite a bit faster than the pascal on a PII) and also came up with some ancillary routines for a sigmoid, division, and natural log. Here they are, comments and tweaks encouraged: // dxxx is a double based version // fxxx is a single based version const ln2=0.693147181; dexpA:double=\$100000/ln2; // 2^20 / ln2 dlnA:double=ln2/\$100000; fexpA:double=\$800000/ln2; // 2^23 / ln2 flnA:double=ln2/\$800000; fexpC=545947; flnC=545947; dexpC=68243; fnexpA:double=-\$800000/ln2; // 2^23 / ln2 fone:single=1; fdivA:single=-1; function dexp(x:single):double; var y:double; asm fld x fmul dexpA fistp dword ptr y[4] mov dword ptr y[0],0 add dword ptr y[4],(\$3FF00000-dExpC) fld y end; function dln(x:double):double; asm sub dword ptr x[4],(\$3FF00000-dExpC); fild dword ptr x[4] fmul dlnA end; function fexp(x:single):double; // ~ 10x faster asm fld x fmul fexpA fistp dword ptr x add dword ptr x,(\$3F800000-fExpC) fld x end; function fexp1(var x:single):double; // special case when x can be prevented from being pushed on the stack // ~14.5x faster asm fld dword ptr [eax] fmul fexpA fistp dword ptr [esp-4] add dword ptr [esp-4],(\$3F800000-fExpC) fld dword ptr [esp-4] end; function fln(x:single):double; // ~ 7x faster // very poor accuracy as x -> 1.0 asm sub dword ptr x, (\$3F800000-flnC) fild dword ptr x fmul flnA end; function fdiv(x:single):single; // ~20 faster in some circumstances var y:single; asm sub dword ptr x, (\$3F800000-flnC) fild dword ptr x fmul fdivA fistp dword ptr y add dword ptr y,(\$3F800000-fExpC) fld y // 1/x end; function sigmoid(x:single):single; // ~1.6x faster than 1/(1+fexp) asm fld x fmul fnexpA fistp dword ptr x add dword ptr x,(\$3F800000-fExpC) fld x fld dword ptr fone fadd st(1),st(0); fdivr end; function sigmoid1(var x:single):single; // ~2.3x faster than 1/(1+fexp) asm fld dword ptr [esp-4] fmul fnexpA fistp dword ptr [esp-4] add dword ptr [esp-4],(\$3F800000-fExpC) fld dword ptr [esp-4] fld dword ptr fone fadd st(1),st(0); fdivr end; function sigmoid0(x:single):single; begin result:=1/(1+fexp(-x)); end; -- Bob Lee High Performance Delphi - http://www.econos.com/optimize/