RFC: removing slow paths in various dbl-64 libm functions

[Joseph Myers <joseph at codesourcery dot com>]

Various dbl-64 libm functions (pow, exp, log, sin, cos, tan, asin, acos, 
atan, atan2) contain slow paths that try to achieve correctly rounded 
results by using extra precision internally - in some cases, a lot of 
extra precision, resulting in the functions being very slow (bugs 5781, 
13932, 17211).  The correct rounding is in general for round-to-nearest 
only (various of those functions set round-to-nearest internally because 
they do things that do not work in other rounding modes).

We now have documented libm accuracy goals which do not involve these 
functions being correctly rounded (and indeed on some architectures other 
implementations are used that are not correctly rounded).  We also have TS 
18661-4 defining names such as crsin that are reserved for 
correctly-rounded (for all rounding modes) versions of those functions.

I'd like to propose that we consider it OK to remove those slow paths, and 
thereby make the functions not correctly rounded, if this would not result 
in errors from those functions of more than 1ulp in round-to-nearest or 
2ulp in other rounding modes (this is stronger than the limits 
libm-test.inc places on libm function errors, but might be a longer-term 
goal for general libm function accuracy).  The error analysis may be done 
on the basis of the existing conditions in those functions for when to use 
the slow paths.  For example, in e_asin.c:

    res = x+t;         /*  res=arcsin(x) according to Taylor series  */
    cor = (x-res)+t;
    if (res == res+1.025*cor) return res;

This is a typical example of code checking whether to use a slow path.  If 
it could not be assumed that the computed result is within 1ulp of the 
exact value, then the res == res+1.025*cor condition would not be valid 
for determining whether the result is necessarily correctly rounded; the 
implied error analysis is that x+t is within 0.5*(1 - 1/1.025)ulp of the 
exact value.  So, given my proposal, this code could just do "return 
res;", with both following slow cases removed.

The slow paths in question include uses of __ieee754_expl in exp and pow 
if USE_LONG_DOUBLE_FOR_MP (defined only in 
sysdeps/powerpc/power4/fpu/Makefile).

This does not affect cases where multiple-precision slow paths are needed 
simply to get results within the normal accuracy goals for libm functions, 
or where they are involved in functions that are fully defined by 
reference to IEEE 754 operations and so *are* required to be correctly 
rounded.

-- 
Joseph S. Myers
joseph@codesourcery.com