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)