Re: [PATCH 1/9] Auxiliary function for reduction modulo 2*pi.

[Adhemerval Zanella <adhemerval.zanella@linaro.org>]

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;
> +}
>