correctly rounded mathematical functions

[Paul Zimmermann <Paul.Zimmermann@inria.fr>]

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/