[PATCH 0/4] Fix huge error for pow detected by Paul Zimmermann

[PATCH 0/4] Fix huge error for pow detected by Paul Zimmermann

[Andoni Arregi]

This patch series fixes in pow the huge error detected by Paul
Zimmermann where x is negative and y is a large odd integer.

There is also an accuracy fix for cases where x is close to 1 and y is
large.

Andoni Arregi (4):
  Fix missing sign for overflow/underflow where x is negative and y is
    large odd integer
  Fix x close to 1, y between 2^31 and 2^64
  Compare j as unsigned
  Replace always true if with else

 newlib/libm/math/e_pow.c | 48 +++++++++++++++++++++-------------------
 1 file changed, 25 insertions(+), 23 deletions(-)

-- 
2.40.0

[PATCH 1/4] Fix missing sign for overflow/underflow where x is negative and y is large odd integer

[Andoni Arregi]

---
 newlib/libm/math/e_pow.c | 37 ++++++++++++++++++-------------------
 1 file changed, 18 insertions(+), 19 deletions(-)

diff --git a/newlib/libm/math/e_pow.c b/newlib/libm/math/e_pow.c
index 258cca8dd..f02f76dc0 100644
--- a/newlib/libm/math/e_pow.c
+++ b/newlib/libm/math/e_pow.c
@@ -107,7 +107,7 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 #endif
 {
 	double z,ax,z_h,z_l,p_h,p_l;
-	double y1,t1,t2,r,s,t,u,v,w;
+	double y1,t1,t2,r,s,sign,t,u,v,w;
 	__int32_t i,j,k,yisint,n;
 	__int32_t hx,hy,ix,iy;
 	__uint32_t lx,ly;
@@ -184,23 +184,26 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 		return z;
 	    }
 	}
-    
+
     /* (x<0)**(non-int) is NaN */
