mailing list of musl libc
 help / color / mirror / code / Atom feed
* [musl] correctly rounded mathematical functions
@ 2022-01-03 13:07 Paul Zimmermann
  2022-01-04  0:04 ` Damian McGuckin
  2022-01-04  2:29 ` Rich Felker
  0 siblings, 2 replies; 6+ messages in thread
From: Paul Zimmermann @ 2022-01-03 13:07 UTC (permalink / raw)
  To: musl; +Cc: sibid, christoph.lauter, Jean-Michel.Muller

      Dear Musl developers,

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

We propose to provide such correctly rounded implementations
for the three IEEE formats (binary32, binary64, binary128) and the
"extended double" format (long double on x86_64).

These implementations will 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);

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.

Our objective is to provide open-source implementations that can be integrated
in the major mathematical libraries (GNU libc, Intel Math Library, AMD Libm,
Redhat Newlib, OpenLibm, Musl, llvm-libc, CUDA, ROCm).

Are developers of Musl interested by such functions?
If so, we could discuss what would be the requirements for integration in
Musl in terms of license, table size, allowed operations.

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

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/

^ permalink raw reply	[flat|nested] 6+ messages in thread

* Re: [musl] correctly rounded mathematical functions
  2022-01-03 13:07 [musl] correctly rounded mathematical functions Paul Zimmermann
@ 2022-01-04  0:04 ` Damian McGuckin
  2022-01-04  2:29 ` Rich Felker
  1 sibling, 0 replies; 6+ messages in thread
From: Damian McGuckin @ 2022-01-04  0:04 UTC (permalink / raw)
  To: musl


Hi Paul,

On Mon, 3 Jan 2022, Paul Zimmermann wrote:

> Are developers of Musl interested by such functions?

Speaking only for myself, yes.

> If so, we could discuss what would be the requirements for integration 
> in Musl in terms of license, table size, allowed operations.

Outside my remit.  I am sure that only as you consider the needs of the 
various targets, will a common set of such requirements become apparent.

Stay safe - Damian

^ permalink raw reply	[flat|nested] 6+ messages in thread

* Re: [musl] correctly rounded mathematical functions
  2022-01-03 13:07 [musl] correctly rounded mathematical functions Paul Zimmermann
  2022-01-04  0:04 ` Damian McGuckin
@ 2022-01-04  2:29 ` Rich Felker
  2022-01-04 14:19   ` Paul Zimmermann
  1 sibling, 1 reply; 6+ messages in thread
From: Rich Felker @ 2022-01-04  2:29 UTC (permalink / raw)
  To: Paul Zimmermann; +Cc: musl, sibid, christoph.lauter, Jean-Michel.Muller

On Mon, Jan 03, 2022 at 02:07:59PM +0100, Paul Zimmermann wrote:
>       Dear Musl developers,
> 
> the current C working draft [1, p392] has reserved names for correctly
> rounded functions (cr_exp, cr_log, cr_sin, ...).
> 
> We propose to provide such correctly rounded implementations
> for the three IEEE formats (binary32, binary64, binary128) and the
> "extended double" format (long double on x86_64).
> 
> These implementations will 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);
> 
> 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.
> 
> Our objective is to provide open-source implementations that can be integrated
> in the major mathematical libraries (GNU libc, Intel Math Library, AMD Libm,
> Redhat Newlib, OpenLibm, Musl, llvm-libc, CUDA, ROCm).
> 
> Are developers of Musl interested by such functions?
> If so, we could discuss what would be the requirements for integration in
> Musl in terms of license, table size, allowed operations.

License: should be licensed under something permissive and
MIT-compatible, preferably explicitly MIT if you'll be
dual-/multi-licensing anyway.

Table size: main consideration is that tables need to be pure
shareable rodata, so no runtime generation of contents or pointers in
contents. Failure to follow this would produce significant cost in
*all processes* even ones that don't use the cr_* functions. I would
hope they'd also fit in a few tens of kB of rodata, but I don't know
how realistic that is.

Allowed operations: I'm not sure what the scope of this question is.
Are you asking if things like changing rounding mode would be okay? Or
something else. musl implements C11 and hopefully future versions of
the language library, but is implemented in C99 (+ very minimal set of
common/"GNU C" extensions), so code to be included in musl shouldn't
depend on new language features or compiler extensions not already
used (at least without checking on them first). We also have
softfloat-only targets and targets with excess precision, so code
should be compatible with those (which means not depending on changing
fenv; see commit 4f3d346bffdf9ed2b1803653643dc31242490944 for an
example of where this came up).

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

Is there a standalone version of the code where it can be read in full
not as a patch to glibc? Is the code being developed in such a way
that it's not potentially a derivative work of the glibc versions?

Rich

^ permalink raw reply	[flat|nested] 6+ messages in thread

* Re: [musl] correctly rounded mathematical functions
  2022-01-04  2:29 ` Rich Felker
