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

[David Schultz]

On Wed, Mar 02, 2005, Bruce Evans wrote:
> On Sat, 26 Feb 2005, David Schultz wrote:
> >>: 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?

IIRC, changing the precision is much slower than changing the
rounding mode.  (Intel had to make changing the rounding mode fast
because of C's broken float->int conversion semantics.)
But yes, you're right that we could do this if we decide to change
the rounding precision on x86.

Ideally, libm functions should round the right way for the current
rounding mode, although most of them probably don't do so.  The
polynomials they use internally, as well as the splitting
tricks to get extra precision, probably only work for round to
nearest.  I was careful about setting the rounding mode
appropriately in fma(), and it would not be hard to fix the rest
of libm.  It turns out that once you get round-to-nearest right,
the other rounding modes are not very hard.  This is because in
round-to-nearest, you have to account for the possibility that
the error term might be something like 0.50000000000001 * 2^-p,
whereas the other rounding modes just require one extra bit of
precision.  In most of the libm routines, you'd probably just have
to set the rounding to round-to-nearest, then change it back to
whatever it was before immediately before the final calculation
(often of the form a+b, where b is the tail of a).

Very few people care about evaluating libm functions in those
other rounding modes, so hopefully the performance impact of doing
so would not be very large.

> We support the extension of changing the
> rounding precision.

gcc does constant folding to 53-bit precision on FreeBSD/i386, so
we clearly don't support it very well.

> "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.

Have you tested what the error is with extended precision?  For
some of the transcendental functions, I wouldn't be surprised if
you get an error greater than 1 double precision ulp.