sysdeps/ieee754/flt-32/s_sincosf.c

   1 /* Compute sine and cosine of argument.
   2    Copyright (C) 2018-2024 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