Text 434, 127 rader
Skriven 2011-01-16 20:13:26 av Robert AH Prins
Ärende: Re: exp bug in Virtual Pascal and bpl70v20
==================================================
Gecko/20101207 Thunderbird/3.1.7
comp.lang.pascal.borland:411
From: Robert AH Prins <spamtrap@prino.org>

On 2011-01-16 14:45, Wolfgang Ehrhardt wrote:
>>
>> 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.
>
> Here is a quote from the Virtual Pascal bug list found at
>
<https://admin.topica.com/lists/virtualpascal@topica.com/read/message.html?sort=a&mid=908867704&start=7>

> [quote]...
> In the old days, F2XM1 was limited to arguments in
> [-0.5, 0.5]. But since the 387, F2XM1 can accecpt
> arguments in [-1, 1].
>
> So, we can split the argument into an integer and a
> fraction part using FRNDINT and the fraction part
> will always be -1<= f<= 1 no matter what rounding
> control. This means we don't have to load/restore the
> FPU control word (CW) which is slow on modern
> OOO FPUs (since FLDCW is a serializing instruction).
> ....
> -- Norbert
> [/quote]
>
> So at least all rounding modes are supposed to work. But I see no
> way to PROVE that even for "round to nearest" the calculated
> fraction part will ALWAYS be in [-1, 1].
>
> Up do now, I have only found errors (more correct: negative exp
> values) for the rmUp rouding mode for arguments near x = k*ln(2)*2^e
> with (k=9,11,18,22,30,33,36,41,44,47,60,63) and values -2<= e<= 10.
>
>>
>> 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 think, the RTLs of Borland (Pascal and Delphi) or Virtual Pascal are
> very badly designed for the math functions. Most (all?) of them return
> extended result but very few (with the possible exceptions sqrt and
> ln, because these can be directly mapped to the FPU instructions)
> return correctly rounded results. A better designed RTL would return
> correctly rounded double results (round to nearest).
>
> Other rounding modes could have relaxed accuracy requirements but
> they should to return totally wrong results.

Two more replies from Norbert Juffa:

=== 1 ===
So far I am unable to reproduce this problem. I don't have a working 
Pascal installation anymore, so I tried with inline assembler inside 
Microsoft C. Maybe I don't have quite the right input operand. What is 
the byte pattern (10 bytes) for x (i.e. 20.0 + 0.7944154167983592825) 
that Pascal produces ?

 From code inspection, it's not quite clear to me how we could wind up 
with a result that is too large (i.e. > 1.0 prior to F2XM1). If the 
rounding mode is "up", each of the FMUL products with ln2 should be 
slightly too large, and subtracting that from the input should result in 
a fraction that is too small, not too large. At least if all operands 
are positive. Having a repro case might at least lead to better 
understanding as to what is going on.

I doubt that there is an easy way to make the code functional regardless 
of rounding mode; certainly analysis would be very tricky as all 
intermediate operations are subject to the direct rounding mode. As I 
said the code was not designed with directed rounding in mind, as I 
recall I designed the code for use with Digital Fortran originally.

=== 2 ===
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.

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)

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 the above is the best than can be done, if you really need 
non-default rounding.

Thanks Norbert!

Robert

-- 
Robert AH Prins
spamtrap(a)prino(d)org

--- Internet Rex 2.31
 * Origin: The gateway at Omicron Theta (1:261/20.999)