[3/9] math: Simplify hypotf implementation

Message ID 20211006180557.933826-4-adhemerval.zanella@linaro.org
State New
Headers show
Series
  • Improve hypot()
Related show

Commit Message

Carlos O'Donell via Libc-alpha Oct. 6, 2021, 6:05 p.m.
Use isnan()/isinf() instead of GET_FLOAT_WORD and interger operations.
There is also no need to check for 0.0.

The file Copyright is also change to use  GPL, the implementation was
complete change by 7c10fd3515f to use double precision instead of
scaling and this change removes all the GET_FLOAT_WORD usage.

Checked on x86_64-linux-gnu.
---
 sysdeps/ieee754/flt-32/e_hypotf.c | 57 +++++++++++++------------------
 1 file changed, 23 insertions(+), 34 deletions(-)

-- 
2.30.2

Comments

Paul Zimmermann Oct. 7, 2021, 9:44 a.m. | #1
Dear Adhemerval,

> Use isnan()/isinf() instead of GET_FLOAT_WORD and interger operations.


interger -> integer

> There is also no need to check for 0.0.

> 

> The file Copyright is also change to use  GPL, the implementation was


change -> changed

> complete change by 7c10fd3515f to use double precision instead of


complete change -> completely changed ?

> scaling and this change removes all the GET_FLOAT_WORD usage.

> 

> Checked on x86_64-linux-gnu.

> ---

>  sysdeps/ieee754/flt-32/e_hypotf.c | 57 +++++++++++++------------------

>  1 file changed, 23 insertions(+), 34 deletions(-)

> 

> diff --git a/sysdeps/ieee754/flt-32/e_hypotf.c b/sysdeps/ieee754/flt-32/e_hypotf.c

> index e770947dc1..6495a91cd4 100644

> --- a/sysdeps/ieee754/flt-32/e_hypotf.c

> +++ b/sysdeps/ieee754/flt-32/e_hypotf.c

> @@ -1,46 +1,35 @@

> -/* e_hypotf.c -- float version of e_hypot.c.

> - */

> +/* Euclidean distance function.  Float/Binary32 version.

> +   Copyright (C) 2012-2021 Free Software Foundation, Inc.

> +   This file is part of the GNU C Library.

>  

> -/*

> - * ====================================================

> - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.

> - *

> - * Developed at SunPro, a Sun Microsystems, Inc. business.

> - * Permission to use, copy, modify, and distribute this

> - * software is freely granted, provided that this notice

> - * is preserved.

> - * ====================================================

> - */

> +   The GNU C Library is free software; you can redistribute it and/or

> +   modify it under the terms of the GNU Lesser General Public

> +   License as published by the Free Software Foundation; either

> +   version 2.1 of the License, or (at your option) any later version.

> +

> +   The GNU C Library is distributed in the hope that it will be useful,

> +   but WITHOUT ANY WARRANTY; without even the implied warranty of

> +   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

> +   Lesser General Public License for more details.

> +

> +   You should have received a copy of the GNU Lesser General Public

> +   License along with the GNU C Library; if not, see

> +   <https://www.gnu.org/licenses/>.  */

>  

>  #include <math.h>

>  #include <math_private.h>

>  #include <libm-alias-finite.h>

>  

>  float

> -__ieee754_hypotf(float x, float y)

> +__ieee754_hypotf (float x, float y)

>  {

> -	double d_x, d_y;

> -	int32_t ha, hb;

> -

> -	GET_FLOAT_WORD(ha,x);

> -	ha &= 0x7fffffff;

> -	GET_FLOAT_WORD(hb,y);

> -	hb &= 0x7fffffff;

> -	if (ha == 0x7f800000 && !issignaling (y))

> -	  return fabsf(x);

> -	else if (hb == 0x7f800000 && !issignaling (x))

> -	  return fabsf(y);

> -	else if (ha > 0x7f800000 || hb > 0x7f800000)

> -	  return fabsf(x) * fabsf(y);

> -	else if (ha == 0)

> -	  return fabsf(y);

> -	else if (hb == 0)

> -	  return fabsf(x);

> -

> -	d_x = (double) x;

> -	d_y = (double) y;

> +  if ((isinf (x) || isinf (y))

> +      && !issignaling (x) && !issignaling (y))

> +    return INFINITY;

> +  if (isnan (x) || isnan (y))

> +    return x + y;


why not test NaN before +/-Inf, to avoid the issignaling test?

> -	return (float) sqrt(d_x * d_x + d_y * d_y);

> +  return sqrt ((double) x * (double) x + (double) y * (double) y);


a double-rounding problem can happen here, but anyway the function is not
correctly rounded

>  }