@ 2022-01-04 14:19   ` Paul Zimmermann
  2022-01-04 21:12     ` Szabolcs Nagy
  0 siblings, 1 reply; 6+ messages in thread
From: Paul Zimmermann @ 2022-01-04 14:19 UTC (permalink / raw)
  To: Rich Felker; +Cc: musl, sibid, christoph.lauter, Jean-Michel.Muller

       Hi Rich,

> License: should be licensed under something permissive and
> MIT-compatible, preferably explicitly MIT if you'll be
> dual-/multi-licensing anyway.

ok, we'll consider multi-licensing with MIT.

> Table size: main consideration is that tables need to be pure
> shareable rodata, so no runtime generation of contents or pointers in
> contents. Failure to follow this would produce significant cost in
> *all processes* even ones that don't use the cr_* functions. I would
> hope they'd also fit in a few tens of kB of rodata, but I don't know
> how realistic that is.

tables so far are "static const"

> Allowed operations: I'm not sure what the scope of this question is.
> Are you asking if things like changing rounding mode would be okay? Or
> something else. musl implements C11 and hopefully future versions of
> the language library, but is implemented in C99 (+ very minimal set of
> common/"GNU C" extensions), so code to be included in musl shouldn't
> depend on new language features or compiler extensions not already
> used (at least without checking on them first). We also have
> softfloat-only targets and targets with excess precision, so code
> should be compatible with those (which means not depending on changing
> fenv; see commit 4f3d346bffdf9ed2b1803653643dc31242490944 for an
> example of where this came up).

we avoid playing with fesetround/fegetround since in general this makes
the code slower

> > 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].
> > 
> > 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/
> 
> Is there a standalone version of the code where it can be read in full
> not as a patch to glibc? Is the code being developed in such a way
> that it's not potentially a derivative work of the glibc versions?

