sysdeps/ieee754/flt-32/s_sinf.c

   1 /* Compute sine 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
  24 #ifndef SINF
  25 # define SINF_FUNC __sinf
  26 #else
  27 # define SINF_FUNC SINF
  28 #endif
  29
  30 /* Chebyshev constants for cos, range -PI/4 - PI/4.  */
  31 static const double C0 = -0x1.ffffffffe98aep-2;
  32 static const double C1 =  0x1.55555545c50c7p-5;
  33 static const double C2 = -0x1.6c16b348b6874p-10;
  34 static const double C3 =  0x1.a00eb9ac43ccp-16;
  35 static const double C4 = -0x1.23c97dd8844d7p-22;
  36
  37 /* Chebyshev constants for sin, range -PI/4 - PI/4.  */
  38 static const double S0 = -0x1.5555555551cd9p-3;
  39 static const double S1 =  0x1.1111110c2688bp-7;
  40 static const double S2 = -0x1.a019f8b4bd1f9p-13;
  41 static const double S3 =  0x1.71d7264e6b5b4p-19;
  42 static const double S4 = -0x1.a947e1674b58ap-26;
  43
  44 /* Chebyshev constants for sin, range 2^-27 - 2^-5.  */
  45 static const double SS0 = -0x1.555555543d49dp-3;
  46 static const double SS1 =  0x1.110f475cec8c5p-7;
  47
  48 /* PI/2 with 98 bits of accuracy.  */
  49 static const double PI_2_hi = -0x1.921fb544p+0;
  50 static const double PI_2_lo = -0x1.0b4611a626332p-34;
  51
  52 static const double SMALL = 0x1p-50; /* 2^-50.  */
  53 static const double inv_PI_4 = 0x1.45f306dc9c883p+0; /* 4/PI.  */
  54
  55 #define FLOAT_EXPONENT_SHIFT 23
  56 #define FLOAT_EXPONENT_BIAS 127
  57
  58 static const double pio2_table[] = {
  59   0 * M_PI_2,
  60   1 * M_PI_2,
  61   2 * M_PI_2,
  62   3 * M_PI_2,
  63   4 * M_PI_2,
  64   5 * M_PI_2
  65 };
  66
  67 static const double invpio4_table[] = {
  68   0x0p+0,
  69   0x1.45f306cp+0,
  70   0x1.c9c882ap-28,
  71   0x1.4fe13a8p-58,
  72   0x1.f47d4dp-85,
  73   0x1.bb81b6cp-112,
  74   0x1.4acc9ep-142,
  75   0x1.0e4107cp-169
  76 };
  77
  78 static const double ones[] = { 1.0, -1.0 };
  79
  80 /* Compute the sine value using Chebyshev polynomials where
  81    THETA is the range reduced absolute value of the input
  82    and it is less than Pi/4,
  83    N is calculated as trunc(|x|/(Pi/4)) + 1 and it is used to decide
  84    whether a sine or cosine approximation is more accurate and
  85    SIGNBIT is used to add the correct sign after the Chebyshev
  86    polynomial is computed.  */
  87 static inline float
  88 reduced (const double theta, const unsigned int n,
  89          const unsigned int signbit)
  90 {
  91   double sx;
  92   const double theta2 = theta * theta;
  93   /* We are operating on |x|, so we need to add back the original
  94      signbit for sinf.  */
  95   double sign;
  96   /* Determine positive or negative primary interval.  */
  97   sign = ones[((n >> 2) & 1) ^ signbit];
  98   /* Are we in the primary interval of sin or cos?  */
  99   if ((n & 2) == 0)
 100     {
 101       /* Here sinf() is calculated using sin Chebyshev polynomial:
 102         x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))).  */
 103       sx = S3 + theta2 * S4;     /* S3+x^2*S4.  */
 104       sx = S2 + theta2 * sx;     /* S2+x^2*(S3+x^2*S4).  */
 105       sx = S1 + theta2 * sx;     /* S1+x^2*(S2+x^2*(S3+x^2*S4)).  */
 106       sx = S0 + theta2 * sx;     /* S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4))).  */
 107       sx = theta + theta * theta2 * sx;
 108     }
 109   else
 110     {
 111      /* Here sinf() is calculated using cos Chebyshev polynomial:
 112         1.0+x^2*(C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4)))).  */
 113       sx = C3 + theta2 * C4;     /* C3+x^2*C4.  */
 114       sx = C2 + theta2 * sx;     /* C2+x^2*(C3+x^2*C4).  */
 115       sx = C1 + theta2 * sx;     /* C1+x^2*(C2+x^2*(C3+x^2*C4)).  */
 116       sx = C0 + theta2 * sx;     /* C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4))).  */
 117       sx = 1.0 + theta2 * sx;
 118     }
 119
 120   /* Add in the signbit and assign the result.  */
 121   return sign * sx;
 122 }
 123
 124 float
 125 SINF_FUNC (float x)
 126 {
 127   double cx;
 128   double theta = x;
 129   double abstheta = fabs (theta);
 130   /* If |x|< Pi/4.  */
 131   if (isless (abstheta, M_PI_4))
 132     {
 133       if (abstheta >= 0x1p-5) /* |x| >= 2^-5.  */
 134         {
 135           const double theta2 = theta * theta;
 136           /* Chebyshev polynomial of the form for sin
 137              x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))).  */
 138           cx = S3 + theta2 * S4;
 139           cx = S2 + theta2 * cx;
 140           cx = S1 + theta2 * cx;
 141           cx = S0 + theta2 * cx;
 142           cx = theta + theta * theta2 * cx;
 143           return cx;
 144         }
 145       else if (abstheta >= 0x1p-27)     /* |x| >= 2^-27.  */
 146         {
 147           /* A simpler Chebyshev approximation is close enough for this range:
 148              for sin: x+x^3*(SS0+x^2*SS1).  */
 149           const double theta2 = theta * theta;
 150           cx = SS0 + theta2 * SS1;
 151           cx = theta + theta * theta2 * cx;
 152           return cx;
 153         }
 154       else
 155         {
 156           /* Handle some special cases.  */
 157           if (theta)
 158             return theta - (theta * SMALL);
 159           else
 160             return theta;
 161         }
 162     }
 163   else                          /* |x| >= Pi/4.  */
 164     {
 165       unsigned int signbit = isless (x, 0);
 166       if (isless (abstheta, 9 * M_PI_4))        /* |x| < 9*Pi/4.  */
 167         {
 168           /* There are cases where FE_UPWARD rounding mode can
 169              produce a result of abstheta * inv_PI_4 == 9,
 170              where abstheta < 9pi/4, so the domain for
 171              pio2_table must go to 5 (9 / 2 + 1).  */
 172           unsigned int n = (abstheta * inv_PI_4) + 1;
 173           theta = abstheta - pio2_table[n / 2];
 174           return reduced (theta, n, signbit);
 175         }
 176       else if (isless (abstheta, INFINITY))
 177         {
 178           if (abstheta < 0x1p+23)     /* |x| < 2^23.  */
 179             {
 180               unsigned int n = ((unsigned int) (abstheta * inv_PI_4)) + 1;
 181               double x = n / 2;
 182               theta = x * PI_2_lo + (x * PI_2_hi + abstheta);
 183               /* Argument reduction needed.  */
 184               return reduced (theta, n, signbit);
 185             }
 186           else                  /* |x| >= 2^23.  */
 187             {
 188               x = fabsf (x);
 189               int exponent;
 190               GET_FLOAT_WORD (exponent, x);
 191               exponent
 192                 = (exponent >> FLOAT_EXPONENT_SHIFT) - FLOAT_EXPONENT_BIAS;
 193               exponent += 3;
 194               exponent /= 28;
 195               double a = invpio4_table[exponent] * x;
 196               double b = invpio4_table[exponent + 1] * x;
 197               double c = invpio4_table[exponent + 2] * x;
 198               double d = invpio4_table[exponent + 3] * x;
 199               uint64_t l = a;
 200               l &= ~0x7;
 201               a -= l;
 202               double e = a + b;
 203               l = e;
 204               e = a - l;
 205               if (l & 1)
 206                 {
 207                   e -= 1.0;
 208                   e += b;
 209                   e += c;
 210                   e += d;
 211                   e *= M_PI_4;
 212                   return reduced (e, l + 1, signbit);
 213                 }
 214               else
 215                 {
 216                   e += b;
 217                   e += c;
 218                   e += d;
 219                   if (e <= 1.0)
 220                     {
 221                       e *= M_PI_4;
 222                       return reduced (e, l + 1, signbit);
 223                     }
 224                   else
 225                     {
 226                       l++;
 227                       e -= 2.0;
 228                       e *= M_PI_4;
 229                       return reduced (e, l + 1, signbit);
 230                     }
 231                 }
 232             }
 233         }
 234       else
 235         {
 236           int32_t ix;
 237           /* High word of x.  */
 238           GET_FLOAT_WORD (ix, abstheta);
 239           /* Sin(Inf or NaN) is NaN.  */
 240           if (ix == 0x7f800000)
 241             __set_errno (EDOM);
 242           return x - x;
 243         }
 244     }
 245 }
 246
 247 #ifndef SINF
 248 libm_alias_float (__sin, sin)
 249 #endif