>  #ifndef __ieee754_hypotf

>  libm_alias_finite (__ieee754_hypotf, __hypotf)

> -- 

> 2.30.2

> 

> 


Paul
Carlos O'Donell via Libc-alpha Oct. 7, 2021, 11:37 a.m. | #2
On 07/10/2021 06:44, Paul Zimmermann wrote:
>        Dear Adhemerval,

> 

>> Use isnan()/isinf() instead of GET_FLOAT_WORD and interger operations.

> 

> interger -> integer


Ack.

> 

>> There is also no need to check for 0.0.

>>

>> The file Copyright is also change to use  GPL, the implementation was

> 

> change -> changed


Ack.

> 

>> complete change by 7c10fd3515f to use double precision instead of

> 

> complete change -> completely changed ?


Ack.

> 

>> scaling and this change removes all the GET_FLOAT_WORD usage.

>>

>> Checked on x86_64-linux-gnu.

>> ---

>>  sysdeps/ieee754/flt-32/e_hypotf.c | 57 +++++++++++++------------------

>>  1 file changed, 23 insertions(+), 34 deletions(-)

>>

>> diff --git a/sysdeps/ieee754/flt-32/e_hypotf.c b/sysdeps/ieee754/flt-32/e_hypotf.c

>> index e770947dc1..6495a91cd4 100644

>> --- a/sysdeps/ieee754/flt-32/e_hypotf.c

>> +++ b/sysdeps/ieee754/flt-32/e_hypotf.c

>> @@ -1,46 +1,35 @@

>> -/* e_hypotf.c -- float version of e_hypot.c.

>> - */

>> +/* Euclidean distance function.  Float/Binary32 version.

>> +   Copyright (C) 2012-2021 Free Software Foundation, Inc.

>> +   This file is part of the GNU C Library.

>>  

>> -/*

>> - * ====================================================

>> - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.

>> - *

>> - * Developed at SunPro, a Sun Microsystems, Inc. business.

>> - * Permission to use, copy, modify, and distribute this

>> - * software is freely granted, provided that this notice

>> - * is preserved.

>> - * ====================================================

>> - */

>> +   The GNU C Library is free software; you can redistribute it and/or

>> +   modify it under the terms of the GNU Lesser General Public

>> +   License as published by the Free Software Foundation; either

>> +   version 2.1 of the License, or (at your option) any later version.

>> +

>> +   The GNU C Library is distributed in the hope that it will be useful,

>> +   but WITHOUT ANY WARRANTY; without even the implied warranty of

>> +   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

>> +   Lesser General Public License for more details.

>> +

>> +   You should have received a copy of the GNU Lesser General Public

>> +   License along with the GNU C Library; if not, see

>> +   <https://www.gnu.org/licenses/>.  */

>>  

>>  #include <math.h>

>>  #include <math_private.h>

>>  #include <libm-alias-finite.h>

>>  

>>  float

>> -__ieee754_hypotf(float x, float y)

>> +__ieee754_hypotf (float x, float y)

