public inbox for libc-alpha@sourceware.org
 help / color / mirror / Atom feed
* [PATCH 1/9] Auxiliary function for reduction modulo 2*pi.
@ 2021-03-19 15:06 Paul Zimmermann
  2021-03-19 15:06 ` [PATCH 2/9] Fix the inaccuracy of j0f and of y0f [BZ #14469 and #14471] Paul Zimmermann
                   ` (8 more replies)
  0 siblings, 9 replies; 17+ messages in thread
From: Paul Zimmermann @ 2021-03-19 15:06 UTC (permalink / raw)
  To: libc-alpha

---
 sysdeps/ieee754/flt-32/reduce_aux.c | 55 +++++++++++++++++++++++++++++
 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;
+}
-- 
2.30.2


^ permalink raw reply	[flat|nested] 17+ messages in thread
end of thread, other threads:[~2021-03-30 18:09 UTC | newest]

Thread overview: 17+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2021-03-19 15:06 [PATCH 1/9] Auxiliary function for reduction modulo 2*pi Paul Zimmermann
2021-03-19 15:06 ` [PATCH 2/9] Fix the inaccuracy of j0f and of y0f [BZ #14469 and #14471] Paul Zimmermann
2021-03-30 16:59   ` Adhemerval Zanella
2021-03-19 15:06 ` [PATCH 3/9] Fix the inaccuracy of j1f (BZ 14470) and y1f (BZ 14472) Paul Zimmermann
2021-03-30 17:27   ` Adhemerval Zanella
2021-03-30 17:46     ` Paul Zimmermann
2021-03-30 18:09       ` Adhemerval Zanella
2021-03-19 15:06 ` [PATCH 4/9] Added new entries for j0/j1/y0/y1 and regenerate ulps Paul Zimmermann
2021-03-30 16:13   ` Adhemerval Zanella
2021-03-19 15:06 ` [PATCH 5/9] ulps update for sparc (from Adhemerval Zanella) Paul Zimmermann
2021-03-19 15:06 ` [PATCH 6/9] ulps update for s390 " Paul Zimmermann
2021-03-19 15:06 ` [PATCH 7/9] ulps update for aarch64 " Paul Zimmermann
2021-03-19 15:06 ` [PATCH 8/9] ulps update for powerpc " Paul Zimmermann
2021-03-19 15:06 ` [PATCH 9/9] ulps update for x86_64 Paul Zimmermann
2021-03-30 12:51 ` [PATCH 1/9] Auxiliary function for reduction modulo 2*pi Adhemerval Zanella
2021-03-30 17:24   ` Paul Zimmermann
2021-03-30 17:28     ` Adhemerval Zanella

This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).