sysdeps/ieee754/flt-32/reduce_aux.h

   1 /* Auxiliary routine for the Bessel functions (j0f, y0f, j1f, y1f).
   2    Copyright (C) 2021-2023 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4
   5    The GNU C Library is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU Lesser General Public
   7    License as published by the Free Software Foundation; either
   8    version 2.1 of the License, or (at your option) any later version.
   9
  10    The GNU C Library is distributed in the hope that it will be useful,
  11    but WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13    Lesser General Public License for more details.
  14
  15    You should have received a copy of the GNU Lesser General Public
  16    License along with the GNU C Library; if not, see
  17    <https://www.gnu.org/licenses/>.  */
  18
  19 #ifndef _MATH_REDUCE_AUX_H
  20 #define _MATH_REDUCE_AUX_H
  21
  22 #include <math.h>
  23 #include <math_private.h>
  24 #include <s_sincosf.h>
  25
  26 /* Return h and update n such that:
  27    Now x - pi/4 - alpha = h + n*pi/2 mod (2*pi).  */
  28 static inline double
  29 reduce_aux (float x, int *n, double alpha)
  30 {
  31   double h;
  32   h = reduce_large (asuint (x), n);
  33   /* Now |x| = h+n*pi/2 mod 2*pi.  */
  34   /* Recover sign.  */
  35   if (x < 0)
  36     {
  37       h = -h;
  38       *n = -*n;
  39     }
  40   /* Subtract pi/4.  */
  41   double piover2 = 0xc.90fdaa22168cp-3;
  42   if (h >= 0)
  43     h -= piover2 / 2;
  44   else
  45     {
  46       h += piover2 / 2;
  47       (*n) --;
  48     }
  49   /* Subtract alpha and reduce if needed mod pi/2.  */
  50   h -= alpha;
  51   if (h > piover2)
  52     {
  53       h -= piover2;
  54       (*n) ++;
  55     }
  56   else if (h < -piover2)
  57     {
  58       h += piover2;
  59       (*n) --;
  60     }
  61   return h;
  62 }
  63
  64 #endif