>>  {

>> -	double d_x, d_y;

>> -	int32_t ha, hb;

>> -

>> -	GET_FLOAT_WORD(ha,x);

>> -	ha &= 0x7fffffff;

>> -	GET_FLOAT_WORD(hb,y);

>> -	hb &= 0x7fffffff;

>> -	if (ha == 0x7f800000 && !issignaling (y))

>> -	  return fabsf(x);

>> -	else if (hb == 0x7f800000 && !issignaling (x))

>> -	  return fabsf(y);

>> -	else if (ha > 0x7f800000 || hb > 0x7f800000)

>> -	  return fabsf(x) * fabsf(y);

>> -	else if (ha == 0)

>> -	  return fabsf(y);

>> -	else if (hb == 0)

>> -	  return fabsf(x);

>> -

>> -	d_x = (double) x;

>> -	d_y = (double) y;

>> +  if ((isinf (x) || isinf (y))

>> +      && !issignaling (x) && !issignaling (y))

>> +    return INFINITY;

>> +  if (isnan (x) || isnan (y))

>> +    return x + y;

> 

> why not test NaN before +/-Inf, to avoid the issignaling test?


Because NaN should be returned iff there is a signaling one, for 
instance hypot (qNaN, Inf) should return Inf.  To test NaN before
Inf we would need to:

  if ((isnan (x) || isnan (y))
      && (issignaling (x) || issignaling (y)))
    return x + y;
  if (isinf (x) || isinf (y))
    return INFINITY;

Which is essentially the same.

> 

>> -	return (float) sqrt(d_x * d_x + d_y * d_y);

>> +  return sqrt ((double) x * (double) x + (double) y * (double) y);

> 

> a double-rounding problem can happen here, but anyway the function is not

> correctly rounded


Yeah, but I think this is preexisting issue in the current approach.

> 

>>  }

>>  #ifndef __ieee754_hypotf

>>  libm_alias_finite (__ieee754_hypotf, __hypotf)

>> -- 

>> 2.30.2

>>

>>

> 

> Paul

>
Paul Zimmermann Oct. 7, 2021, 12:08 p.m. | #3
Dear Adhemerval,

> >> +  if ((isinf (x) || isinf (y))

> >> +      && !issignaling (x) && !issignaling (y))

> >> +    return INFINITY;

> >> +  if (isnan (x) || isnan (y))

> >> +    return x + y;

> > 

> > why not test NaN before +/-Inf, to avoid the issignaling test?

> 

> Because NaN should be returned iff there is a signaling one, for 

> instance hypot (qNaN, Inf) should return Inf.


thank you, I forgot that special case.

Paul

Patch

diff --git a/sysdeps/ieee754/flt-32/e_hypotf.c b/sysdeps/ieee754/flt-32/e_hypotf.c
index e770947dc1..6495a91cd4 100644
--- a/sysdeps/ieee754/flt-32/e_hypotf.c
+++ b/sysdeps/ieee754/flt-32/e_hypotf.c
@@ -1,46 +1,35 @@ 
-/* e_hypotf.c -- float version of e_hypot.c.
- */
+/* Euclidean distance function.  Float/Binary32 version.
+   Copyright (C) 2012-2021 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
 
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <https://www.gnu.org/licenses/>.  */
 
 #include <math.h>
 #include <math_private.h>
 #include <libm-alias-finite.h>
 
 float
-__ieee754_hypotf(float x, float y)
+__ieee754_hypotf (float x, float y)
 {
-	double d_x, d_y;
-	int32_t ha, hb;
-
-	GET_FLOAT_WORD(ha,x);
-	ha &= 0x7fffffff;
-	GET_FLOAT_WORD(hb,y);
-	hb &= 0x7fffffff;
-	if (ha == 0x7f800000 && !issignaling (y))
-	  return fabsf(x);
-	else if (hb == 0x7f800000 && !issignaling (x))
-	  return fabsf(y);
-	else if (ha > 0x7f800000 || hb > 0x7f800000)
-	  return fabsf(x) * fabsf(y);
-	else if (ha == 0)
-	  return fabsf(y);
-	else if (hb == 0)
-	  return fabsf(x);
-
-	d_x = (double) x;
-	d_y = (double) y;
+  if ((isinf (x) || isinf (y))
+      && !issignaling (x) && !issignaling (y))
+    return INFINITY;
+  if (isnan (x) || isnan (y))
+    return x + y;
 
-	return (float) sqrt(d_x * d_x + d_y * d_y);
+  return sqrt ((double) x * (double) x + (double) y * (double) y);
 }
 #ifndef __ieee754_hypotf
 libm_alias_finite (__ieee754_hypotf, __hypotf)