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