regression in Newlib 4.2.0

regression in Newlib 4.2.0

[Paul Zimmermann]

       Hi,

I noticed a regression in Newlib 4.2.0:

Checking exp with newlib-4.2.0.20211231
MPFR library: 4.1.0
MPFR header:  4.1.0 (based on 4.1.0)
libm wrong by up to 1.94e+04 ulp(s) for x=0x1.62e302p+6
exp      gives inf
mpfr_exp gives 0x1.ff691ep+127
Total: errors=17982998 (0.42%) errors2=151 maxerr=1.94e+04 ulp(s)

The largest error for expf in Newlib 4.1.0 was 0.911 ulps
(https://members.loria.fr/PZimmermann/papers/accuracy.pdf).

Best regards,
Paul

[PATCH] Fix expf overflow limit

[Andoni Arregi]

Correct the overflow limit in the variable o_threshold to be consistent
with the FLT_UWORD_LOG_MAX variable used by the internal implementation
of the expf algorithm itself.
The u_threshold variable has also been modified to be written in the
same format.

Note that this fix improves the situation but does not completely
correct the inconsistencies regarding the overflow and underflow limits
between the expf wrapper (wf_exp.c) and the expf algorithm itself
(ef_exp.c).

Currently these limits are different for the
_FLT_LARGEST_EXPONENT_IS_NORMAL and _FLT_NO_DENORMALS cases as well as
for the case where __OBSOLETE_MATH is not defined (only for the
underflow limit in this case).
---
 newlib/libm/math/wf_exp.c | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/newlib/libm/math/wf_exp.c b/newlib/libm/math/wf_exp.c
index 38cacd95a..075b5bb0a 100644
--- a/newlib/libm/math/wf_exp.c
+++ b/newlib/libm/math/wf_exp.c
@@ -26,8 +26,8 @@ static const float
 #else
 static float
 #endif
-o_threshold=  8.8721679688e+01,  /* 0x42b17180 */
-u_threshold= -1.0397208405e+02;  /* 0xc2cff1b5 */
+o_threshold=  0x1.62e42ep+06,  /* 0x42b17217 */
+u_threshold= -0x1.9fe36ap+06;  /* 0xc2cff1b5 */
 
 #ifdef __STDC__
 	float expf(float x)		/* wrapper expf */
-- 
2.35.1

Re: regression in Newlib 4.2.0

[Andoni Arregi]

Dear Paul,

although we could fix the problem you spotted, we don't really got to
the root cause of the regression, i.e., this is not something that was
introduced in Newlib 4.2.0 and was not present in 4.1.0. wf_expf.c
didn't change in the meantime at all.

The only explanation I have is that it may have to do with the fact
that what you observed happens with the default compilation of Newlib.
If you compile with -D_IEEE_LIBM, then you will not see this
regression, as the wrappers for error handling disappear and with them
the problem.
It is a pity that for such a detailed analysis as the one you do, such
a small change in the build mechanism yields you completely different
implementations of the mathematical functions. I don't know whether
Newlib 4.1.0 was built by default with _IEEE_LIBM defined.

By the way, as we finished the work we have been carrying out for the
European Space Agency to qualify a mathematical library for critical
systems (which we based on the Newlib libm), we will now start
submitting a bunch of improvement patches for the Newlib libm.

We presented the LibmCS for the first time during the current Flight
Software Workshop 2022
(http://flightsoftware.jhuapl.edu/workshop/FSW2022)
https://www.youtube.com/watch?v=UcRefVlVL_A&list=PLK-T7jljJ6zbk0PshIwofgE2hiakBufJF&index=21

Best regards
  
On Fri, 2022-02-04 at 05:36 +0100, Paul Zimmermann wrote:
>        Hi,
> 
> I noticed a regression in Newlib 4.2.0:
> 
> Checking exp with newlib-4.2.0.20211231
> MPFR library: 4.1.0
> MPFR header:  4.1.0 (based on 4.1.0)
> libm wrong by up to 1.94e+04 ulp(s) for x=0x1.62e302p+6
> exp      gives inf
> mpfr_exp gives 0x1.ff691ep+127
> Total: errors=17982998 (0.42%) errors2=151 maxerr=1.94e+04 ulp(s)
> 
> The largest error for expf in Newlib 4.1.0 was 0.911 ulps
> (https://members.loria.fr/PZimmermann/papers/accuracy.pdf).
> 
> Best regards,
> Paul

-- 
Andoni Arregi
Geschäftsführer
 
GTD GmbH
Ravensburger Str. 32a, 88677 Markdorf
T: +49 7544 96440 22 | M: +49 151 65620499 | F: +49 7544 96440 29
http://www.gtd-gmbh.de andoni.arregi@gtd-gmbh.de

Re: [PATCH] Fix expf overflow limit

[Corinna Vinschen]

On Feb  9 14:52, Andoni Arregi wrote:
> Correct the overflow limit in the variable o_threshold to be consistent
> with the FLT_UWORD_LOG_MAX variable used by the internal implementation
> of the expf algorithm itself.
> The u_threshold variable has also been modified to be written in the
> same format.
> 
> Note that this fix improves the situation but does not completely
> correct the inconsistencies regarding the overflow and underflow limits
> between the expf wrapper (wf_exp.c) and the expf algorithm itself
> (ef_exp.c).
> 
> Currently these limits are different for the
> _FLT_LARGEST_EXPONENT_IS_NORMAL and _FLT_NO_DENORMALS cases as well as
> for the case where __OBSOLETE_MATH is not defined (only for the
> underflow limit in this case).
> ---
>  newlib/libm/math/wf_exp.c | 4 ++--
>  1 file changed, 2 insertions(+), 2 deletions(-)
> 
> diff --git a/newlib/libm/math/wf_exp.c b/newlib/libm/math/wf_exp.c
> index 38cacd95a..075b5bb0a 100644
> --- a/newlib/libm/math/wf_exp.c
> +++ b/newlib/libm/math/wf_exp.c
> @@ -26,8 +26,8 @@ static const float
>  #else
>  static float
>  #endif
> -o_threshold=  8.8721679688e+01,  /* 0x42b17180 */
> -u_threshold= -1.0397208405e+02;  /* 0xc2cff1b5 */
> +o_threshold=  0x1.62e42ep+06,  /* 0x42b17217 */
> +u_threshold= -0x1.9fe36ap+06;  /* 0xc2cff1b5 */
>  
>  #ifdef __STDC__
>  	float expf(float x)		/* wrapper expf */
> -- 
> 2.35.1

Pushed.


Thanks,
Corinna

Re: [PATCH] Fix expf overflow limit

[Paul Zimmermann]

       Dear Andoni,

I confirm that this fixes the expf large errors. With commit 8f2bd6f I now get
a maximal error of 0.911 ulps:

Checking exp with newlib-4.2.0.20211231
Using RedHat newlib
MPFR library: 4.1.0       
MPFR header:  4.1.0 (based on 4.1.0)
Checking function expf with MPFR_RNDN
libm wrong by up to 9.11e-01 ulp(s) [1] for x=-0x1.2e35b6p+5
exp      gives 0x1.6a2544p-55
mpfr_exp gives 0x1.6a2546p-55
Total: errors=17982848 (0.42%) errors2=0 maxerr=9.11e-01 ulp(s)

Thank you,
Paul

PS: in this case and in the two other issues I reported (about lgammaf and tgammaf),
I have to apologize. This was not a regression in Newlib 4.2.0, but a bug in our
checking program which did not report such issues. With the new version of our
program, we get similar results with Newlib 4.1.0 and Newlib 4.2.0.

> From: Andoni Arregi <andoni.arregui@gtd-gmbh.de>
> Cc: Paul Zimmermann <Paul.Zimmermann@inria.fr>
> Date: Wed, 09 Feb 2022 14:52:07 +0100
> Organization: GTD GmbH
> User-Agent: Evolution 3.42.3 
> 
> Correct the overflow limit in the variable o_threshold to be consistent
> with the FLT_UWORD_LOG_MAX variable used by the internal implementation
> of the expf algorithm itself.
> The u_threshold variable has also been modified to be written in the
> same format.
> 
> Note that this fix improves the situation but does not completely
> correct the inconsistencies regarding the overflow and underflow limits
> between the expf wrapper (wf_exp.c) and the expf algorithm itself
> (ef_exp.c).
> 
> Currently these limits are different for the
> _FLT_LARGEST_EXPONENT_IS_NORMAL and _FLT_NO_DENORMALS cases as well as
> for the case where __OBSOLETE_MATH is not defined (only for the
> underflow limit in this case).
> ---
>  newlib/libm/math/wf_exp.c | 4 ++--
>  1 file changed, 2 insertions(+), 2 deletions(-)
> 
> diff --git a/newlib/libm/math/wf_exp.c b/newlib/libm/math/wf_exp.c
> index 38cacd95a..075b5bb0a 100644
> --- a/newlib/libm/math/wf_exp.c
> +++ b/newlib/libm/math/wf_exp.c
> @@ -26,8 +26,8 @@ static const float
>  #else
>  static float
>  #endif
> -o_threshold=  8.8721679688e+01,  /* 0x42b17180 */
> -u_threshold= -1.0397208405e+02;  /* 0xc2cff1b5 */
> +o_threshold=  0x1.62e42ep+06,  /* 0x42b17217 */
> +u_threshold= -0x1.9fe36ap+06;  /* 0xc2cff1b5 */
>  
>  #ifdef __STDC__
>  	float expf(float x)		/* wrapper expf */
> -- 
> 2.35.1
> 
> 

regression in Newlib 4.2.0

[Paul Zimmermann]

       Hi,

another regression in Newlib 4.2.0 (with 4.1.0 the largest error for tgammaf
was 239 ulps according to
https://members.loria.fr/PZimmermann/papers/accuracy.pdf):

Checking tgamma with newlib-4.2.0.20211231
MPFR library: 4.1.0       
MPFR header:  4.1.0 (based on 4.1.0)
Checking function tgammaf with MPFR_RNDN
libm wrong by up to 1.55e+26 ulp(s) [18446744073709551615] for x=-0x1.000002p-65
tgamma      gives -inf
mpfr_tgamma gives -0x1.fffffcp+64
Total: errors=2967689031 (69.37%) errors2=2773050476 maxerr=1.55e+26 ulp(s)

Paul