1  /* Compute sine and cosine of argument.
       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 <errno.h>
      20  #include <stdint.h>
      21  #include <math.h>
      22  #include <math-barriers.h>
      23  #include <libm-alias-float.h>
      24  #include "math_config.h"
      25  #include "s_sincosf.h"
      26  
      27  #ifndef SECTION
      28  # define SECTION
      29  #endif
      30  
      31  #ifndef SINCOSF
      32  # define SINCOSF_FUNC __sincosf
      33  #else
      34  # define SINCOSF_FUNC SINCOSF
      35  #endif
      36  
      37  /* Fast sincosf implementation.  Worst-case ULP is 0.5607, maximum relative
      38     error is 0.5303 * 2^-23.  A single-step range reduction is used for
      39     small values.  Large inputs have their range reduced using fast integer
      40     arithmetic.  */
      41  void
      42  SECTION
      43  SINCOSF_FUNC (float y, float *sinp, float *cosp)
      44  {
      45    double x = y;
      46    double s;
      47    int n;
      48    const sincos_t *p = &__sincosf_table[0];
      49  
      50    if (abstop12 (y) < abstop12 (pio4))
      51      {
      52        double x2 = x * x;
      53  
      54        if (__glibc_unlikely (abstop12 (y) < abstop12 (0x1p-12f)))
      55        {
      56  	/* Force underflow for tiny y.  */
      57  	if (__glibc_unlikely (abstop12 (y) < abstop12 (0x1p-126f)))
      58  	  math_force_eval ((float)x2);
      59  	*sinp = y;
      60  	*cosp = 1.0f;
      61  	return;
      62        }
      63  
      64        sincosf_poly (x, x2, p, 0, sinp, cosp);
      65      }
      66    else if (abstop12 (y) < abstop12 (120.0f))
      67      {
      68        x = reduce_fast (x, p, &n);
      69  
      70        /* Setup the signs for sin and cos.  */
      71        s = p->sign[n & 3];
      72  
      73        if (n & 2)
      74  	p = &__sincosf_table[1];
      75  
      76        sincosf_poly (x * s, x * x, p, n, sinp, cosp);
      77      }
      78    else if (__glibc_likely (abstop12 (y) < abstop12 (INFINITY)))
      79      {
      80        uint32_t xi = asuint (y);
      81        int sign = xi >> 31;
      82  
      83        x = reduce_large (xi, &n);
      84  
      85        /* Setup signs for sin and cos - include original sign.  */
      86        s = p->sign[(n + sign) & 3];
      87  
      88        if ((n + sign) & 2)
      89  	p = &__sincosf_table[1];
      90  
      91        sincosf_poly (x * s, x * x, p, n, sinp, cosp);
      92      }
      93    else
      94      {
      95        /* Return NaN if Inf or NaN for both sin and cos.  */
      96        *sinp = *cosp = y - y;
      97  #if WANT_ERRNO
      98        /* Needed to set errno for +-Inf, the add is a hack to work
      99  	 around a gcc register allocation issue: just passing y
     100  	 affects code generation in the fast path (PR86901).  */
     101        __math_invalidf (y + y);
     102  #endif
     103      }
     104  }
     105  
     106  #ifndef SINCOSF
     107  libm_alias_float (__sincos, sincos)
     108  #endif