correctly rounded mathematical functions

[Paul Zimmermann]

       Dear Joseph,

thank you for your valuable comments.

As a preamble I should say that we will provide "on the shelf" implementations
of correctly rounded functions to all libms (not just glibc), but it will be
the responsability of the libm developers to integrate and maintain them
(however, we will take care of bugs reported against our original version).

> Date: Tue, 4 Jan 2022 19:15:19 +0000
> From: Joseph Myers <joseph@codesourcery.com>
> 
> On Mon, 3 Jan 2022, Paul Zimmermann wrote:
> 
> > We propose to add such correctly rounded functions to the GNU libc,
> > for the three IEEE formats (binary32, binary64, binary128) and the
> > "extended double" format (long double on x86_64).
> 
> Do you have some way of determining the worst cases for correct rounding 
> for binary128 with feasible resource usage, for functions where it's 
> necessary to determine such worst cases by exhaustive search (which is 
> most of them)?

so far, no. But there is some progress in this domain, see for example
https://hal.inria.fr/hal-03240179.

> > cr_exp/cr_log/cr_sin/... functions. Our goal is nevertheless to get the
> > best performance possible, and in some cases our implementation outperforms
> > the GNU libc one.
> 
> If the correctly rounded implementation is faster than the current one, it 
> might make sense to build it twice for glibc, so the non-correctly-rounded 
> function names get the same implementation as the correctly-rounded names 
> but without any additional slower-case checks or computations that are 
> only needed for being correctly rounded and not for results within 1ulp of 
> the correctly rounded value.

yes that would make sense. Note that some of our sample CR functions are faster
than the glibc non-CR ones, for example Alexei Sibidanov measured a latency of
about 75 cycles for cr_tanf, compared to 150 for the glibc 2.33 function on an
i7-8750H.

> (In general, I expect that many of the less accurate libm functions could 
> be reimplemented to get results within 1ulp of the correctly rounded 
> results while becoming faster at the same time, even if correctly rounded 
> versions of those functions are hard.)

our "fast path" codes will deliver results with error < 0.501 ulp (for rounding
to nearest), since in 99.9% of the cases they should give the correct rounding.
This would already be more accurate than all current libraries
(https://members.loria.fr/PZimmermann/papers/accuracy.pdf).

> > We have started to work on two functions (cbrt and acos), for which we
> > provide presumably correctly rounded implementations (up to the knowledge
> > of hard-to-round cases) [2]. [These implementations do replace the current
> 
> For cbrt I suppose you can produce functions that are correctly rounded, 
> without needing exhaustive search, based on elementary results on rational 
> approximation that imply you don't need more than about three times the 
> precision to determine the correctly rounded result - is that what you're 
> doing here?  (Unless there are any effective versions of Thue-Siegel-Roth 
> that give a bound that's actually better for cbrt for the specific 
> precisions relevant here, rather than just asymptotically.)

no, we perform an accurate step with precision about twice the target
precision, and filter out the remaining known hard-to-round cases.

> Note that cr_cbrt is not actually one of the names reserved in the current 
> C23 draft (N2731) (presumably because the cr_* names are for the 
> recommended correctly rounded operations in IEEE 754, and cbrt isn't one 
> of those).  So there may be two separate questions here: do we want cr_* 
> functions reserved in C23, and do we want other cr_* variants of math.h 
> functions for which C23 doesn't reserve the name (such as cr_cbrt or 
> cr_tgamma)?

this is a good question for the CFP group!

> A key question for adding such functions is what impact they have on glibc 
> maintainability.  At present libm functions are generally provided across 
> all supported types - the headers and build system need to know which 
> types are supported, and sometimes which types have the same format as 
> which other types, but they don't generally need to deal with issues such 
> as "function X is supported for long double in some glibc configurations 
> but not others" (and thus to a large extent the headers don't need to know 
> details of e.g. which format long double has; such information is kept as 
> local as possible in a few places in the installed headers).  So it's 
> important to understand how things would be complicated by having such 
> functions.

I see one solution: declare cr_xxx as an alias for xxx for all functions,
and create corresponding bugzilla issues for those which are not (yet) CR.
As new CR functions are integrated within glibc, the corresponding bugs
will be declared as fixed. One could imagine macros like CR_SIN_IS_CR
to tell the user that the cr_sin function is really correctly rounded.

