XPost: comp.lang.pascal.misc   
   From: spamtrap@prino.org   
      
   On 2011-01-16 00:37, Robert AH Prins wrote:   
   > 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.   
      
   His initial reply:   
      
   Robert:   
      
   Math library functions targetting x87 were typically designed with the   
   assumption that the dynamic rounding control of the x87 FPU at entry is   
   set to round-to-nearest-or-even. That is certainly the specification to   
   which I designed this code, so the results observed are a consequence of   
   operating my code outside its design specifications.   
      
   Does documentation for Borland Pascal (or Virtual Pascal) state anywhere   
   that common math library routines are guaranteed to work correctly if a   
   user changes the x87 rounding mode away from the default setting, prior   
   to invoking such a function? If so, what is the stated definition of the   
   correct behavior in such a case?   
   For example, the specification may state that the result returned must   
   be the same for all rounding modes, which may require each function to   
   use a wrapper that saves the current rounding mode, sets it to round to   
   nearest, then restores it at the end. Obviously this would be slow.   
      
   I will take a quick look to see if there is a small modification to the   
   code that gets results closer to what the users expects, but no promises.   
      
   I'll get back to you by Sunday night.   
      
   -- Norbert   
      
      
   --   
   Robert AH Prins   
   spamtrap(a)prino(d)org   
      
   --- SoupGate-Win32 v1.05   
    * Origin: you cannot sedate... all the things you hate (1:229/2)