[PATCH] Fix hypotf missing mask in hi+lo decomposition

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

[PATCH] Fix hypotf missing mask in hi+lo decomposition

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


Reply | Threaded
Open this post in threaded view
|

Re: [PATCH] Fix hypotf missing mask in hi+lo decomposition

Sourceware - newlib list mailing list
Fabian Schriever <[hidden email]> 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

signature.asc (847 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: [PATCH] Fix hypotf missing mask in hi+lo decomposition

Sourceware - newlib list mailing list
In reply to this post by Fabian Schriever
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

signature.asc (849 bytes) Download Attachment