Fix hypotf missing mask in hi+lo decomposition

Message ID 20200319153408.1413-1-fabian.schriever@gtd-gmbh.de
State Accepted
Commit 6b0c1e7cc8e4b0568f00c8444fff8da67688add3
Headers show
Series
  • Fix hypotf missing mask in hi+lo decomposition
Related show

Commit Message

Fabian Schriever March 19, 2020, 3:34 p.m.
Add the missing mask for the decomposition of hi+lo which caused some
errors of 1-2 ULP.

This change is taken over from FreeBSD:
https://github.com/freebsd/freebsd/commit/95436ce20dab5a34ba46373410b96411b1734578

Additionally I've removed some variable assignments which were never
read before being overwritten again in the next 2 lines.
---
 newlib/libm/math/ef_hypot.c | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

-- 
2.24.1.windows.2

Comments

Jeff Johnston via Newlib March 19, 2020, 10:27 p.m. | #1
Fabian Schriever <fabian.schriever@gtd-gmbh.de> writes:

> Add the missing mask for the decomposition of hi+lo which caused some

> errors of 1-2 ULP.


> -	    SET_FLOAT_WORD(t1,ha+0x00800000L);

> +	    SET_FLOAT_WORD(t1,(ha+0x00800000L)&0xfffff000UL);


Looks correct to me. I note that glibc has switched to using double for
float paths, so there's no comparable code there.

-- 
-keith
Jeff Johnston via Newlib March 20, 2020, 9:16 a.m. | #2
On Mar 19 16:34, Fabian Schriever wrote:
> Add the missing mask for the decomposition of hi+lo which caused some

> errors of 1-2 ULP.

> 

> This change is taken over from FreeBSD:

> https://github.com/freebsd/freebsd/commit/95436ce20dab5a34ba46373410b96411b1734578

> 

> Additionally I've removed some variable assignments which were never

> read before being overwritten again in the next 2 lines.

> ---

>  newlib/libm/math/ef_hypot.c | 4 ++--

>  1 file changed, 2 insertions(+), 2 deletions(-)

> 

> diff --git a/newlib/libm/math/ef_hypot.c b/newlib/libm/math/ef_hypot.c

> index 9368eb41c..a70c92b88 100644

> --- a/newlib/libm/math/ef_hypot.c

> +++ b/newlib/libm/math/ef_hypot.c

> @@ -29,7 +29,7 @@

>  	ha &= 0x7fffffffL;

>  	GET_FLOAT_WORD(hb,y);

>  	hb &= 0x7fffffffL;

> -	if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}

> +	if(hb > ha) { j = ha; ha = hb; hb = j; }

>  	SET_FLOAT_WORD(a,ha);	/* a <- |a| */

>  	SET_FLOAT_WORD(b,hb);	/* b <- |b| */

>  	if((ha-hb)>0xf000000L) {return a+b;} /* x/y > 2**30 */

> @@ -72,7 +72,7 @@

>  	    a  = a+a;

>  	    SET_FLOAT_WORD(y1,hb&0xfffff000L);

>  	    y2 = b - y1;

> -	    SET_FLOAT_WORD(t1,ha+0x00800000L);

> +	    SET_FLOAT_WORD(t1,(ha+0x00800000L)&0xfffff000UL);

>  	    t2 = a - t1;

>  	    w  = __ieee754_sqrtf(t1*y1-(w*(-w)-(t1*y2+t2*b)));

>  	}

> -- 

> 2.24.1.windows.2

> 


Pushed.


Thanks,
Corinna

-- 
Corinna Vinschen
Cygwin Maintainer
Red Hat

Patch

diff --git a/newlib/libm/math/ef_hypot.c b/newlib/libm/math/ef_hypot.c
index 9368eb41c..a70c92b88 100644
--- a/newlib/libm/math/ef_hypot.c
+++ b/newlib/libm/math/ef_hypot.c
@@ -29,7 +29,7 @@ 
 	ha &= 0x7fffffffL;
 	GET_FLOAT_WORD(hb,y);
 	hb &= 0x7fffffffL;
-	if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
+	if(hb > ha) { j = ha; ha = hb; hb = j; }
 	SET_FLOAT_WORD(a,ha);	/* a <- |a| */
 	SET_FLOAT_WORD(b,hb);	/* b <- |b| */
 	if((ha-hb)>0xf000000L) {return a+b;} /* x/y > 2**30 */
@@ -72,7 +72,7 @@ 
 	    a  = a+a;
 	    SET_FLOAT_WORD(y1,hb&0xfffff000L);
 	    y2 = b - y1;
-	    SET_FLOAT_WORD(t1,ha+0x00800000L);
+	    SET_FLOAT_WORD(t1,(ha+0x00800000L)&0xfffff000UL);
 	    t2 = a - t1;
 	    w  = __ieee754_sqrtf(t1*y1-(w*(-w)-(t1*y2+t2*b)));
 	}