cvs commit: src/lib/msun/src e_expf.c

[Bruce Evans]

On Sat, 26 Feb 2005, David Schultz wrote:

> On Sat, Feb 26, 2005, M. Warner Losh wrote:
>> In message: <20050226023149.GA63314 at VARK.MIT.EDU>
>>             David Schultz <das at FreeBSD.ORG> writes:
>> : On Thu, Feb 24, 2005, David Schultz wrote:
>> : >   ...
>> : >   Log:
>> : >   Revert rev 1.8, which causes small (e.g. 2 ulp) errors for some
>> : >   inputs.  The trouble with replacing two floats with a double is that
>> : >   the latter has 6 extra bits of precision, which actually hurts
>> : >   accuracy in many cases.  All of the constants are optimal when float
>> : >   arithmetic is used, and would need to be recomputed to do this right.
>> :
>> : This is related to a good reason why we can't switch the default
>> : precision on i386 to extended.  Many of the functions in libm use
>> : minimax approximations, which are ``optimal'' approximations in
>> : the sense that their maximum error over all in-range inputs is the
>> : smallest possible (unless more terms are used).  These approximations
>> : take rounding error into account, so when the machine precision is
>> : increased, they're no longer optimal and the error in the approximation
>> : can increase significantly.  ...

I didn't know that there was such a good reason.

The library could switch to the appropriate mode.  That wuld be
inefficient, but it may need to be done in some cases anyway when the
user has changed the FP environment.  I'm not familiar with all the
details of what C99 allows.  It allows or requires changes to the
rounding mode to affect most functions, but shouldn't this only affect
rounding of the result and not break algorithms that depending on
certain rounding internally?  We support the extension of changing the
rounding precision.  This has a much larger effect on the algorithms.
i387/e_exp.S already changes the rounding precision to 64 bits, since
it needs that and the FreeBSD default is 53 bits.  I only recently
remembered changing it do do this, and found a 10 year old patch set
with macros to switch the (i387) mode for all of fdlibm.

>> I guess it comes down to which is more important: long double working
>> outside of the double range, or small errors in these functions.  I
>> think the former is, but will defer to those with greater floating
>> point foo.  Since to make long double working by setting extended mode
>> in the application would also break things in the manner you describe.

"small" (1-ulp) errors in these functions are 2^11 times larger than
1-ulp errors in extended precision, so you get few benefits from
extended precision if you call one of these functions and it is
broken.  Even if it isn't broken, then it will have 0.5-ulp errors.
Mainly code that never calls a math function or calls only long double
math functions that actually exist and actually are accurate to nearly
1 extended precision ulp can benefit.

> ...
> On the bright side, very few people have noticed the problem in
> the last decade, and once amd64 takes over the world, nobody will
> care anymore at all.  ;-)

And they will mostly use doubles and not long doubles, and thus get
only pure double precision.

Bruce