Message ID | 20211006180557.933826-4-adhemerval.zanella@linaro.org |
---|---|
State | New |
Headers | show |
Series | Improve hypot() | expand |
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
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 >
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
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)