yes there are several standalone versions of the code (entries marked "full"
on https://homepages.loria.fr/PZimmermann/CORE-MATH/).

Best regards,
Paul

^ permalink raw reply	[flat|nested] 6+ messages in thread

* Re: [musl] correctly rounded mathematical functions
  2022-01-04 14:19   ` Paul Zimmermann
@ 2022-01-04 21:12     ` Szabolcs Nagy
  2022-01-06  8:55       ` Paul Zimmermann
  0 siblings, 1 reply; 6+ messages in thread
From: Szabolcs Nagy @ 2022-01-04 21:12 UTC (permalink / raw)
  To: Paul Zimmermann
  Cc: Rich Felker, musl, sibid, christoph.lauter, Jean-Michel.Muller

* Paul Zimmermann <Paul.Zimmermann@inria.fr> [2022-01-04 15:19:48 +0100]:
> > Is there a standalone version of the code where it can be read in full
> > not as a patch to glibc? Is the code being developed in such a way
> > that it's not potentially a derivative work of the glibc versions?
> 
> yes there are several standalone versions of the code (entries marked "full"
> on https://homepages.loria.fr/PZimmermann/CORE-MATH/).

i only looked at cr_asinf:

  float
  cr_asinf (float x)
  {
    /* deal here with NaN, +Inf and -Inf */

yeah that can be tricky and different across math libs.
(errno vs no errno, wrapper to handle special cases vs no wrapper etc)

to minimize branches you may want to merge it with the |x|>1 check below.

    double absx = fabsf (x), y;
    if (absx > 1)
      return sqrt (-1.0); /* NaN */
    ...

this reminds me that musl does not use compiler builtins at build time,
which means fabsf, sqrt are extern calls, which means this will not be
a leaf function (note: sqrt is not a tail call here because of the
implicit conversion, i think it should be sqrtf).

i prefer to tail call a function with descriptive name when handling
errors, but it does not always work out well: e.g. underflow can be
tricky if you want to ensure that an exact subnormal result does not
raise it (e.g. powf(2, -140)), but inexact does.

^ permalink raw reply	[flat|nested] 6+ messages in thread

* Re: [musl] correctly rounded mathematical functions
  2022-01-04 21:12     ` Szabolcs Nagy
@ 2022-01-06  8:55       ` Paul Zimmermann
  0 siblings, 0 replies; 6+ messages in thread
From: Paul Zimmermann @ 2022-01-06  8:55 UTC (permalink / raw)
  To: Szabolcs Nagy; +Cc: dalias, musl, sibid, christoph.lauter, Jean-Michel.Muller

       Dear Szabolcs,

> Date: Tue, 4 Jan 2022 22:12:02 +0100
> From: Szabolcs Nagy <nsz@port70.net>
>
> i only looked at cr_asinf:

thank you for looking! This is much appreciated.

>   float
>   cr_asinf (float x)
>   {
>     /* deal here with NaN, +Inf and -Inf */
> 
> yeah that can be tricky and different across math libs.
> (errno vs no errno, wrapper to handle special cases vs no wrapper etc)

indeed, and that will be specific to each libm

> to minimize branches you may want to merge it with the |x|>1 check below.
> 
>     double absx = fabsf (x), y;
>     if (absx > 1)
>       return sqrt (-1.0); /* NaN */
>     ...

I tried to merge the detection of NaN, Inf with the case |x| > 1 using
the glibc macro asuint, but it was not faster (at least for glibc make bench).

> this reminds me that musl does not use compiler builtins at build time,
> which means fabsf, sqrt are extern calls, which means this will not be
> a leaf function (note: sqrt is not a tail call here because of the
> implicit conversion, i think it should be sqrtf).

then probably for musl it would be faster to have the following:

   double absx = (x > 0) ? x : -x;

we can probably define macros for this, say CORE_MATH_FABS, that would
be instantiated differently for each libm.

> i prefer to tail call a function with descriptive name when handling
> errors, but it does not always work out well: e.g. underflow can be
> tricky if you want to ensure that an exact subnormal result does not
> raise it (e.g. powf(2, -140)), but inexact does.

yes, sqrt (-1.0) was just a portable way to generate NaN. This code
will evolve over time with the feedback we get.

Best regards,
Paul

^ permalink raw reply	[flat|nested] 6+ messages in thread

end of thread, other threads:[~2022-01-06  8:56 UTC | newest]

Thread overview: 6+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2022-01-03 13:07 [musl] correctly rounded mathematical functions Paul Zimmermann
2022-01-04  0:04 ` Damian McGuckin
2022-01-04  2:29 ` Rich Felker
2022-01-04 14:19   ` Paul Zimmermann
2022-01-04 21:12     ` Szabolcs Nagy
2022-01-06  8:55       ` Paul Zimmermann

Code repositories for project(s) associated with this inbox:

	https://git.vuxu.org/mirror/musl/

This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for NNTP newsgroup(s).