> Also note when adding new functions to glibc that there are a 
> lot of different variations between floating-point support on different 
> architectures to allow for, including:
> 
> * before-rounding and after-rounding tininess detection;

we probably will have to take care of this, but I guess there will be very few
cases per function, and they can be dealt in the "exceptional" cases if needed.

> * architectures may or may not have hardware fma support;

we will take care our code runs with or without fma.

> * architectures might have just round-to-nearest, or just a subset of C99 
> rounding modes, or all the C99 rounding modes;

our functions will be CR for the four IEEE rounding modes, thus if there is
only round-to-nearest, this should not be a problem (we won't use fesetround
nor fegetround)

> * architectures might or might not support floating-point exceptions.

it is not yet clear to me, but we might leave that to the libm developers,
since some libms might not care about exceptions

> So some important questions to consider for understanding how glibc might 
> be complicated by adding these functions include:
> 
> * How will you arrange for installed headers to declare the supported 
> functions for those types with formats for which the correctly rounded 
> functions are supported but not for types for which they are not supported 
> - what will the conditionals in bits/mathcalls.h (for example) look like?

see my proposal above

> * Likewise, what will things look like in the makefiles to arrange for the 
> functions to be built and tested when supported but not otherwise?

idem

> * When you have correctly rounded functions for x86 "extended" (ldbl-96 in 
> glibc), do those also work properly for the m68k variant of that format 
> (which handles biased exponent 0 like other biased exponents, so has 
> normal and subnormal exponent ranges going one exponent smaller), or do 
> you avoid providing the cr_* functions for that format variant?

at the moment, unless we get more manpower, we don't plan to provide the
m68k variants.

> * Are there any dependencies on features of the architecture (for example, 
> do any functions depend on fma support, or on the architecture supporting 
> exceptions or all the IEEE rounding modes)?

some of our sample implementations depend on fma support or specific
intrinsics, but we plan to provide a "generic" implementation for all
functions.

> * How do you arrange for test inputs in auto-libm-test-in to be shared for 
> the cr_* functions (while presumably needing separate auto-libm-test-out-* 
> and lib-test-*.inc files for them, since gen-auto-libm-tests does some 
> things differently for functions required to be correctly rounded)?

this question is more for the glibc developers (if any) who are interested
to integrate CR functions, I let them answer.

> * What configurations are you doing execution testing on?  (Build testing 
> with build-many-glibcs.py is also generally a good idea for new libm 
> functions to e.g. make sure you got all the ABI baseline updates right, 
> given all the variations in supported types between different 
> architectures, but execution testing should also be done for enough 
> architectures to cover all the main variations in what code is used on 
> what architectures.)

we are at the very beginning of the project, and it will take years to
complete, but we plan to test our functions on several configurations,
and surely we'll also get feedback from the libm developers.

> * Where applicable, do you handle both before-rounding and after-rounding 
> tininess detection (that matters for e.g. cr_sin (DBL_MIN), which is tiny 
> with before-rounding tininess detection but not with after-rounding 
> tininess detection)?  The semantics of cr_* include getting this right - 
> raising the correct exceptions, including this case where it's 
> architecture-specific what exceptions should be raised.  (If explicit 
> conditionals in the code are needed, use glibc's tininess.h 
> TININESS_AFTER_ROUNDING macro; see sysdeps/ieee754/dbl-64/s_fma.c for 
> example.)

as said above, I guess we can deal with these cases one by one
(so far we don't take care of underflow exceptions)

> * Does the analysis that shows functions to be correctly rounded cover 
> FE_TONEARESTFROMZERO as well as the four C99 rounding modes (i.e., does it 
> cover all the current IEEE rounding modes)?  We don't currently have any 
> support for FE_TONEARESTFROMZERO, but the RISC-V architecture supports it, 
> and we don't want these functions to be an obstacle to adding support for 
> FE_TONEARESTFROMZERO on RISC-V in future.

so far we did check only for the four C99 rounding modes, but once
FE_TONEARESTFROMZERO is supported say in glibc, it should be easy to
run the hard-to-round cases with that rounding mode.

Best regards,
Paul