-    /* REDHAT LOCAL: This used to be
-	if((((hx>>31)+1)|yisint)==0) return (x-x)/(x-x);
-       but ANSI C says a right shift of a signed negative quantity is
-       implementation defined.  */
-	if(((((__uint32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x);
+
+	n = ((uint32_t)hx >> 31U) - 1U;
+	if ((n | yisint) == 0) return (x-x)/(x-x);
+
+	sign = one; /* (sign of result -ve**odd) = -1 else = 1 */
+	if ((n | (yisint - 1)) == 0) {
+	    sign = -one;    /* (-ve)**(odd int) */
+        }
 
     /* |y| is huge */
 	if(iy>0x41e00000) { /* if |y| > 2**31 */
-	    if(iy>0x43f00000){	/* if |y| > 2**64, must o/uflow */
+	    if(iy>0x43f00000){	/* if |y| > 2**64, must o/uflow and y is an even integer */
 		if(ix<=0x3fefffff) return (hy<0)? __math_oflow(0):__math_uflow(0);
 		if(ix>=0x3ff00000) return (hy>0)? __math_oflow(0):__math_uflow(0);
 	    }
 	/* over/underflow if x is not close to one */
-	    if(ix<0x3fefffff) return (hy<0)? __math_oflow(0):__math_uflow(0);
-	    if(ix>0x3ff00000) return (hy>0)? __math_oflow(0):__math_uflow(0);
+	    if(ix<0x3fefffff) return (hy<0)? __math_oflow(sign<0):__math_uflow(sign<0);
+	    if(ix>0x3ff00000) return (hy>0)? __math_oflow(sign<0):__math_uflow(sign<0);
 	/* now |1-x| is tiny <= 2**-20, suffice to compute 
 	   log(x) by x-x^2/2+x^3/3-x^4/4 */
 	    t = ax-1;		/* t has 20 trailing zeros */
@@ -260,10 +263,6 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 	    t2 = z_l-(((t1-t)-dp_h[k])-z_h);
 	}
 
-	s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
-	if(((((__uint32_t)hx>>31)-1)|(yisint-1))==0)
-	    s = -one;/* (-ve)**(odd int) */
-
     /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
 	y1  = y;
 	SET_LOW_WORD(y1,0);
@@ -273,15 +272,15 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 	EXTRACT_WORDS(j,i,z);
 	if (j>=0x40900000) {				/* z >= 1024 */
 	    if(((j-0x40900000)|i)!=0)			/* if z > 1024 */
-		return __math_oflow(s<0);			/* overflow */
+		return __math_oflow(sign<0);			/* overflow */
 	    else {
-		if(p_l+ovt>z-p_h) return __math_oflow(s<0);	/* overflow */
+		if(p_l+ovt>z-p_h) return __math_oflow(sign<0);	/* overflow */
 	    }
 	} else if((j&0x7fffffff)>=0x4090cc00 ) {	/* z <= -1075 */
 	    if(((j-0xc090cc00)|i)!=0) 		/* z < -1075 */
-		return __math_uflow(s<0);		/* underflow */
+		return __math_uflow(sign<0);		/* underflow */
 	    else {
-		if(p_l<=z-p_h) return __math_uflow(s<0);	/* underflow */
+		if(p_l<=z-p_h) return __math_uflow(sign<0);	/* underflow */
 	    }
 	}
     /*
@@ -313,7 +312,7 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 	j += (n<<20);
 	if((j>>20)<=0) z = scalbn(z,(int)n);	/* subnormal output */
 	else SET_HIGH_WORD(z,j);
-	return s*z;
+	return sign*z;
 }
 
 #endif /* defined(_DOUBLE_IS_32BITS) */
-- 
2.40.0

[PATCH 2/4] Fix x close to 1, y between 2^31 and 2^64

[Andoni Arregi]

E.g. known errors:
x = 0x1.000002c5e2e99p+0, y = 0x1.c9eee35374af6p+31 had an error of 639
ULP and is now correctly rounded.
x = 0x1.fffffd2e3e669p-1, y = 0x1.344c9823eb66cp+32 had an error of 428
ULP and is now correctly rounded.
---
 newlib/libm/math/e_pow.c | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/newlib/libm/math/e_pow.c b/newlib/libm/math/e_pow.c
index f02f76dc0..f24a98f6e 100644
--- a/newlib/libm/math/e_pow.c
+++ b/newlib/libm/math/e_pow.c
@@ -196,8 +196,8 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
         }
 
     /* |y| is huge */
-	if(iy>0x41e00000) { /* if |y| > 2**31 */
-	    if(iy>0x43f00000){	/* if |y| > 2**64, must o/uflow and y is an even integer */
+	if(iy>0x42000000) { /* if |y| > ~2**33 (does not regard mantissa) */
+	    if(iy>0x43f00000){	/* if |y| > ~2**64, must o/uflow and y is an even integer */
 		if(ix<=0x3fefffff) return (hy<0)? __math_oflow(0):__math_uflow(0);
 		if(ix>=0x3ff00000) return (hy>0)? __math_oflow(0):__math_uflow(0);
 	    }
-- 
2.40.0

[PATCH 3/4] Compare j as unsigned

[Andoni Arregi]

j is int32_t and thus j<<31 is undefined if j==1.

Taken from FreeBSD commit bdd8abc6d6a93ce3ab8ad5db716222ee3110c4a3
---
 newlib/libm/math/e_pow.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/newlib/libm/math/e_pow.c b/newlib/libm/math/e_pow.c
index f24a98f6e..0e2251f10 100644
--- a/newlib/libm/math/e_pow.c
+++ b/newlib/libm/math/e_pow.c
@@ -141,7 +141,7 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
 		k = (iy>>20)-0x3ff;	   /* exponent */
 		if(k>20) {
 		    j = ly>>(52-k);
-		    if((j<<(52-k))==ly) yisint = 2-(j&1);
+		    if(((uint32_t)j<<(52-k))==ly) yisint = 2-(j&1);
 		} else if(ly==0) {
 		    j = iy>>(20-k);
 		    if((j<<(20-k))==iy) yisint = 2-(j&1);
-- 
2.40.0

[PATCH 4/4] Replace always true if with else

[Andoni Arregi]

---
 newlib/libm/math/e_pow.c | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/newlib/libm/math/e_pow.c b/newlib/libm/math/e_pow.c
index 0e2251f10..b21726007 100644
--- a/newlib/libm/math/e_pow.c
+++ b/newlib/libm/math/e_pow.c
@@ -198,8 +198,11 @@ ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
     /* |y| is huge */
 	if(iy>0x42000000) { /* if |y| > ~2**33 (does not regard mantissa) */
 	    if(iy>0x43f00000){	/* if |y| > ~2**64, must o/uflow and y is an even integer */
-		if(ix<=0x3fefffff) return (hy<0)? __math_oflow(0):__math_uflow(0);
-		if(ix>=0x3ff00000) return (hy>0)? __math_oflow(0):__math_uflow(0);
+		if(ix<=0x3fefffff) { /* |x| < 1 */
+		    return (hy<0)? __math_oflow(0):__math_uflow(0);
+		} else {             /* |x| >= 1 */
+		    return (hy>0)? __math_oflow(0):__math_uflow(0);
+		}
 	    }
 	/* over/underflow if x is not close to one */
 	    if(ix<0x3fefffff) return (hy<0)? __math_oflow(sign<0):__math_uflow(sign<0);
-- 
2.40.0

Re: [PATCH 0/4] Fix huge error for pow detected by Paul Zimmermann

[Paul Zimmermann]

       Dear Andoni,

this kind of feedback is very much appreciated.

I will test your patches next Monday.

Best regards,
Paul

> From: Andoni Arregi <andoni.arregui@gtd-gmbh.de>
> Cc: Andoni Arregi <andoni.arregui@gtd-gmbh.de>
> Date: Wed, 12 Apr 2023 17:34:41 +0200
> 
> This patch series fixes in pow the huge error detected by Paul
> Zimmermann where x is negative and y is a large odd integer.
> 
> There is also an accuracy fix for cases where x is close to 1 and y is
> large.
> 
> Andoni Arregi (4):
>   Fix missing sign for overflow/underflow where x is negative and y is
>     large odd integer
>   Fix x close to 1, y between 2^31 and 2^64
>   Compare j as unsigned
>   Replace always true if with else
> 
>  newlib/libm/math/e_pow.c | 48 +++++++++++++++++++++-------------------
>  1 file changed, 25 insertions(+), 23 deletions(-)
> 
> -- 
> 2.40.0
> 
> 
> 

Re: [PATCH 0/4] Fix huge error for pow detected by Paul Zimmermann

[Paul Zimmermann]

       Hi Andoni,

I had time to check this patch series on top of be2749c.

I confirm the issue I reported is fixed.

In addition, the maximal known error, which was 636 ulps
apart from the above issue, has dropped to less than 1 ulp:

pow 0 -1 0x1.4c3064d46aef3p-851,-0x1.cef4118ed50c2p-10 [0.893] 0.892119 0.8921192373561432

Thank you, great work!

Paul

> From: Andoni Arregi <andoni.arregui@gtd-gmbh.de>
> Cc: Andoni Arregi <andoni.arregui@gtd-gmbh.de>
> Date: Wed, 12 Apr 2023 17:34:41 +0200
> 
> This patch series fixes in pow the huge error detected by Paul
> Zimmermann where x is negative and y is a large odd integer.
> 
> There is also an accuracy fix for cases where x is close to 1 and y is
> large.
> 
> Andoni Arregi (4):
>   Fix missing sign for overflow/underflow where x is negative and y is
>     large odd integer
>   Fix x close to 1, y between 2^31 and 2^64
>   Compare j as unsigned
>   Replace always true if with else
> 
>  newlib/libm/math/e_pow.c | 48 +++++++++++++++++++++-------------------
>  1 file changed, 25 insertions(+), 23 deletions(-)
> 
> -- 
> 2.40.0
> 
> 
> 

Re: [PATCH 0/4] Fix huge error for pow detected by Paul Zimmermann

[Andoni Arregi]

Dear Paul,

thank *you* for detecting the errors/inaccuracies first (questions are as important as the answers) and then
independently confirming that the patch solves them.

You know how tricky this is and how easy it is to fix something here and break it there. That's why such an independent
verification in addition to our tests is so important :-)

By the way, most of the submitted work is product of the effort of our colleague Fabian Schriever.

Best regards

On Thu, 2023-04-13 at 14:16 +0200, Paul Zimmermann wrote:
>        Hi Andoni,
> 
> I had time to check this patch series on top of be2749c.
> 
> I confirm the issue I reported is fixed.
> 
> In addition, the maximal known error, which was 636 ulps
> apart from the above issue, has dropped to less than 1 ulp:
> 
> pow 0 -1 0x1.4c3064d46aef3p-851,-0x1.cef4118ed50c2p-10 [0.893] 0.892119 0.8921192373561432
> 
> Thank you, great work!
> 
> Paul
> 
> > From: Andoni Arregi <andoni.arregui@gtd-gmbh.de>
> > Cc: Andoni Arregi <andoni.arregui@gtd-gmbh.de>
> > Date: Wed, 12 Apr 2023 17:34:41 +0200
> > 
> > This patch series fixes in pow the huge error detected by Paul
> > Zimmermann where x is negative and y is a large odd integer.
> > 
> > There is also an accuracy fix for cases where x is close to 1 and y is
> > large.
> > 
> > Andoni Arregi (4):
> >   Fix missing sign for overflow/underflow where x is negative and y is
> >     large odd integer
> >   Fix x close to 1, y between 2^31 and 2^64
> >   Compare j as unsigned
> >   Replace always true if with else
> > 
> >  newlib/libm/math/e_pow.c | 48 +++++++++++++++++++++-------------------
> >  1 file changed, 25 insertions(+), 23 deletions(-)
> > 
> > -- 
> > 2.40.0
> > 
> > 
> > 

-- 
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 0/4] Fix huge error for pow detected by Paul Zimmermann

[Jeff Johnston]

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

Thanks Andoni and Paul.  Patches pushed to master.

-- Jeff J.

On Thu, Apr 13, 2023 at 10:42 AM Andoni Arregi <andoni.arregui@gtd-gmbh.de>
wrote:

> Dear Paul,
>
> thank *you* for detecting the errors/inaccuracies first (questions are as
> important as the answers) and then
> independently confirming that the patch solves them.
>
> You know how tricky this is and how easy it is to fix something here and
> break it there. That's why such an independent
> verification in addition to our tests is so important :-)
>
> By the way, most of the submitted work is product of the effort of our
> colleague Fabian Schriever.
>
> Best regards
>
> On Thu, 2023-04-13 at 14:16 +0200, Paul Zimmermann wrote:
> >        Hi Andoni,
> >
> > I had time to check this patch series on top of be2749c.
> >
> > I confirm the issue I reported is fixed.
> >
> > In addition, the maximal known error, which was 636 ulps
> > apart from the above issue, has dropped to less than 1 ulp:
> >
> > pow 0 -1 0x1.4c3064d46aef3p-851,-0x1.cef4118ed50c2p-10 [0.893] 0.892119
> 0.8921192373561432
> >
> > Thank you, great work!
> >
> > Paul
> >
> > > From: Andoni Arregi <andoni.arregui@gtd-gmbh.de>
> > > Cc: Andoni Arregi <andoni.arregui@gtd-gmbh.de>
> > > Date: Wed, 12 Apr 2023 17:34:41 +0200
> > >
> > > This patch series fixes in pow the huge error detected by Paul
> > > Zimmermann where x is negative and y is a large odd integer.
> > >
> > > There is also an accuracy fix for cases where x is close to 1 and y is
> > > large.
> > >
> > > Andoni Arregi (4):
> > >   Fix missing sign for overflow/underflow where x is negative and y is
> > >     large odd integer
> > >   Fix x close to 1, y between 2^31 and 2^64
> > >   Compare j as unsigned
> > >   Replace always true if with else
> > >
> > >  newlib/libm/math/e_pow.c | 48 +++++++++++++++++++++-------------------
> > >  1 file changed, 25 insertions(+), 23 deletions(-)
> > >
> > > --
> > > 2.40.0
> > >
> > >
> > >
>
> --
> 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
>
>