Gecko/20101207 Thunderbird/3.1.7   
   comp.lang.pascal.misc:163   
   From: Robert AH Prins    
      
   On 2011-01-15 16:54, Wolfgang Ehrhardt wrote:   
   > 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?   
      
   I've forwarded your code to the grandmaster of everything FPU, let's    
   wait and see.   
      
   Robert   
   --    
   Robert AH Prins   
   spamtrap(a)prino(d)org   
      
   --- Internet Rex 2.31   
    * Origin: The gateway at Omicron Theta (1:261/20.999)