correctly rounded mathematical functions

correctly rounded mathematical functions

[Paul Zimmermann]

      Dear Newlib 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 Newlib interested by such functions?
If so, we could discuss what would be the requirements for integration in
Newlib 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/

Re: correctly rounded mathematical functions

[Joel Sherrill]

On Mon, Jan 3, 2022 at 6:58 AM Paul Zimmermann <Paul.Zimmermann@inria.fr> wrote:
>
>       Dear Newlib 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 Newlib interested by such functions?
> If so, we could discuss what would be the requirements for integration in
> Newlib in terms of license, table size, allowed operations.

Speaking from the RTEMS perspective, we are very interested
in the addition of more POSIX and C Standard Library methods.
We have been having GSoC students add them where possible
for a few years.

The license has to be permissive and should have no advertising
although based on COPYING.NEWLIB, some must have a BSD
advertising clause still. Possibly those need review.

https://sourceware.org/git/?p=newlib-cygwin.git;a=blob_plain;f=COPYING.NEWLIB;hb=HEAD

As long as the method's footprint only impacts applications that use
it, there shouldn't be a huge concern on size but the target domain is
mostly smaller single process, no-MMU embedded systems. That is
ignoring Cygwin which has fewer constraints.  If you end up adding
megabytes, that's going to be bad in general.



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

Great! What's the size profile on those?

What's the minimum compiler or C language version to build these?

--joel

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

[Paul Zimmermann]

       Dear Joel,

thank you for your answer.

> From: Joel Sherrill <joel@rtems.org>
> Date: Mon, 3 Jan 2022 14:41:20 -0600
> 
> On Mon, Jan 3, 2022 at 6:58 AM Paul Zimmermann <Paul.Zimmermann@inria.fr> wrote:
> >
> >       Dear Newlib 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 Newlib interested by such functions?
> > If so, we could discuss what would be the requirements for integration in
> > Newlib in terms of license, table size, allowed operations.
> 
> Speaking from the RTEMS perspective, we are very interested
> in the addition of more POSIX and C Standard Library methods.
> We have been having GSoC students add them where possible
> for a few years.

great!

> The license has to be permissive and should have no advertising
> although based on COPYING.NEWLIB, some must have a BSD
> advertising clause still. Possibly those need review.
> 
> https://sourceware.org/git/?p=newlib-cygwin.git;a=blob_plain;f=COPYING.NEWLIB;hb=HEAD

since some Newlib files are already under LGPL v3+, for example
include/cgen/basic-ops.h, I guess the current
license we have (LGPL v3+) should be ok for Newlib.

> As long as the method's footprint only impacts applications that use
> it, there shouldn't be a huge concern on size but the target domain is
> mostly smaller single process, no-MMU embedded systems. That is
> ignoring Cygwin which has fewer constraints.  If you end up adding
> megabytes, that's going to be bad in general.

we'll try to have reasonable table sizes. And maybe for some functions
we'll provide several implementations, with different table sizes.

> > 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].
> 
> Great! What's the size profile on those?

the binary32 cbrt code is quite small:

-rw-r--r-- 1 zimmerma caramba   2766 Dec 23 14:27 cbrtf.c

the binary128 code is larger:

-rw-r--r-- 1 zimmerma caramba  38450 Dec 22 18:06 cbrtf128-v2.patch

> What's the minimum compiler or C language version to build these?

good question. So far we work with current gcc versions (for example 11.2.0)
but did not exercise with older version. We don't use advanced C feature,
except the "glibc patches" use some GNU libc macros (for example mul_split
to split the product of two doubles into a+b, fast_two_sum, ...).

Best regards,
Paul

Re: correctly rounded mathematical functions

[Joel Sherrill]

