comp.lang.pascal.misc:162
From: WE@completely.invalid (Wolfgang Ehrhardt)

Virtual Pascal and Borland Pascal with Robert Prins' bpl70v20 have an
enhanced exp routine based on a code snippet from Norbert Juffa. This
is more accurate than Borland's standard Pascal/Delphi exp or older VP
routines.

The key idea is to calculate exp(x) = 2^(log2(e)*x - 1) + 1 with extra
precision for frac(log2(e)*x) and then use the x87 instruction f2xm1
to compute 2^frac(log2(e)*x)-1.  f2xm1 requires that the abs value of
its argument must be <= 1.

Unfortunately the Juffa's code miserably fails for some x values if
the "round up toward +infinity" is used. Example:

program t_expbug;
{$N+}
const
  rmUp: word = $1B72;
var
  x,y: extended;
begin
  asm
    fldcw [rmUp]
  end;
  x := 20.0 + 0.7944154167983592825;
  y := exp(x);
  writeln(x:18:14, y:27, 1073741824.0:27);
  x := 2*x;
  y := exp(x);
  writeln(x:18:14, y:27, 1152921504606846976.0:27);
end.

The output is for both Virtual Pascal 2.1.279 and BP7/bpl70v20:

20.79441541679836 -2.32830643653869629E-10  1.07374182400000000E+09
41.58883083359672 -5.00000000000000000E-01  1.15292150460684698E+18

A fix is to test frac() against 1 and correct. I guess this is
faster than temporarily change the rounding control to "round to
nearest". 

Here is a corrected version (suppressing most comments in order to
display properly in newsreaders)

const
 ln2_hi: array[0..7] of byte = ($00,$00,$E0,$FE,$42,$2E,$E6,$3F);
 ln2_lo: array[0..7] of byte = ($76,$3C,$79,$35,$EF,$39,$EA,$3D);

function exp(x: extended): extended;assembler;{&Frame-}{&Uses none}
  {-Accurate and debugged exp}
asm
  {Based on Norbert Juffa's exp}
  fld     [x]
  fldl2e
  fmul    st,st(1)
  frndint
  fld     qword ptr [ln2_hi]
  fmul    st,st(1)
  fsubp   st(2),st
  fld     qword ptr [ln2_lo]
  fmul    st, st(1)
  fsubp   st(2),st
  fxch    st(1)
  fldl2e
  fmulp   st(1),st

  {It may happen (especially for rounding modes other than
   "round to nearest")  that |frac(z)| > 1. In this case the
   result of f2xm1 is undefined. The next lines will test
   frac(z) and truncate it to +1 or -1 if necessary.}

  fld     st
  fabs
  fld1
  fcompp
  fstsw   ax
  sahf
  jae     @@1
  fldz
  fcompp
  sahf
  fld1
  jae     @@1
  fchs
@@1:

  f2xm1
  fld1
  faddp   st(1),st
  fscale
  fstp    st(1)
  fwait
end;


Any comments or improvements?
Wolfgang

-- 
Do not use the e-mail address from the header! You can get a
valid address from http://home.netsurf.de/wolfgang.ehrhardt
(Free open source Crypto, CRC/Hash, MPArith for Pascal/Delphi)

--- Internet Rex 2.31
 * Origin: The gateway at Omicron Theta (1:261/20.999)