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

Message ID 20210319150626.1206905-1-Paul.Zimmermann@inria.fr
State New
Headers show
Series
  • [1/9] Auxiliary function for reduction modulo 2*pi.
Related show

Commit Message

Paul Zimmermann March 19, 2021, 3:06 p.m.
---
 sysdeps/ieee754/flt-32/reduce_aux.c | 55 +++++++++++++++++++++++++++++
 1 file changed, 55 insertions(+)
 create mode 100644 sysdeps/ieee754/flt-32/reduce_aux.c

-- 
2.30.2

Comments

Adhemerval Zanella via Libc-alpha March 30, 2021, 12:51 p.m. | #1
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;

> +}

>
Paul Zimmermann March 30, 2021, 5:24 p.m. | #2
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
Adhemerval Zanella via Libc-alpha March 30, 2021, 5:28 p.m. | #3
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?

Patch

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