On Tue, Jan 4, 2022 at 3:44 AM Paul Zimmermann <Paul.Zimmermann@inria.fr> wrote:
>
>        Dear Joel,
>
> thank you for your answer.
>
> > From: Joel Sherrill <joel@rtems.org>
> > Date: Mon, 3 Jan 2022 14:41:20 -0600
> >
> > On Mon, Jan 3, 2022 at 6:58 AM Paul Zimmermann <Paul.Zimmermann@inria.fr> wrote:
> > >
> > >       Dear Newlib 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 Newlib interested by such functions?
> > > If so, we could discuss what would be the requirements for integration in
> > > Newlib in terms of license, table size, allowed operations.
> >
> > Speaking from the RTEMS perspective, we are very interested
> > in the addition of more POSIX and C Standard Library methods.
> > We have been having GSoC students add them where possible
> > for a few years.
>
> great!
>
> > The license has to be permissive and should have no advertising
> > although based on COPYING.NEWLIB, some must have a BSD
> > advertising clause still. Possibly those need review.
> >
> > https://sourceware.org/git/?p=newlib-cygwin.git;a=blob_plain;f=COPYING.NEWLIB;hb=HEAD
>
> since some Newlib files are already under LGPL v3+, for example
> include/cgen/basic-ops.h, I guess the current
> license we have (LGPL v3+) should be ok for Newlib.

I think those files are all under a linux/ or rdos/ subdirectory and do not
impact any other target unless they use iconvdata/. I grep'ed and didn't
see anything else. Jeff or Corrina should speak up and provide a
definitive answer.

> > As long as the method's footprint only impacts applications that use
> > it, there shouldn't be a huge concern on size but the target domain is
> > mostly smaller single process, no-MMU embedded systems. That is
> > ignoring Cygwin which has fewer constraints.  If you end up adding
> > megabytes, that's going to be bad in general.
>
> we'll try to have reasonable table sizes. And maybe for some functions
> we'll provide several implementations, with different table sizes.

Great! It's hard to comment on size without something to complain about. :)

> > > 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].
> >
> > Great! What's the size profile on those?
>
> the binary32 cbrt code is quite small:
>
> -rw-r--r-- 1 zimmerma caramba   2766 Dec 23 14:27 cbrtf.c
>
> the binary128 code is larger:
>
> -rw-r--r-- 1 zimmerma caramba  38450 Dec 22 18:06 cbrtf128-v2.patch

How is binary128 implemented on target CPUs without native support?
Is this something GCC will include soft methods for?

> > What's the minimum compiler or C language version to build these?
>
> good question. So far we work with current gcc versions (for example 11.2.0)
> but did not exercise with older version. We don't use advanced C feature,
> except the "glibc patches" use some GNU libc macros (for example mul_split
> to split the product of two doubles into a+b, fast_two_sum, ...).

I didn't expect anything too fancy. Will they compile with -std=c99, c11? newlib
generally provides all methods but the header file guards protect the prototype
visibility.

glibc specific macros can be reimplemented. I assume that's going to be an
issue on other non-glibc libraries.  Hopefully that's not too hard to deal with.

--joel

> Best regards,
> Paul

Re: correctly rounded mathematical functions

[Paul Zimmermann]

       Hi Joel,

> How is binary128 implemented on target CPUs without native support?
> Is this something GCC will include soft methods for?

we use the _Float128 type. For internal computations we avoid it as much as
possible, using integer types, which yields a speedup of up to 10 with respect
to the glibc code.

> I didn't expect anything too fancy. Will they compile with -std=c99, c11? newlib
> generally provides all methods but the header file guards protect the prototype
> visibility.

I'll check if the code compiles with -std=c99, c11.

Best regards,
Paul

Re: correctly rounded mathematical functions

[Jeff Johnston]

See below.

On Tue, Jan 4, 2022 at 10:23 AM Joel Sherrill <joel@rtems.org> wrote:

> On Tue, Jan 4, 2022 at 3:44 AM Paul Zimmermann <Paul.Zimmermann@inria.fr>
> wrote:
> >
> >        Dear Joel,
> >
> > thank you for your answer.
> >
> > > From: Joel Sherrill <joel@rtems.org>
> > > Date: Mon, 3 Jan 2022 14:41:20 -0600
> > >
> > > On Mon, Jan 3, 2022 at 6:58 AM Paul Zimmermann <
> Paul.Zimmermann@inria.fr> wrote:
> > > >
> > > >       Dear Newlib 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 Newlib interested by such functions?
> > > > If so, we could discuss what would be the requirements for
> integration in
> > > > Newlib in terms of license, table size, allowed operations.
> > >
> > > Speaking from the RTEMS perspective, we are very interested
> > > in the addition of more POSIX and C Standard Library methods.
> > > We have been having GSoC students add them where possible
> > > for a few years.
> >
> > great!
> >
> > > The license has to be permissive and should have no advertising
> > > although based on COPYING.NEWLIB, some must have a BSD
> > > advertising clause still. Possibly those need review.
> > >
> > >
> https://sourceware.org/git/?p=newlib-cygwin.git;a=blob_plain;f=COPYING.NEWLIB;hb=HEAD
> >
> > since some Newlib files are already under LGPL v3+, for example
> > include/cgen/basic-ops.h, I guess the current
> > license we have (LGPL v3+) should be ok for Newlib.
>
> I think those files are all under a linux/ or rdos/ subdirectory and do not
> impact any other target unless they use iconvdata/. I grep'ed and didn't
> see anything else. Jeff or Corrina should speak up and provide a
> definitive answer.
>
>
Shared newlib code cannot have an LGPL license.  Our Linux code is for a
shared library
and uses code from glibc (including depending on glibc headers) so it is
permitted, but the code
is not shared with other platforms.

