sysdeps/ieee754/dbl-64/s_sin.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2024 Free Software Foundation, Inc.
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU Lesser General Public License as published by
   8  * the Free Software Foundation; either version 2.1 of the License, or
   9  * (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU  Lesser General Public License
  17  * along with this program; if not, see <https://www.gnu.org/licenses/>.
  18  */
  19 /****************************************************************************/
  20 /*                                                                          */
  21 /* MODULE_NAME:usncs.c                                                      */
  22 /*                                                                          */
  23 /* FUNCTIONS: usin                                                          */
  24 /*            ucos                                                          */
  25 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h  usncs.h                     */
  26 /*               branred.c sincos.tbl                                       */
  27 /*                                                                          */
  28 /* An ultimate sin and cos routine. Given an IEEE double machine number x   */
  29 /* it computes sin(x) or cos(x) with ~0.55 ULP.                             */
  30 /* Assumption: Machine arithmetic operations are performed in               */
  31 /* round to nearest mode of IEEE 754 standard.                              */
  32 /*                                                                          */
  33 /****************************************************************************/
  34
  35
  36 #include <errno.h>
  37 #include <float.h>
  38 #include "endian.h"
  39 #include "mydefs.h"
  40 #include "usncs.h"
  41 #include <math.h>
  42 #include <math_private.h>
  43 #include <fenv_private.h>
  44 #include <math-underflow.h>
  45 #include <libm-alias-double.h>
  46 #include <fenv.h>
  47
  48 /* Helper macros to compute sin of the input values.  */
  49 #define POLYNOMIAL2(xx) ((((s5 * (xx) + s4) * (xx) + s3) * (xx) + s2) * (xx))
  50
  51 #define POLYNOMIAL(xx) (POLYNOMIAL2 (xx) + s1)
  52
  53 /* The computed polynomial is a variation of the Taylor series expansion for
  54    sin(x):
  55
  56    x - x^3/3! + x^5/5! - x^7/7! + x^9/9! - dx*x^2/2 + dx
  57
  58    The constants s1, s2, s3, etc. are pre-computed values of 1/3!, 1/5! and so
  59    on.  The result is returned to LHS.  */
  60 #define TAYLOR_SIN(xx, x, dx) \
  61 ({                                                                            \
  62   double t = ((POLYNOMIAL (xx)  * (x) - 0.5 * (dx))  * (xx) + (dx));          \
  63   double res = (x) + t;                                                       \
  64   res;                                                                        \
  65 })
  66
  67 #define SINCOS_TABLE_LOOKUP(u, sn, ssn, cs, ccs) \
  68 ({                                                                            \
  69   int4 k = u.i[LOW_HALF] << 2;                                                \
  70   sn = __sincostab.x[k];                                                      \
  71   ssn = __sincostab.x[k + 1];                                                 \
  72   cs = __sincostab.x[k + 2];                                                  \
  73   ccs = __sincostab.x[k + 3];                                                 \
  74 })
  75
  76 #ifndef SECTION
  77 # define SECTION
  78 #endif
  79
  80 extern const union
  81 {
  82   int4 i[880];
  83   double x[440];
  84 } __sincostab attribute_hidden;
  85
  86 static const double
  87   sn3 = -1.66666666666664880952546298448555E-01,
  88   sn5 = 8.33333214285722277379541354343671E-03,
  89   cs2 = 4.99999999999999999999950396842453E-01,
  90   cs4 = -4.16666666666664434524222570944589E-02,
  91   cs6 = 1.38888874007937613028114285595617E-03;
  92
  93 int __branred (double x, double *a, double *aa);
  94
  95 /* Given a number partitioned into X and DX, this function computes the cosine
  96    of the number by combining the sin and cos of X (as computed by a variation
  97    of the Taylor series) with the values looked up from the sin/cos table to
  98    get the result.  */
  99 static __always_inline double
 100 do_cos (double x, double dx)
 101 {
 102   mynumber u;
 103
 104   if (x < 0)
 105     dx = -dx;
 106
 107   u.x = big + fabs (x);
 108   x = fabs (x) - (u.x - big) + dx;
 109
 110   double xx, s, sn, ssn, c, cs, ccs, cor;
 111   xx = x * x;
 112   s = x + x * xx * (sn3 + xx * sn5);
 113   c = xx * (cs2 + xx * (cs4 + xx * cs6));
 114   SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
 115   cor = (ccs - s * ssn - cs * c) - sn * s;
 116   return cs + cor;
 117 }
 118
 119 /* Given a number partitioned into X and DX, this function computes the sine of
 120    the number by combining the sin and cos of X (as computed by a variation of
 121    the Taylor series) with the values looked up from the sin/cos table to get
 122    the result.  */
 123 static __always_inline double
 124 do_sin (double x, double dx)
 125 {
 126   double xold = x;
 127   /* Max ULP is 0.501 if |x| < 0.126, otherwise ULP is 0.518.  */
 128   if (fabs (x) < 0.126)
 129     return TAYLOR_SIN (x * x, x, dx);
 130
 131   mynumber u;
 132
 133   if (x <= 0)
 134     dx = -dx;
 135   u.x = big + fabs (x);
 136   x = fabs (x) - (u.x - big);
 137
 138   double xx, s, sn, ssn, c, cs, ccs, cor;
 139   xx = x * x;
 140   s = x + (dx + x * xx * (sn3 + xx * sn5));
 141   c = x * dx + xx * (cs2 + xx * (cs4 + xx * cs6));
 142   SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
 143   cor = (ssn + s * ccs - sn * c) + cs * s;
 144   return copysign (sn + cor, xold);
 145 }
 146
 147 /* Reduce range of x to within PI/2 with abs (x) < 105414350.  The high part
 148    is written to *a, the low part to *da.  Range reduction is accurate to 136
 149    bits so that when x is large and *a very close to zero, all 53 bits of *a
 150    are correct.  */
 151 static __always_inline int4
 152 reduce_sincos (double x, double *a, double *da)
 153 {
 154   mynumber v;
 155
 156   double t = (x * hpinv + toint);
 157   double xn = t - toint;
 158   v.x = t;
 159   double y = (x - xn * mp1) - xn * mp2;
 160   int4 n = v.i[LOW_HALF] & 3;
 161
 162   double b, db, t1, t2;
 163   t1 = xn * pp3;
 164   t2 = y - t1;
 165   db = (y - t2) - t1;
 166
 167   t1 = xn * pp4;
 168   b = t2 - t1;
 169   db += (t2 - b) - t1;
 170
 171   *a = b;
 172   *da = db;
 173   return n;
 174 }
 175
 176 /* Compute sin or cos (A + DA) for the given quadrant N.  */
 177 static __always_inline double
 178 do_sincos (double a, double da, int4 n)
 179 {
 180   double retval;
 181
 182   if (n & 1)
 183     /* Max ULP is 0.513.  */
 184     retval = do_cos (a, da);
 185   else
 186     /* Max ULP is 0.501 if xx < 0.01588, otherwise ULP is 0.518.  */
 187     retval = do_sin (a, da);
 188
 189   return (n & 2) ? -retval : retval;
 190 }
 191
 192
 193 /*******************************************************************/
 194 /* An ultimate sin routine. Given an IEEE double machine number x  */
 195 /* it computes the rounded value of sin(x).                        */
 196 /*******************************************************************/
 197 #ifndef IN_SINCOS
 198 double
 199 SECTION
 200 __sin (double x)
 201 {
 202   double t, a, da;
 203   mynumber u;
 204   int4 k, m, n;
 205   double retval = 0;
 206
 207   SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
 208
 209   u.x = x;
 210   m = u.i[HIGH_HALF];
 211   k = 0x7fffffff & m;           /* no sign           */
 212   if (k < 0x3e500000)           /* if x->0 =>sin(x)=x */
 213     {
 214       math_check_force_underflow (x);
 215       retval = x;
 216     }
 217 /*--------------------------- 2^-26<|x|< 0.855469---------------------- */
 218   else if (k < 0x3feb6000)
 219     {
 220       /* Max ULP is 0.548.  */
 221       retval = do_sin (x, 0);
 222     }                           /*   else  if (k < 0x3feb6000)    */
 223
 224 /*----------------------- 0.855469  <|x|<2.426265  ----------------------*/
 225   else if (k < 0x400368fd)
 226     {
 227       t = hp0 - fabs (x);
 228       /* Max ULP is 0.51.  */
 229       retval = copysign (do_cos (t, hp1), x);
 230     }                           /*   else  if (k < 0x400368fd)    */
 231
 232 /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/
 233   else if (k < 0x419921FB)
 234     {
 235       n = reduce_sincos (x, &a, &da);
 236       retval = do_sincos (a, da, n);
 237     }                           /*   else  if (k <  0x419921FB )    */
 238
 239 /* --------------------105414350 <|x| <2^1024------------------------------*/
 240   else if (k < 0x7ff00000)
 241     {
 242       n = __branred (x, &a, &da);
 243       retval = do_sincos (a, da, n);
 244     }
 245 /*--------------------- |x| > 2^1024 ----------------------------------*/
 246   else
 247     {
 248       if (k == 0x7ff00000 && u.i[LOW_HALF] == 0)
 249         __set_errno (EDOM);
 250       retval = x / x;
 251     }
 252
 253   return retval;
 254 }
 255
 256
 257 /*******************************************************************/
 258 /* An ultimate cos routine. Given an IEEE double machine number x  */
 259 /* it computes the rounded value of cos(x).                        */
 260 /*******************************************************************/
 261
 262 double
 263 SECTION
 264 __cos (double x)
 265 {
 266   double y, a, da;
 267   mynumber u;
 268   int4 k, m, n;
 269
 270   double retval = 0;
 271
 272   SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
 273
 274   u.x = x;
 275   m = u.i[HIGH_HALF];
 276   k = 0x7fffffff & m;
 277
 278   /* |x|<2^-27 => cos(x)=1 */
 279   if (k < 0x3e400000)
 280     retval = 1.0;
 281
 282   else if (k < 0x3feb6000)
 283     {                           /* 2^-27 < |x| < 0.855469 */
 284       /* Max ULP is 0.51.  */
 285       retval = do_cos (x, 0);
 286     }                           /*   else  if (k < 0x3feb6000)    */
 287
 288   else if (k < 0x400368fd)
 289     { /* 0.855469  <|x|<2.426265  */ ;
 290       y = hp0 - fabs (x);
 291       a = y + hp1;
 292       da = (y - a) + hp1;
 293       /* Max ULP is 0.501 if xx < 0.01588 or 0.518 otherwise.
 294          Range reduction uses 106 bits here which is sufficient.  */
 295       retval = do_sin (a, da);
 296     }                           /*   else  if (k < 0x400368fd)    */
 297
 298   else if (k < 0x419921FB)
 299     {                           /* 2.426265<|x|< 105414350 */
 300       n = reduce_sincos (x, &a, &da);
 301       retval = do_sincos (a, da, n + 1);
 302     }                           /*   else  if (k <  0x419921FB )    */
 303
 304   /* 105414350 <|x| <2^1024 */
 305   else if (k < 0x7ff00000)
 306     {
 307       n = __branred (x, &a, &da);
 308       retval = do_sincos (a, da, n + 1);
 309     }
 310
 311   else
 312     {
 313       if (k == 0x7ff00000 && u.i[LOW_HALF] == 0)
 314         __set_errno (EDOM);
 315       retval = x / x;           /* |x| > 2^1024 */
 316     }
 317
 318   return retval;
 319 }
 320
 321 #ifndef __cos
 322 libm_alias_double (__cos, cos)
 323 #endif
 324 #ifndef __sin
 325 libm_alias_double (__sin, sin)
 326 #endif
 327
 328 #endif