sysdeps/ieee754/flt-32/s_sincosf.c

   1 /* Compute sine and cosine of argument.
   2    Copyright (C) 2017 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    <http://www.gnu.org/licenses/>.  */
  18
  19 #include <errno.h>
  20 #include <math.h>
  21 #include <math_private.h>
  22 #include <libm-alias-float.h>
  23 #include "s_sincosf.h"
  24
  25 #ifndef SINCOSF
  26 # define SINCOSF_FUNC __sincosf
  27 #else
  28 # define SINCOSF_FUNC SINCOSF
  29 #endif
  30
  31 void
  32 SINCOSF_FUNC (float x, float *sinx, float *cosx)
  33 {
  34   double cx;
  35   double theta = x;
  36   double abstheta = fabs (theta);
  37   /* If |x|< Pi/4.  */
  38   if (isless (abstheta, M_PI_4))
  39     {
  40       if (abstheta >= 0x1p-5) /* |x| >= 2^-5.  */
  41         {
  42           const double theta2 = theta * theta;
  43           /* Chebyshev polynomial of the form for sin and cos.  */
  44           cx = C3 + theta2 * C4;
  45           cx = C2 + theta2 * cx;
  46           cx = C1 + theta2 * cx;
  47           cx = C0 + theta2 * cx;
  48           cx = 1.0 + theta2 * cx;
  49           *cosx = cx;
  50           cx = S3 + theta2 * S4;
  51           cx = S2 + theta2 * cx;
  52           cx = S1 + theta2 * cx;
  53           cx = S0 + theta2 * cx;
  54           cx = theta + theta * theta2 * cx;
  55           *sinx = cx;
  56         }
  57       else if (abstheta >= 0x1p-27)     /* |x| >= 2^-27.  */
  58         {
  59           /* A simpler Chebyshev approximation is close enough for this range:
  60              for sin: x+x^3*(SS0+x^2*SS1)
  61              for cos: 1.0+x^2*(CC0+x^3*CC1).  */
  62           const double theta2 = theta * theta;
  63           cx = CC0 + theta * theta2 * CC1;
  64           cx = 1.0 + theta2 * cx;
  65           *cosx = cx;
  66           cx = SS0 + theta2 * SS1;
  67           cx = theta + theta * theta2 * cx;
  68           *sinx = cx;
  69         }
  70       else
  71         {
  72           /* Handle some special cases.  */
  73           if (theta)
  74             *sinx = theta - (theta * SMALL);
  75           else
  76             *sinx = theta;
  77           *cosx = 1.0 - abstheta;
  78         }
  79     }
  80   else                          /* |x| >= Pi/4.  */
  81     {
  82       unsigned int signbit = isless (x, 0);
  83       if (isless (abstheta, 9 * M_PI_4))        /* |x| < 9*Pi/4.  */
  84         {
  85           /* There are cases where FE_UPWARD rounding mode can
  86              produce a result of abstheta * inv_PI_4 == 9,
  87              where abstheta < 9pi/4, so the domain for
  88              pio2_table must go to 5 (9 / 2 + 1).  */
  89           unsigned int n = (abstheta * inv_PI_4) + 1;
  90           theta = abstheta - pio2_table[n / 2];
  91           *sinx = reduced_sin (theta, n, signbit);
  92           *cosx = reduced_cos (theta, n);
  93         }
  94       else if (isless (abstheta, INFINITY))
  95         {
  96           if (abstheta < 0x1p+23)     /* |x| < 2^23.  */
  97             {
  98               unsigned int n = ((unsigned int) (abstheta * inv_PI_4)) + 1;
  99               double x = n / 2;
 100               theta = (abstheta - x * PI_2_hi) - x * PI_2_lo;
 101               /* Argument reduction needed.  */
 102               *sinx = reduced_sin (theta, n, signbit);
 103               *cosx = reduced_cos (theta, n);
 104             }
 105           else                  /* |x| >= 2^23.  */
 106             {
 107               x = fabsf (x);
 108               int exponent;
 109               GET_FLOAT_WORD (exponent, x);
 110               exponent
 111                 = (exponent >> FLOAT_EXPONENT_SHIFT) - FLOAT_EXPONENT_BIAS;
 112               exponent += 3;
 113               exponent /= 28;
 114               double a = invpio4_table[exponent] * x;
 115               double b = invpio4_table[exponent + 1] * x;
 116               double c = invpio4_table[exponent + 2] * x;
 117               double d = invpio4_table[exponent + 3] * x;
 118               uint64_t l = a;
 119               l &= ~0x7;
 120               a -= l;
 121               double e = a + b;
 122               l = e;
 123               e = a - l;
 124               if (l & 1)
 125                 {
 126                   e -= 1.0;
 127                   e += b;
 128                   e += c;
 129                   e += d;
 130                   e *= M_PI_4;
 131                   *sinx = reduced_sin (e, l + 1, signbit);
 132                   *cosx = reduced_cos (e, l + 1);
 133                 }
 134               else
 135                 {
 136                   e += b;
 137                   e += c;
 138                   e += d;
 139                   if (e <= 1.0)
 140                     {
 141                       e *= M_PI_4;
 142                       *sinx = reduced_sin (e, l + 1, signbit);
 143                       *cosx = reduced_cos (e, l + 1);
 144                     }
 145                   else
 146                     {
 147                       l++;
 148                       e -= 2.0;
 149                       e *= M_PI_4;
 150                       *sinx = reduced_sin (e, l + 1, signbit);
 151                       *cosx = reduced_cos (e, l + 1);
 152                     }
 153                 }
 154             }
 155         }
 156       else
 157         {
 158           int32_t ix;
 159           /* High word of x.  */
 160           GET_FLOAT_WORD (ix, abstheta);
 161           /* sin/cos(Inf or NaN) is NaN.  */
 162           *sinx = *cosx = x - x;
 163           if (ix == 0x7f800000)
 164             __set_errno (EDOM);
 165         }
 166     }
 167 }
 168
 169 #ifndef SINCOSF
 170 libm_alias_float (__sincos, sincos)
 171 #endif