Message ID | 20210319150626.1206905-1-Paul.Zimmermann@inria.fr |
---|---|

State | New |

Headers | show |

Series | [1/9] Auxiliary function for reduction modulo 2*pi. | expand |

On 19/03/2021 12:06, Paul Zimmermann wrote: > --- > sysdeps/ieee754/flt-32/reduce_aux.c | 55 +++++++++++++++++++++++++++++ I think it is better to move it to a header with proper include guards, since it is usual way of defining internal static inline functions (just rename to reduce_aux.h and add #ifndef guards). I also tried to move it to its own TU to see what kind of code size gain it would yield. For aarch64 I am seeing a reduction of 488 bytes, so if you might consider it if this is not a performance-wise routine. LGTM with the header change. Reviewed-by: Adhemerval Zanella <adhemerval.zanella@linaro.org> > 1 file changed, 55 insertions(+) > create mode 100644 sysdeps/ieee754/flt-32/reduce_aux.c > > diff --git a/sysdeps/ieee754/flt-32/reduce_aux.c b/sysdeps/ieee754/flt-32/reduce_aux.c > new file mode 100644 > index 0000000000..412b4d22cb > --- /dev/null > +++ b/sysdeps/ieee754/flt-32/reduce_aux.c > @@ -0,0 +1,55 @@ > +/* Auxiliary routine for the Bessel functions (j0f, y0f, j1f, y1f). > + Copyright (C) 2021 Free Software Foundation, Inc. > + This file is part of the GNU C Library. > + > + 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/>. */ > + > +/* Return h and update n such that: > + Now x - pi/4 - alpha = h + n*pi/2 mod (2*pi). */ > +static inline double > +reduce_aux (float x, int *n, double alpha) > +{ > + double h; > + h = reduce_large (asuint (x), n); > + /* Now |x| = h+n*pi/2 mod 2*pi. */ > + /* Recover sign. */ > + if (x < 0) > + { > + h = -h; > + *n = -*n; > + } > + /* Subtract pi/4. */ > + double piover2 = 0xc.90fdaa22168cp-3; > + if (h >= 0) > + h -= piover2 / 2; > + else > + { > + h += piover2 / 2; > + (*n) --; > + } > + /* Subtract alpha and reduce if needed mod pi/2. */ > + h -= alpha; > + if (h > piover2) > + { > + h -= piover2; > + (*n) ++; > + } > + else if (h < -piover2) > + { > + h += piover2; > + (*n) --; > + } > + return h; > +} >

thank you Adhemerval for your review. I have just submitted a cumulated patch of this series (I spend too much time maintaining such a series of patches when some to need to be updated and some not). The only changes are the following ones. For reduce_aux.c: I have moved it to a file reduce_aux.h with include guards (btw I noticed that s_sincosf.h has not include guards). I have added your patch with xfail entries for powerpc (it is strange I did not hit that since I've tested on powerpc too). Best regards, Paul

On 30/03/2021 14:24, Paul Zimmermann wrote: > thank you Adhemerval for your review. I have just submitted a cumulated patch > of this series (I spend too much time maintaining such a series of patches > when some to need to be updated and some not). > > The only changes are the following ones. > > For reduce_aux.c: I have moved it to a file reduce_aux.h with include guards > (btw I noticed that s_sincosf.h has not include guards). > > I have added your patch with xfail entries for powerpc (it is strange I did > not hit that since I've tested on powerpc too). Right, so this new patch supersedes this whole set?

diff --git a/sysdeps/ieee754/flt-32/reduce_aux.c b/sysdeps/ieee754/flt-32/reduce_aux.c new file mode 100644 index 0000000000..412b4d22cb --- /dev/null +++ b/sysdeps/ieee754/flt-32/reduce_aux.c @@ -0,0 +1,55 @@ +/* Auxiliary routine for the Bessel functions (j0f, y0f, j1f, y1f). + Copyright (C) 2021 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + 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/>. */ + +/* Return h and update n such that: + Now x - pi/4 - alpha = h + n*pi/2 mod (2*pi). */ +static inline double +reduce_aux (float x, int *n, double alpha) +{ + double h; + h = reduce_large (asuint (x), n); + /* Now |x| = h+n*pi/2 mod 2*pi. */ + /* Recover sign. */ + if (x < 0) + { + h = -h; + *n = -*n; + } + /* Subtract pi/4. */ + double piover2 = 0xc.90fdaa22168cp-3; + if (h >= 0) + h -= piover2 / 2; + else + { + h += piover2 / 2; + (*n) --; + } + /* Subtract alpha and reduce if needed mod pi/2. */ + h -= alpha; + if (h > piover2) + { + h -= piover2; + (*n) ++; + } + else if (h < -piover2) + { + h += piover2; + (*n) --; + } + return h; +}