correctly rounded mathematical functions

correctly rounded mathematical functions

[Paul Zimmermann]

The current C working draft [1, p392] has reserved names for correctly
rounded functions (cr_exp, cr_log, cr_sin, ...).

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

These implementations should be correctly rounded for all rounding modes,
for example one could do the following to emulate interval arithmetic:

   fesetround (FE_DOWNWARD);
   y_lo = cr_exp (x_lo);
   fesetround (FE_UPWARD);
   y_hi = cr_exp (x_hi);

These functions will not replace the current functions (exp, log, sin, ...).
Users who want a fast implementation will call the exp/log/sin/... functions,
users who want a correctly rounded function and thus reproducible results
(whatever the hardware, compiler or operating system) will use the
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.

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
glibc ones, since it was simpler to demonstrate our skills.]

Christoph Lauter
Jean-Michel Muller
Alexei Sibidanov
Paul Zimmermann

[1] http://www.open-std.org/jtc1/sc22/wg14/www/docs/n2596.pdf
[2] https://homepages.loria.fr/PZimmermann/CORE-MATH/

Re: correctly rounded mathematical functions

[Joseph Myers]

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)?

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

(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.)

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

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)?


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

* architectures may or may not have hardware fma support;

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

* architectures might or might not support floating-point 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?

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

* 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?

* 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)?

* 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)?

* 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.)

* 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.)

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

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: 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

Re: correctly rounded mathematical functions

[Joseph Myers]

On Wed, 5 Jan 2022, Paul Zimmermann wrote:

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

That is contrary to recent practice for e.g. new Linux kernel syscalls, 
where we no longer tend to add emulations if they have problems such as 
race conditions meaning they are not a good emulation of the semantics of 
the syscall.  Presumably if someone uses cr_* it's because they need 
correct rounding, just as if they use pselect it's because they need to 
avoid a race (cf. bug 9813 criticising the emulation added before we 
tended to avoid such emulations).

Now, in the syscall case, the functions produce an ENOSYS error when the 
syscall is unavailable.  But that's not such a good idea for libm 
functions, given that users don't expect such functions to fail (other 
than domain/range/pole errors) and so won't check for such an error - and 
even if you use the stubs mechanism, that only produces a link-time 
warning, and while autoconf knows about not detecting as a available a 
function glibc defines only as a stub, other configuration systems might 
not.  So I think for these functions, not defining or declaring them at 
all when we don't have an implementation with the right semantics is the 
better approach - but it also introduces complexity in the installed 
headers.


One more consideration to mention: the cr_* names reserved in draft C23 
include some for functions that are new in C23 (from TS 18661-4) and not 
currently (except for the exp10 functions) supported by glibc.  It would 
be reasonable to say that we don't add cr_<func> before we have the 
underlying, not-correctly-rounded, function <func> (for all floating-point 
types and formats for all architectures) - so no cr_exp2m1 before we have 
exp2m1, for example.

I don't expect that to be a major issue, given that (a) most of those new 
functions can be implemented reasonably, if not as fast or precise as 
ideal, with fairly straightforward generic implementations in terms of 
other functions, (b) I expect the C99 functions are a priority for your 
correctly-rounded functions work over the newer ones and (c) when I get 
time I hope to implement those new functions for glibc if no-one else has 
done them by then, just as with other new C23 features.

Those new functions do have a soft dependency on MPFR support in order to 
generate expected test results (it's not impossible to add them without 
the MPFR support, but it means adding a local emulation in 
gen-auto-libm-tests), but most of that MPFR support is in current git MPFR 
(to be 4.2) - the main exception being that rootn with negative integer 
argument has no corresponding MPFR function.  (If gen-auto-libm-tests is 
to be usable on 32-bit hosts, which probably isn't that important, there 
would also be the issue that the MPFR functions for rootn and compoundn 
use unsigned long / long as their integer arguments but the C23 functions 
can have any integer argument in the range of signed long long int.)

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: correctly rounded mathematical functions

[Paul Zimmermann]

       Dear Joseph,

> Date: Wed, 5 Jan 2022 18:44:46 +0000
> From: Joseph Myers <joseph@codesourcery.com>
> 
> On Wed, 5 Jan 2022, Paul Zimmermann wrote:
> 
> > 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.
> 
> That is contrary to recent practice for e.g. new Linux kernel syscalls, 
> where we no longer tend to add emulations if they have problems such as 
> race conditions meaning they are not a good emulation of the semantics of 
> the syscall.  Presumably if someone uses cr_* it's because they need 
> correct rounding, just as if they use pselect it's because they need to 
> avoid a race (cf. bug 9813 criticising the emulation added before we 
> tended to avoid such emulations).

this was just one proposal, feel free to propose a better one!

> Now, in the syscall case, the functions produce an ENOSYS error when the 
> syscall is unavailable.  But that's not such a good idea for libm 
> functions, given that users don't expect such functions to fail (other 
> than domain/range/pole errors) and so won't check for such an error - and 
> even if you use the stubs mechanism, that only produces a link-time 
> warning, and while autoconf knows about not detecting as a available a 
> function glibc defines only as a stub, other configuration systems might 
> not.  So I think for these functions, not defining or declaring them at 
> all when we don't have an implementation with the right semantics is the 
> better approach - but it also introduces complexity in the installed 
> headers.
> 
> 
> One more consideration to mention: the cr_* names reserved in draft C23 
> include some for functions that are new in C23 (from TS 18661-4) and not 
> currently (except for the exp10 functions) supported by glibc.  It would 
> be reasonable to say that we don't add cr_<func> before we have the 
> underlying, not-correctly-rounded, function <func> (for all floating-point 
> types and formats for all architectures) - so no cr_exp2m1 before we have 
> exp2m1, for example.

