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)