1  /* Used by sinf, cosf and sincosf functions.
       2     Copyright (C) 2018-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  #include <stdint.h>
      20  #include <math.h>
      21  #include "math_config.h"
      22  #include <sincosf_poly.h>
      23  
      24  /* 2PI * 2^-64.  */
      25  static const double pi63 = 0x1.921FB54442D18p-62;
      26  /* PI / 4.  */
      27  static const float pio4 = 0x1.921FB6p-1f;
      28  
      29  /* Polynomial data (the cosine polynomial is negated in the 2nd entry).  */
      30  extern const sincos_t __sincosf_table[2] attribute_hidden;
      31  
      32  /* Table with 4/PI to 192 bit precision.  */
      33  extern const uint32_t __inv_pio4[] attribute_hidden;
      34  
      35  /* Top 12 bits of the float representation with the sign bit cleared.  */
      36  static inline uint32_t
      37  abstop12 (float x)
      38  {
      39    return (asuint (x) >> 20) & 0x7ff;
      40  }
      41  
      42  /* Fast range reduction using single multiply-subtract.  Return the modulo of
      43     X as a value between -PI/4 and PI/4 and store the quadrant in NP.
      44     The values for PI/2 and 2/PI are accessed via P.  Since PI/2 as a double
      45     is accurate to 55 bits and the worst-case cancellation happens at 6 * PI/4,
      46     the result is accurate for |X| <= 120.0.  */
      47  static inline double
      48  reduce_fast (double x, const sincos_t *p, int *np)
      49  {
      50    double r;
      51  #if TOINT_INTRINSICS
      52    /* Use fast round and lround instructions when available.  */
      53    r = x * p->hpi_inv;
      54    *np = converttoint (r);
      55    return x - roundtoint (r) * p->hpi;
      56  #else
      57    /* Use scaled float to int conversion with explicit rounding.
      58       hpi_inv is prescaled by 2^24 so the quadrant ends up in bits 24..31.
      59       This avoids inaccuracies introduced by truncating negative values.  */
      60    r = x * p->hpi_inv;
      61    int n = ((int32_t)r + 0x800000) >> 24;
      62    *np = n;
      63    return x - n * p->hpi;
      64  #endif
      65  }
      66  
      67  /* Reduce the range of XI to a multiple of PI/2 using fast integer arithmetic.
      68     XI is a reinterpreted float and must be >= 2.0f (the sign bit is ignored).
      69     Return the modulo between -PI/4 and PI/4 and store the quadrant in NP.
      70     Reduction uses a table of 4/PI with 192 bits of precision.  A 32x96->128 bit
      71     multiply computes the exact 2.62-bit fixed-point modulo.  Since the result
      72     can have at most 29 leading zeros after the binary point, the double
      73     precision result is accurate to 33 bits.  */
      74  static inline double
      75  reduce_large (uint32_t xi, int *np)
      76  {
      77    const uint32_t *arr = &__inv_pio4[(xi >> 26) & 15];
      78    int shift = (xi >> 23) & 7;
      79    uint64_t n, res0, res1, res2;
      80  
      81    xi = (xi & 0xffffff) | 0x800000;
      82    xi <<= shift;
      83  
      84    res0 = xi * arr[0];
      85    res1 = (uint64_t)xi * arr[4];
      86    res2 = (uint64_t)xi * arr[8];
      87    res0 = (res2 >> 32) | (res0 << 32);
      88    res0 += res1;
      89  
      90    n = (res0 + (1ULL << 61)) >> 62;
      91    res0 -= n << 62;
      92    double x = (int64_t)res0;
      93    *np = n;
      94    return x * pi63;
      95  }