Re: [PATCH] Optimized generic expf and exp2f

[Joseph Myers <joseph at codesourcery dot com>]

On Tue, 5 Sep 2017, Szabolcs Nagy wrote:

> Error checks are inline and actual error handling is
> in separate functions that are tail called and expected
> to be reusable when other math functions are implemented
> without wrappers. (expf, __expf, __ieee754_expf symbols
> are aliases, _LIB_VERSION is not checked, errno is set
> unconditionally according to POSIX rules.)

The _LIB_VERSION and matherr handling is part of the ABI for the existing 
symbol versions of expf and exp2f.  Thus, if you wish to use integrated 
error handling that does not handle _LIB_VERSION and matherr the same way 
as the existing code, you need a new symbol version (with the old version 
then probably keeping the existing wrapper but having it wrap the new 
implementation - duplicate errno setting for existing binaries should be 
fine).

> One issue with the argument reduction is that it only
> works right with nearest rounded rint, otherwise the
> interval is [0,2c] or [-2c,0] instead of [-c,c]. The
> polynomial is optimized for [-c,c] but it has sufficent
> extra precision that it gives acceptable results on
> [-2c,2c] too assuming users are less interested in
> non-nearest rounded precision, however this means some
> glibc ulp error limits will need adjustment.

By ulp error limits do you mean the entries in libm-test-ulps, or the 
global max_valid_error limit in libm-test-support.c which no 
libm-test-ulps entry is allowed to exceed?

> 	* sysdeps/ieee754/flt-32/w_expf_compat.c: Move to...
> 	* math/w_expf_compat.c: ... here.

I'd expect all the w_exp{,f,l}_compat.c files to move at the same time 
(modulo probably needing to add an ldbl-opt version to deal with the long 
double versioning in the ldbl-64-128 and ldbl-128ibm versions).  I think 
they are all in fact generic and are only in sysdeps directories because 
older versions of them used to hardcode bounds on arguments to the exp 
functions that gave finite nonzero results for a particular format.

> 	* sysdeps/x86_64/fpu/w_expf_compat.c: New file.

Is this something to do with x86_64 having its own expf implementations?

It seems i386, ia64, m68k, powerpc64 and x86_64 all have their own 
implementations of expf, exp2f or both (sometimes multiarch, sometimes the 
multiarch variants may have a fallback to using the generic C version, or 
using it built with particular options).  I'd expect an expf replacement 
patch like this one to explain explicitly how it affects those 
architectures.  If on any architecture the answer is that the new C 
versions are not used at all, then I'd expect that architecture to get its 
own dummy versions of math_errf.c and e_exp2f_data.c to avoid building 
unused code into libm.  Of course, if an architecture uses either the 
generic expf or the generic exp2f under any circumstances, it needs the 
relevant support code built into libm when there is code referencing it in 
libm.

Now, I do not make assertions about the performance merits of any of those 
architecture-specific variants; it's entirely possible that some of them 
should be removed if their performance is inferior to the performance of 
this C version (once the C version has been appropriately configured for 
that architecture) (or replaced by building it with specific options to 
generate multiarch variants, if that is beneficial on a particular 
architecture).

> +  if (__builtin_expect (abstop >= top12 (128.0f), 0))

We use __glibc_unlikely in glibc.

> +  if (__builtin_expect (abstop >= top12 (88.0f), 0))

Likewise.

> +#ifndef WANT_ROUNDING
> +/* Correct special case results in non-nearest rounding modes.  */
> +#define WANT_ROUNDING 1
> +#endif

glibc style uses "# define" inside #if, etc.

> +static inline int
> +ieee_2008_issignaling (float x)
> +{
> +  uint32_t ix = asuint (x);
> +  ix ^= 0x00400000; /* IEEE 754-2008 snan bit.  */
> +  return 2 * ix > 2u * 0x7fc00000;
> +}

This doesn't seem to be used, but if you need issignaling tests in future 
functions (powf?), you need to respect HIGH_ORDER_BIT_IS_SET_FOR_SNAN from 
nan-high-order-bit.h.

> +#ifdef __GNUC__
> +#define HIDDEN __attribute__ ((__visibility__ ("hidden")))
> +#define NOINLINE __attribute__ ((noinline))

We don't generally want __GNUC__ conditionals in glibc code, unless it's 
actually shared verbatim with an external repository such as gnulib, or 
it's an installed header.  So attribute_hidden can be used 
unconditionally.

-- 
Joseph S. Myers
joseph@codesourcery.com