p.lang.pascal.borland:412
From: WE@completely.invalid (Wolfgang Ehrhardt)

On Sun, 16 Jan 2011 20:13:26 +0000, Robert AH Prins
 wrote:
>Two more replies from Norbert Juffa:
...
>I got repro. The failing input is 4003_A65AF678_54B28211. With rounding 
>to positive infinity, multiplying the input by L2E gives 30.0 + 1 ulp. 
>FRNDINT then rounds this up further to 31.0. Things are already 
>irreparably screwed up at this point. In any event, the input to F2XM1 
>becomes -(1.0 + 1 ulp), which is clearly out of bounds. On my CPU it 
>returns the unmodified input operand in this case.

Meanwhile I found that the unmodified code has a lot more errors for
rmTruncate with negative x arguments: here are the smallest values
found (also with hex dumps for better reference), the positive
arguments give wrong results with rmUp the negative with rmTruncate.

        x               exp(x)                 x as LSB Hex
-6.238324625040 -1.27054942088145051E-0021  C001C7A05AF6CC0968E2
-7.624618986159 -3.17637355220362627E-0022  C001F3FCE0F4C07D474D
-8.317766166719 -5.29395592033937712E-0023  C002851591F9DD5B9B41
-9.010913347279 -7.94093388050906568E-0023  C002902CB3795A7892DC

 6.238324625040 -1.11022302462515654E-0016  4001C7A05AF6CC0968E1
 7.624618986159 -4.44089209850062616E-0016  4001F3FCE0F4C07D474C
12.476649250079 -1.13686837721616030E-0013  4002C7A05AF6CC0968E1
15.249237972319 -1.81898940354585648E-0012  4002F3FCE0F4C07D474C

>The only reasonable workaround for this that I can see is to reduce the 
>input to F2XM1 further, to make sure we are definitely inside the 
>interval supported by the hardware. In other words, replace F2MX1 by the 
>following instruction sequence:
>
>  fmul  dword ptr [half]   ; 0.5*frac(z) int(z)
>  f2mx1                    ; 2^(0.5*frac(z))-1 int(z)
>  fld   st(0)              ; 2^(0.5*frac(z))-1 2^(0.5*frac(z))-1 int(z)
>  fadd  dword ptr [two]    ; 2^(0.5*frac(z))+1 2^(0.5*frac(z))-1 int(z)
>  fmulp st(1),st           ; 2^frac(z)-1 int(z)
>
This code does not show the errors and seems to work ok (Note that
there is a small type f2mx1 must be changed to f2xm1; and the
constants should be single). The priced to pay is a small decrease in
accuracy, here the result for 20000 random arguments with standard
round to nearest (the error is calculated with my MPArith routines):

 ------ Old version ------
Test of amath.exp
at 10000 random values in [-100.0000000000 .. 100.0000000000]
RMS = 0.24, max rel = 0.71 eps at 
 x(ext) =  7.24700653776128467E+0001 = $400590F0AC68BFA87B80
 y(ext) =  2.97405843003663398E+0031 = $4067BBB0815EA0230490
 y(mpf) = +2.97405843003663397839071382253654488039888365036E+31

Test of amath.exp
at 10000 random values in [100.0000000000 .. 10000.0000000000]
RMS = 0.24, max rel = 0.72 eps at 
 x(ext) =  4.80456952900561126E+0003 = $400B96248E653929CA47
 y(ext) =  3.96309394729273526E+2086 = $5B12B8A5DD42BAF8E736
 y(mpf) = +3.96309394729273526538566268155478642708748950161E+2086

 ------ New version ------
Test of amath.exp
at 10000 random values in [-100.0000000000 .. 100.0000000000]
RMS = 0.25, max rel = 0.83 eps at
 x(ext) =  1.83770128884148676E+0001 = $400393041F554F4C8000
 y(ext) =  9.57271857222288938E+0007 = $4019B695CA371C7FC4DE
 y(mpf) = +9.57271857222288937714622943945524368866640537394E+7

Test of amath.exp
at 10000 random values in [100.0000000000 .. 10000.0000000000]
RMS = 0.24, max rel = 0.83 eps at
 x(ext) =  5.98844546289204336E+0003 = $400BBB23904ED947432E
 y(ext) =  5.60815118133691577E+2600 = $61BEB51753447EF41B4A
 y(mpf) = +5.60815118133691576592733747663153598713873826468E+2600

>If there is consensus that Pascal math functions shall work properly 
>when the dynamic rounding mode is not round-to-nearest-even, I believe 
>this is the way to go (least impact on execution time).
>
>I dug out my old Pascal manuals and can't find anything about the 
>prescribed behavior of math functions when a user dials in some rounding 
>mode different from the default. Given the naive approach to exp() 
>computation in Borland's original library, it is quite possible that it 
>happened to work as desired in round-up mode. My question is: is that to 
>be considered an implementation artifact or a feature?

I guess it is a pure accident, especially as Borland Pascal has no
library functions for getting/setting rounding modes. The situation is
a bit different for Virtual Pascal and Delphi. But IMO at least the
Delphi math unit designer/implementers did not understand there task.
My AMath unit is an attempt to overcome these Delphi shortcomings. 

Robert: Please forward a big "Thank you" to Norbert.

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)