> > As long as the method's footprint only impacts applications that use
> > > it, there shouldn't be a huge concern on size but the target domain is
> > > mostly smaller single process, no-MMU embedded systems. That is
> > > ignoring Cygwin which has fewer constraints.  If you end up adding
> > > megabytes, that's going to be bad in general.
> >
> > we'll try to have reasonable table sizes. And maybe for some functions
> > we'll provide several implementations, with different table sizes.
>
> Great! It's hard to comment on size without something to complain about. :)
>
> > > > 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].
> > >
> > > Great! What's the size profile on those?
> >
> > the binary32 cbrt code is quite small:
> >
> > -rw-r--r-- 1 zimmerma caramba   2766 Dec 23 14:27 cbrtf.c
> >
> > the binary128 code is larger:
> >
> > -rw-r--r-- 1 zimmerma caramba  38450 Dec 22 18:06 cbrtf128-v2.patch
>
> How is binary128 implemented on target CPUs without native support?
> Is this something GCC will include soft methods for?
>
> > > What's the minimum compiler or C language version to build these?
> >
> > good question. So far we work with current gcc versions (for example
> 11.2.0)
> > but did not exercise with older version. We don't use advanced C feature,
> > except the "glibc patches" use some GNU libc macros (for example
> mul_split
> > to split the product of two doubles into a+b, fast_two_sum, ...).
>
> I didn't expect anything too fancy. Will they compile with -std=c99, c11?
> newlib
> generally provides all methods but the header file guards protect the
> prototype
> visibility.
>
> glibc specific macros can be reimplemented. I assume that's going to be an
> issue on other non-glibc libraries.  Hopefully that's not too hard to deal
> with.
>
> --joel
>
> > Best regards,
> > Paul
>
>

Re: correctly rounded mathematical functions

[Keith Packard]

[-- Attachment #1: Type: text/plain, Size: 721 bytes --]

Joel Sherrill <joel@rtems.org> writes:

> I think those files are all under a linux/ or rdos/ subdirectory and do not
> impact any other target unless they use iconvdata/. I grep'ed and didn't
> see anything else. Jeff or Corrina should speak up and provide a
> definitive answer.

Yeah, I've forked newlib into 'picolibc' and have removed all of the
non-BSD licensed files to make the license status of that project
clearly BSD only. All of the non-BSD files are used to build the library
only for specific targets; the general library code is all BSD licensed.

Something BSD-ish would preserve the license status of newlib built for
embedded targets, and would also be usable by picolibc.

-- 
-keith

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]

Re: correctly rounded mathematical functions

[Paul Zimmermann]

       Hi Jeff,

> Shared newlib code cannot have an LGPL license.  Our Linux code is for a
> shared library
> and uses code from glibc (including depending on glibc headers) so it is
> permitted, but the code
> is not shared with other platforms.

ok. Would the MIT license be ok for newlib? We are open to multi-licensing,
but we'll try to minimize the number of licenses (the MIT license would be
ok for Musl).

Best regards,
Paul

Re: correctly rounded mathematical functions

[Paul Zimmermann]

       Hi Keith,

> Yeah, I've forked newlib into 'picolibc' and have removed all of the
> non-BSD licensed files to make the license status of that project
> clearly BSD only. All of the non-BSD files are used to build the library
> only for specific targets; the general library code is all BSD licensed.
> 
> Something BSD-ish would preserve the license status of newlib built for
> embedded targets, and would also be usable by picolibc.