yes that make sense.

> I don't expect that to be a major issue, given that (a) most of those new 
> functions can be implemented reasonably, if not as fast or precise as 
> ideal, with fairly straightforward generic implementations in terms of 
> other functions, (b) I expect the C99 functions are a priority for your 
> correctly-rounded functions work over the newer ones and (c) when I get 
> time I hope to implement those new functions for glibc if no-one else has 
> done them by then, just as with other new C23 features.

yes our priority will be the C99 functions.

> Those new functions do have a soft dependency on MPFR support in order to 
> generate expected test results (it's not impossible to add them without 
> the MPFR support, but it means adding a local emulation in 
> gen-auto-libm-tests), but most of that MPFR support is in current git MPFR 
> (to be 4.2) - the main exception being that rootn with negative integer 
> argument has no corresponding MPFR function.  (If gen-auto-libm-tests is 
> to be usable on 32-bit hosts, which probably isn't that important, there 
> would also be the issue that the MPFR functions for rootn and compoundn 
> use unsigned long / long as their integer arguments but the C23 functions 
> can have any integer argument in the range of signed long long int.)

mpfr_rootn_si should not be difficult to add. I'll see with the MPFR
developers :-)

Best regards,
Paul

Re: correctly rounded mathematical functions

[Carlos O'Donell]

On 1/3/22 07:48, Paul Zimmermann wrote:
> The current C working draft [1, p392] has reserved names for correctly
> rounded functions (cr_exp, cr_log, cr_sin, ...).
> 
> 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).
> 
> These implementations should be correctly rounded for all rounding modes,
> for example one could do the following to emulate interval arithmetic:
> 
>    fesetround (FE_DOWNWARD);
>    y_lo = cr_exp (x_lo);
>    fesetround (FE_UPWARD);
>    y_hi = cr_exp (x_hi);
> 
> These functions will not replace the current functions (exp, log, sin, ...).
> Users who want a fast implementation will call the exp/log/sin/... functions,
> users who want a correctly rounded function and thus reproducible results
> (whatever the hardware, compiler or operating system) will use the
> 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.
> 
> 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
> glibc ones, since it was simpler to demonstrate our skills.]
> 
> Christoph Lauter
> Jean-Michel Muller
> Alexei Sibidanov
> Paul Zimmermann

In general I think correctly rounded functions are something we would want to support.

For a very long time we've considered that there are really 3 sets of users:

(a) Loosely correct answers provided quickly (high ULPs errors allowed, >9 ULPs)
(b) Correct answers (no errors >9 ULPs)
(c) Correctly rounded answers (0 ULPs)

Note: Rounding modes matter.

To date we have been removing the mutli-precision fallback paths and driving glibc's
implementation towards (b) where possible (with some IBM double-double exceptions).

I would like to have (c) included in glibc for selection at link-time depending on
the needs of the users, but as Joseph points out, the costs to maintenance are real
and need to be discussed.
 
> [1] http://www.open-std.org/jtc1/sc22/wg14/www/docs/n2596.pdf
> [2] https://homepages.loria.fr/PZimmermann/CORE-MATH/
> 


-- 
Cheers,
Carlos.

Re: correctly rounded mathematical functions

[Carlos O'Donell]

On 1/5/22 11:03, Paul Zimmermann wrote:
> 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.

This isn't a good idea given the expectations that users will have for calling the
cr_* functions. Either they are defined and have the expected semantics or they are
not defined. However, this makes the semantics of the ABI a bit difficult, we would
probably want the ABI to appear as implementations are made available, and that can
be on a per-arch basis.

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

In these cases then we would not define the functions, and they would not be available
for use. This poses a porting hazard for the cr_* functions, which is unfortunate, but
not the worst. I'm not opposed to per-arch variations in ABI.

-- 
Cheers,
Carlos.

Re: correctly rounded mathematical functions

[Paul Zimmermann]

       Dear Carlos,

> For a very long time we've considered that there are really 3 sets of users:
> 
> (a) Loosely correct answers provided quickly (high ULPs errors allowed, >9 ULPs)
> (b) Correct answers (no errors >9 ULPs)
> (c) Correctly rounded answers (0 ULPs)
> 
> Note: Rounding modes matter.

I would say that rounding modes only matter for (c), since for (a) and (b)
most users will use the default rounding-to-nearest-ties-to-even mode, and
if not, another rounding mode will maybe add a few ulps to the error, which
should not matter for them.

However for (c), users who want to perform interval arithmetic for example
will want correct rounding for all rounding modes.

Progress of the CORE-MATH project can be followed on
https://homepages.loria.fr/PZimmermann/CORE-MATH/.
Take a look at the tanf code and the associated performance graph:
our correctly rounded code is up to 3 times faster than GNU libc.
(and also outperforms the Intel Math Library).

We plan to release our code under the MIT license, I believe this is
compatible for inclusion into GNU libc.

Best regards,
Paul