then I guess the MIT license would be ok for both newlib and picolibc,
since it is more permissive than the BSD license.

Best regards,
Paul

Re: correctly rounded mathematical functions

[Joel Sherrill]

I hate to ask since I know this is unrelated but what is the best source for
properly licensed and accurate long double methods. I am still hoping that
part of an RTEMS GSoC project is to merge them from somewhere into
newlib.

And I'm thrilled at all the autoconf/automake work. That has really
been a challenge
for students.

On Wed, Jan 5, 2022 at 5:39 AM Paul Zimmermann <Paul.Zimmermann@inria.fr> wrote:
>
>        Hi Keith,
>
> > Yeah, I've forked newlib into 'picolibc' and have removed all of the
> > non-BSD licensed files to make the license status of that project
> > clearly BSD only. All of the non-BSD files are used to build the library
> > only for specific targets; the general library code is all BSD licensed.
> >
> > Something BSD-ish would preserve the license status of newlib built for
> > embedded targets, and would also be usable by picolibc.
>
> then I guess the MIT license would be ok for both newlib and picolibc,
> since it is more permissive than the BSD license.
>
> Best regards,
> Paul

Re: correctly rounded mathematical functions

[Keith Packard]

[-- Attachment #1: Type: text/plain, Size: 233 bytes --]

Paul Zimmermann <Paul.Zimmermann@inria.fr> writes:

> then I guess the MIT license would be ok for both newlib and picolibc,
> since it is more permissive than the BSD license.

Yup, MIT is fully BSD compatible.

-- 
-keith

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]

Re: correctly rounded mathematical functions

[Paul Zimmermann]

       Dear Joel,

> I hate to ask since I know this is unrelated but what is the best source for
> properly licensed and accurate long double methods. I am still hoping that
> part of an RTEMS GSoC project is to merge them from somewhere into
> newlib.

as far as I know, four libraries provide long double functions:
GNU libc, Intel Math Library, OpenLibm and Musl. See Table 11 from
https://members.loria.fr/PZimmermann/papers/accuracy.pdf.

And of course we also plan to provide (correctly rounded) long double
functions, acos and cbrt are already available on the CORE-MATH page.

> And I'm thrilled at all the autoconf/automake work. That has really
> been a challenge
> for students.

Best regards,
Paul

Re: correctly rounded mathematical functions

[Joel Sherrill]

On Thu, Jan 6, 2022 at 4:12 AM Paul Zimmermann <Paul.Zimmermann@inria.fr> wrote:
>
>        Dear Joel,
>
> > I hate to ask since I know this is unrelated but what is the best source for
> > properly licensed and accurate long double methods. I am still hoping that
> > part of an RTEMS GSoC project is to merge them from somewhere into
> > newlib.
>
> as far as I know, four libraries provide long double functions:
> GNU libc, Intel Math Library, OpenLibm and Musl. See Table 11 from
> https://members.loria.fr/PZimmermann/papers/accuracy.pdf.

Is the FreeBSD long double support related to any of those libraries?

https://github.com/freebsd/freebsd-src/tree/main/lib/msun

Any comments on their implementation?

> And of course we also plan to provide (correctly rounded) long double
> functions, acos and cbrt are already available on the CORE-MATH page.
>
> > And I'm thrilled at all the autoconf/automake work. That has really
> > been a challenge
> > for students.
>
> Best regards,
> Paul

Re: correctly rounded mathematical functions

[Paul Zimmermann]

       Dear Joel,

> > as far as I know, four libraries provide long double functions:
> > GNU libc, Intel Math Library, OpenLibm and Musl. See Table 11 from
> > https://members.loria.fr/PZimmermann/papers/accuracy.pdf.
> 
> Is the FreeBSD long double support related to any of those libraries?
> 
> https://github.com/freebsd/freebsd-src/tree/main/lib/msun

after a quick look at the cosl code
(https://github.com/freebsd/freebsd-src/blob/main/lib/msun/ld80/k_cosl.c)
it seems to be the same as the OpenLibm implementation
(openlibm-0.7.5/ld80/k_cosl.c).

> Any comments on their implementation?

the OpenLibm implementation has huge errors (up to 533 ulps for powl),
but they have fixed some of the issues I have reported to them.

See also https://github.com/JuliaMath/openlibm/issues/236.

Best regards,
Paul