sysdeps/ieee754/dbl-64/s_sin.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2015 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 <http://www.gnu.org/licenses/>.
  18  */
  19 /****************************************************************************/
  20 /*                                                                          */
  21 /* MODULE_NAME:usncs.c                                                      */
  22 /*                                                                          */
  23 /* FUNCTIONS: usin                                                          */
  24 /*            ucos                                                          */
  25 /*            slow                                                          */
  26 /*            slow1                                                         */
  27 /*            slow2                                                         */
  28 /*            sloww                                                         */
  29 /*            sloww1                                                        */
  30 /*            sloww2                                                        */
  31 /*            bsloww                                                        */
  32 /*            bsloww1                                                       */
  33 /*            bsloww2                                                       */
  34 /*            cslow2                                                        */
  35 /*            csloww                                                        */
  36 /*            csloww1                                                       */
  37 /*            csloww2                                                       */
  38 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h  usncs.h                     */
  39 /*               branred.c sincos32.c dosincos.c mpa.c                      */
  40 /*               sincos.tbl                                                 */
  41 /*                                                                          */
  42 /* An ultimate sin and  routine. Given an IEEE double machine number x       */
  43 /* it computes the correctly rounded (to nearest) value of sin(x) or cos(x) */
  44 /* Assumption: Machine arithmetic operations are performed in               */
  45 /* round to nearest mode of IEEE 754 standard.                              */
  46 /*                                                                          */
  47 /****************************************************************************/
  48
  49
  50 #include <errno.h>
  51 #include "endian.h"
  52 #include "mydefs.h"
  53 #include "usncs.h"
  54 #include "MathLib.h"
  55 #include <math.h>
  56 #include <math_private.h>
  57 #include <fenv.h>
  58
  59 /* Helper macros to compute sin of the input values.  */
  60 #define POLYNOMIAL2(xx) ((((s5 * (xx) + s4) * (xx) + s3) * (xx) + s2) * (xx))
  61
  62 #define POLYNOMIAL(xx) (POLYNOMIAL2 (xx) + s1)
  63
  64 /* The computed polynomial is a variation of the Taylor series expansion for
  65    sin(a):
  66
  67    a - a^3/3! + a^5/5! - a^7/7! + a^9/9! + (1 - a^2) * da / 2
  68
  69    The constants s1, s2, s3, etc. are pre-computed values of 1/3!, 1/5! and so
  70    on.  The result is returned to LHS and correction in COR.  */
  71 #define TAYLOR_SIN(xx, a, da, cor) \
  72 ({                                                                            \
  73   double t = ((POLYNOMIAL (xx)  * (a) - 0.5 * (da))  * (xx) + (da));          \
  74   double res = (a) + t;                                                       \
  75   (cor) = ((a) - res) + t;                                                    \
  76   res;                                                                        \
  77 })
  78
  79 /* This is again a variation of the Taylor series expansion with the term
  80    x^3/3! expanded into the following for better accuracy:
  81
  82    bb * x ^ 3 + 3 * aa * x * x1 * x2 + aa * x1 ^ 3 + aa * x2 ^ 3
  83
  84    The correction term is dx and bb + aa = -1/3!
  85    */
  86 #define TAYLOR_SLOW(x0, dx, cor) \
  87 ({                                                                            \
  88   static const double th2_36 = 206158430208.0;  /*    1.5*2**37   */          \
  89   double xx = (x0) * (x0);                                                    \
  90   double x1 = ((x0) + th2_36) - th2_36;                                       \
  91   double y = aa * x1 * x1 * x1;                                               \
  92   double r = (x0) + y;                                                        \
  93   double x2 = ((x0) - x1) + (dx);                                             \
  94   double t = (((POLYNOMIAL2 (xx) + bb) * xx + 3.0 * aa * x1 * x2)             \
  95               * (x0)  + aa * x2 * x2 * x2 + (dx));                            \
  96   t = (((x0) - r) + y) + t;                                                   \
  97   double res = r + t;                                                         \
  98   (cor) = (r - res) + t;                                                      \
  99   res;                                                                        \
 100 })
 101
 102 #define SINCOS_TABLE_LOOKUP(u, sn, ssn, cs, ccs) \
 103 ({                                                                            \
 104   int4 k = u.i[LOW_HALF] << 2;                                                \
 105   sn = __sincostab.x[k];                                                      \
 106   ssn = __sincostab.x[k + 1];                                                 \
 107   cs = __sincostab.x[k + 2];                                                  \
 108   ccs = __sincostab.x[k + 3];                                                 \
 109 })
 110
 111 #ifndef SECTION
 112 # define SECTION
 113 #endif
 114
 115 extern const union
 116 {
 117   int4 i[880];
 118   double x[440];
 119 } __sincostab attribute_hidden;
 120
 121 static const double
 122   sn3 = -1.66666666666664880952546298448555E-01,
 123   sn5 = 8.33333214285722277379541354343671E-03,
 124   cs2 = 4.99999999999999999999950396842453E-01,
 125   cs4 = -4.16666666666664434524222570944589E-02,
 126   cs6 = 1.38888874007937613028114285595617E-03;
 127
 128 static const double t22 = 0x1.8p22;
 129
 130 void __dubsin (double x, double dx, double w[]);
 131 void __docos (double x, double dx, double w[]);
 132 double __mpsin (double x, double dx, bool reduce_range);
 133 double __mpcos (double x, double dx, bool reduce_range);
 134 static double slow (double x);
 135 static double slow1 (double x);
 136 static double slow2 (double x);
 137 static double sloww (double x, double dx, double orig);
 138 static double sloww1 (double x, double dx, double orig, int m);
 139 static double sloww2 (double x, double dx, double orig, int n);
 140 static double bsloww (double x, double dx, double orig, int n);
 141 static double bsloww1 (double x, double dx, double orig, int n);
 142 static double bsloww2 (double x, double dx, double orig, int n);
 143 int __branred (double x, double *a, double *aa);
 144 static double cslow2 (double x);
 145 static double csloww (double x, double dx, double orig);
 146 static double csloww1 (double x, double dx, double orig, int m);
 147 static double csloww2 (double x, double dx, double orig, int n);
 148
 149 /* Given a number partitioned into U and X such that U is an index into the
 150    sin/cos table, this macro computes the cosine of the number by combining
 151    the sin and cos of X (as computed by a variation of the Taylor series) with
 152    the values looked up from the sin/cos table to get the result in RES and a
 153    correction value in COR.  */
 154 static double
 155 do_cos (mynumber u, double x, double *corp)
 156 {
 157   double xx, s, sn, ssn, c, cs, ccs, res, cor;
 158   xx = x * x;
 159   s = x + x * xx * (sn3 + xx * sn5);
 160   c = xx * (cs2 + xx * (cs4 + xx * cs6));
 161   SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
 162   cor = (ccs - s * ssn - cs * c) - sn * s;
 163   res = cs + cor;
 164   cor = (cs - res) + cor;
 165   *corp = cor;
 166   return res;
 167 }
 168
 169 /* A more precise variant of DO_COS where the number is partitioned into U, X
 170    and DX.  EPS is the adjustment to the correction COR.  */
 171 static double
 172 do_cos_slow (mynumber u, double x, double dx, double eps, double *corp)
 173 {
 174   double xx, y, x1, x2, e1, e2, res, cor;
 175   double s, sn, ssn, c, cs, ccs;
 176   xx = x * x;
 177   s = x * xx * (sn3 + xx * sn5);
 178   c = x * dx + xx * (cs2 + xx * (cs4 + xx * cs6));
 179   SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
 180   x1 = (x + t22) - t22;
 181   x2 = (x - x1) + dx;
 182   e1 = (sn + t22) - t22;
 183   e2 = (sn - e1) + ssn;
 184   cor = (ccs - cs * c - e1 * x2 - e2 * x) - sn * s;
 185   y = cs - e1 * x1;
 186   cor = cor + ((cs - y) - e1 * x1);
 187   res = y + cor;
 188   cor = (y - res) + cor;
 189   if (cor > 0)
 190     cor = 1.0005 * cor + eps;
 191   else
 192     cor = 1.0005 * cor - eps;
 193   *corp = cor;
 194   return res;
 195 }
 196
 197 /* Given a number partitioned into U and X and DX such that U is an index into
 198    the sin/cos table, this macro computes the sine of the number by combining
 199    the sin and cos of X (as computed by a variation of the Taylor series) with
 200    the values looked up from the sin/cos table to get the result in RES and a
 201    correction value in COR.  */
 202 static double
 203 do_sin (mynumber u, double x, double dx, double *corp)
 204 {
 205   double xx, s, sn, ssn, c, cs, ccs, cor, res;
 206   xx = x * x;
 207   s = x + (dx + x * xx * (sn3 + xx * sn5));
 208   c = x * dx + xx * (cs2 + xx * (cs4 + xx * cs6));
 209   SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
 210   cor = (ssn + s * ccs - sn * c) + cs * s;
 211   res = sn + cor;
 212   cor = (sn - res) + cor;
 213   *corp = cor;
 214   return res;
 215 }
 216
 217 /* A more precise variant of res = do_sin where the number is partitioned into U, X
 218    and DX.  EPS is the adjustment to the correction COR.  */
 219 static double
 220 do_sin_slow (mynumber u, double x, double dx, double eps, double *corp)
 221 {
 222   double xx, y, x1, x2, c1, c2, res, cor;
 223   double s, sn, ssn, c, cs, ccs;
 224   xx = x * x;
 225   s = x * xx * (sn3 + xx * sn5);
 226   c = xx * (cs2 + xx * (cs4 + xx * cs6));
 227   SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
 228   x1 = (x + t22) - t22;
 229   x2 = (x - x1) + dx;
 230   c1 = (cs + t22) - t22;
 231   c2 = (cs - c1) + ccs;
 232   cor = (ssn + s * ccs + cs * s + c2 * x + c1 * x2 - sn * x * dx) - sn * c;
 233   y = sn + c1 * x1;
 234   cor = cor + ((sn - y) + c1 * x1);
 235   res = y + cor;
 236   cor = (y - res) + cor;
 237   if (cor > 0)
 238     cor = 1.0005 * cor + eps;
 239   else
 240     cor = 1.0005 * cor - eps;
 241   *corp = cor;
 242   return res;
 243 }
 244
 245 /* Reduce range of X and compute sin of a + da.  K is the amount by which to
 246    rotate the quadrants.  This allows us to use the same routine to compute cos
 247    by simply rotating the quadrants by 1.  */
 248 static inline double
 249 __always_inline
 250 reduce_and_compute (double x, unsigned int k)
 251 {
 252   double retval = 0, a, da;
 253   unsigned int n = __branred (x, &a, &da);
 254   k = (n + k) % 4;
 255   switch (k)
 256     {
 257       case 0:
 258         if (a * a < 0.01588)
 259           retval = bsloww (a, da, x, n);
 260         else
 261           retval = bsloww1 (a, da, x, n);
 262         break;
 263       case 2:
 264         if (a * a < 0.01588)
 265           retval = bsloww (-a, -da, x, n);
 266         else
 267           retval = bsloww1 (-a, -da, x, n);
 268         break;
 269
 270       case 1:
 271       case 3:
 272         retval = bsloww2 (a, da, x, n);
 273         break;
 274     }
 275   return retval;
 276 }
 277
 278 /*******************************************************************/
 279 /* An ultimate sin routine. Given an IEEE double machine number x   */
 280 /* it computes the correctly rounded (to nearest) value of sin(x)  */
 281 /*******************************************************************/
 282 double
 283 SECTION
 284 __sin (double x)
 285 {
 286   double xx, res, t, cor, y, s, c, sn, ssn, cs, ccs, xn, a, da, db, eps, xn1,
 287     xn2;
 288   mynumber u, v;
 289   int4 k, m, n;
 290   double retval = 0;
 291
 292   SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
 293
 294   u.x = x;
 295   m = u.i[HIGH_HALF];
 296   k = 0x7fffffff & m;           /* no sign           */
 297   if (k < 0x3e500000)           /* if x->0 =>sin(x)=x */
 298     retval = x;
 299  /*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/
 300   else if (k < 0x3fd00000)
 301     {
 302       xx = x * x;
 303       /* Taylor series.  */
 304       t = POLYNOMIAL (xx) * (xx * x);
 305       res = x + t;
 306       cor = (x - res) + t;
 307       retval = (res == res + 1.07 * cor) ? res : slow (x);
 308     }                           /*  else  if (k < 0x3fd00000)    */
 309 /*---------------------------- 0.25<|x|< 0.855469---------------------- */
 310   else if (k < 0x3feb6000)
 311     {
 312       u.x = (m > 0) ? big + x : big - x;
 313       y = (m > 0) ? x - (u.x - big) : x + (u.x - big);
 314       xx = y * y;
 315       s = y + y * xx * (sn3 + xx * sn5);
 316       c = xx * (cs2 + xx * (cs4 + xx * cs6));
 317       SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
 318       if (m <= 0)
 319         {
 320           sn = -sn;
 321           ssn = -ssn;
 322         }
 323       cor = (ssn + s * ccs - sn * c) + cs * s;
 324       res = sn + cor;
 325       cor = (sn - res) + cor;
 326       retval = (res == res + 1.096 * cor) ? res : slow1 (x);
 327     }                           /*   else  if (k < 0x3feb6000)    */
 328
 329 /*----------------------- 0.855469  <|x|<2.426265  ----------------------*/
 330   else if (k < 0x400368fd)
 331     {
 332
 333       y = (m > 0) ? hp0 - x : hp0 + x;
 334       if (y >= 0)
 335         {
 336           u.x = big + y;
 337           y = (y - (u.x - big)) + hp1;
 338         }
 339       else
 340         {
 341           u.x = big - y;
 342           y = (-hp1) - (y + (u.x - big));
 343         }
 344       res = do_cos (u, y, &cor);
 345       retval = (res == res + 1.020 * cor) ? ((m > 0) ? res : -res) : slow2 (x);
 346     }                           /*   else  if (k < 0x400368fd)    */
 347
 348 /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/
 349   else if (k < 0x419921FB)
 350     {
 351       t = (x * hpinv + toint);
 352       xn = t - toint;
 353       v.x = t;
 354       y = (x - xn * mp1) - xn * mp2;
 355       n = v.i[LOW_HALF] & 3;
 356       da = xn * mp3;
 357       a = y - da;
 358       da = (y - a) - da;
 359       eps = fabs (x) * 1.2e-30;
 360
 361       switch (n)
 362         {                       /* quarter of unit circle */
 363         case 0:
 364         case 2:
 365           xx = a * a;
 366           if (n)
 367             {
 368               a = -a;
 369               da = -da;
 370             }
 371           if (xx < 0.01588)
 372             {
 373               /* Taylor series.  */
 374               res = TAYLOR_SIN (xx, a, da, cor);
 375               cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps;
 376               retval = (res == res + cor) ? res : sloww (a, da, x);
 377             }
 378           else
 379             {
 380               if (a > 0)
 381                 m = 1;
 382               else
 383                 {
 384                   m = 0;
 385                   a = -a;
 386                   da = -da;
 387                 }
 388               u.x = big + a;
 389               y = a - (u.x - big);
 390               res = do_sin (u, y, da, &cor);
 391               cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps;
 392               retval = ((res == res + cor) ? ((m) ? res : -res)
 393                         : sloww1 (a, da, x, m));
 394             }
 395           break;
 396
 397         case 1:
 398         case 3:
 399           if (a < 0)
 400             {
 401               a = -a;
 402               da = -da;
 403             }
 404           u.x = big + a;
 405           y = a - (u.x - big) + da;
 406           res = do_cos (u, y, &cor);
 407           cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps;
 408           retval = ((res == res + cor) ? ((n & 2) ? -res : res)
 409                     : sloww2 (a, da, x, n));
 410           break;
 411         }
 412     }                           /*   else  if (k <  0x419921FB )    */
 413
 414 /*---------------------105414350 <|x|< 281474976710656 --------------------*/
 415   else if (k < 0x42F00000)
 416     {
 417       t = (x * hpinv + toint);
 418       xn = t - toint;
 419       v.x = t;
 420       xn1 = (xn + 8.0e22) - 8.0e22;
 421       xn2 = xn - xn1;
 422       y = ((((x - xn1 * mp1) - xn1 * mp2) - xn2 * mp1) - xn2 * mp2);
 423       n = v.i[LOW_HALF] & 3;
 424       da = xn1 * pp3;
 425       t = y - da;
 426       da = (y - t) - da;
 427       da = (da - xn2 * pp3) - xn * pp4;
 428       a = t + da;
 429       da = (t - a) + da;
 430       eps = 1.0e-24;
 431
 432       switch (n)
 433         {
 434         case 0:
 435         case 2:
 436           xx = a * a;
 437           if (n)
 438             {
 439               a = -a;
 440               da = -da;
 441             }
 442           if (xx < 0.01588)
 443             {
 444               /* Taylor series.  */
 445               res = TAYLOR_SIN (xx, a, da, cor);
 446               cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps;
 447               retval = (res == res + cor) ? res : bsloww (a, da, x, n);
 448             }
 449           else
 450             {
 451               double t;
 452               if (a > 0)
 453                 {
 454                   m = 1;
 455                   t = a;
 456                   db = da;
 457                 }
 458               else
 459                 {
 460                   m = 0;
 461                   t = -a;
 462                   db = -da;
 463                 }
 464               u.x = big + t;
 465               y = t - (u.x - big);
 466               res = do_sin (u, y, db, &cor);
 467               cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps;
 468               retval = ((res == res + cor) ? ((m) ? res : -res)
 469                         : bsloww1 (a, da, x, n));
 470             }
 471           break;
 472
 473         case 1:
 474         case 3:
 475           if (a < 0)
 476             {
 477               a = -a;
 478               da = -da;
 479             }
 480           u.x = big + a;
 481           y = a - (u.x - big) + da;
 482           res = do_cos (u, y, &cor);
 483           cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps;
 484           retval = ((res == res + cor) ? ((n & 2) ? -res : res)
 485                     : bsloww2 (a, da, x, n));
 486           break;
 487         }
 488     }                           /*   else  if (k <  0x42F00000 )   */
 489
 490 /* -----------------281474976710656 <|x| <2^1024----------------------------*/
 491   else if (k < 0x7ff00000)
 492     retval = reduce_and_compute (x, 0);
 493
 494 /*--------------------- |x| > 2^1024 ----------------------------------*/
 495   else
 496     {
 497       if (k == 0x7ff00000 && u.i[LOW_HALF] == 0)
 498         __set_errno (EDOM);
 499       retval = x / x;
 500     }
 501
 502   return retval;
 503 }
 504
 505
 506 /*******************************************************************/
 507 /* An ultimate cos routine. Given an IEEE double machine number x   */
 508 /* it computes the correctly rounded (to nearest) value of cos(x)  */
 509 /*******************************************************************/
 510
 511 double
 512 SECTION
 513 __cos (double x)
 514 {
 515   double y, xx, res, t, cor, xn, a, da, db, eps, xn1,
 516     xn2;
 517   mynumber u, v;
 518   int4 k, m, n;
 519
 520   double retval = 0;
 521
 522   SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
 523
 524   u.x = x;
 525   m = u.i[HIGH_HALF];
 526   k = 0x7fffffff & m;
 527
 528   /* |x|<2^-27 => cos(x)=1 */
 529   if (k < 0x3e400000)
 530     retval = 1.0;
 531
 532   else if (k < 0x3feb6000)
 533     {                           /* 2^-27 < |x| < 0.855469 */
 534       y = fabs (x);
 535       u.x = big + y;
 536       y = y - (u.x - big);
 537       res = do_cos (u, y, &cor);
 538       retval = (res == res + 1.020 * cor) ? res : cslow2 (x);
 539     }                           /*   else  if (k < 0x3feb6000)    */
 540
 541   else if (k < 0x400368fd)
 542     { /* 0.855469  <|x|<2.426265  */ ;
 543       y = hp0 - fabs (x);
 544       a = y + hp1;
 545       da = (y - a) + hp1;
 546       xx = a * a;
 547       if (xx < 0.01588)
 548         {
 549           res = TAYLOR_SIN (xx, a, da, cor);
 550           cor = (cor > 0) ? 1.02 * cor + 1.0e-31 : 1.02 * cor - 1.0e-31;
 551           retval = (res == res + cor) ? res : csloww (a, da, x);
 552         }
 553       else
 554         {
 555           if (a > 0)
 556             {
 557               m = 1;
 558             }
 559           else
 560             {
 561               m = 0;
 562               a = -a;
 563               da = -da;
 564             }
 565           u.x = big + a;
 566           y = a - (u.x - big);
 567           res = do_sin (u, y, da, &cor);
 568           cor = (cor > 0) ? 1.035 * cor + 1.0e-31 : 1.035 * cor - 1.0e-31;
 569           retval = ((res == res + cor) ? ((m) ? res : -res)
 570                     : csloww1 (a, da, x, m));
 571         }
 572
 573     }                           /*   else  if (k < 0x400368fd)    */
 574
 575
 576   else if (k < 0x419921FB)
 577     {                           /* 2.426265<|x|< 105414350 */
 578       t = (x * hpinv + toint);
 579       xn = t - toint;
 580       v.x = t;
 581       y = (x - xn * mp1) - xn * mp2;
 582       n = v.i[LOW_HALF] & 3;
 583       da = xn * mp3;
 584       a = y - da;
 585       da = (y - a) - da;
 586       eps = fabs (x) * 1.2e-30;
 587
 588       switch (n)
 589         {
 590         case 1:
 591         case 3:
 592           xx = a * a;
 593           if (n == 1)
 594             {
 595               a = -a;
 596               da = -da;
 597             }
 598           if (xx < 0.01588)
 599             {
 600               res = TAYLOR_SIN (xx, a, da, cor);
 601               cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps;
 602               retval = (res == res + cor) ? res : csloww (a, da, x);
 603             }
 604           else
 605             {
 606               if (a > 0)
 607                 {
 608                   m = 1;
 609                 }
 610               else
 611                 {
 612                   m = 0;
 613                   a = -a;
 614                   da = -da;
 615                 }
 616               u.x = big + a;
 617               y = a - (u.x - big);
 618               res = do_sin (u, y, da, &cor);
 619               cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps;
 620               retval = ((res == res + cor) ? ((m) ? res : -res)
 621                         : csloww1 (a, da, x, m));
 622             }
 623           break;
 624
 625         case 0:
 626         case 2:
 627           if (a < 0)
 628             {
 629               a = -a;
 630               da = -da;
 631             }
 632           u.x = big + a;
 633           y = a - (u.x - big) + da;
 634           res = do_cos (u, y, &cor);
 635           cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps;
 636           retval = ((res == res + cor) ? ((n) ? -res : res)
 637                     : csloww2 (a, da, x, n));
 638           break;
 639         }
 640     }                           /*   else  if (k <  0x419921FB )    */
 641
 642   else if (k < 0x42F00000)
 643     {
 644       t = (x * hpinv + toint);
 645       xn = t - toint;
 646       v.x = t;
 647       xn1 = (xn + 8.0e22) - 8.0e22;
 648       xn2 = xn - xn1;
 649       y = ((((x - xn1 * mp1) - xn1 * mp2) - xn2 * mp1) - xn2 * mp2);
 650       n = v.i[LOW_HALF] & 3;
 651       da = xn1 * pp3;
 652       t = y - da;
 653       da = (y - t) - da;
 654       da = (da - xn2 * pp3) - xn * pp4;
 655       a = t + da;
 656       da = (t - a) + da;
 657       eps = 1.0e-24;
 658
 659       switch (n)
 660         {
 661         case 1:
 662         case 3:
 663           xx = a * a;
 664           if (n == 1)
 665             {
 666               a = -a;
 667               da = -da;
 668             }
 669           if (xx < 0.01588)
 670             {
 671               res = TAYLOR_SIN (xx, a, da, cor);
 672               cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps;
 673               retval = (res == res + cor) ? res : bsloww (a, da, x, n);
 674             }
 675           else
 676             {
 677               double t;
 678               if (a > 0)
 679                 {
 680                   m = 1;
 681                   t = a;
 682                   db = da;
 683                 }
 684               else
 685                 {
 686                   m = 0;
 687                   t = -a;
 688                   db = -da;
 689                 }
 690               u.x = big + t;
 691               y = t - (u.x - big);
 692               res = do_sin (u, y, db, &cor);
 693               cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps;
 694               retval = ((res == res + cor) ? ((m) ? res : -res)
 695                         : bsloww1 (a, da, x, n));
 696             }
 697           break;
 698
 699         case 0:
 700         case 2:
 701           if (a < 0)
 702             {
 703               a = -a;
 704               da = -da;
 705             }
 706           u.x = big + a;
 707           y = a - (u.x - big) + da;
 708           res = do_cos (u, y, &cor);
 709           cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps;
 710           retval = ((res == res + cor) ? ((n) ? -res : res)
 711                     : bsloww2 (a, da, x, n));
 712           break;
 713         }
 714     }                           /*   else  if (k <  0x42F00000 )    */
 715
 716   /* 281474976710656 <|x| <2^1024 */
 717   else if (k < 0x7ff00000)
 718     retval = reduce_and_compute (x, 1);
 719
 720   else
 721     {
 722       if (k == 0x7ff00000 && u.i[LOW_HALF] == 0)
 723         __set_errno (EDOM);
 724       retval = x / x;           /* |x| > 2^1024 */
 725     }
 726
 727   return retval;
 728 }
 729
 730 /************************************************************************/
 731 /*  Routine compute sin(x) for  2^-26 < |x|< 0.25 by  Taylor with more   */
 732 /* precision  and if still doesn't accurate enough by mpsin   or dubsin */
 733 /************************************************************************/
 734
 735 static double
 736 SECTION
 737 slow (double x)
 738 {
 739   double res, cor, w[2];
 740   res = TAYLOR_SLOW (x, 0, cor);
 741   if (res == res + 1.0007 * cor)
 742     return res;
 743   else
 744     {
 745       __dubsin (fabs (x), 0, w);
 746       if (w[0] == w[0] + 1.000000001 * w[1])
 747         return (x > 0) ? w[0] : -w[0];
 748       else
 749         return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false);
 750     }
 751 }
 752
 753 /*******************************************************************************/
 754 /* Routine compute sin(x) for 0.25<|x|< 0.855469 by __sincostab.tbl and Taylor */
 755 /* and if result still doesn't accurate enough by mpsin   or dubsin            */
 756 /*******************************************************************************/
 757
 758 static double
 759 SECTION
 760 slow1 (double x)
 761 {
 762   mynumber u;
 763   double w[2], y, cor, res;
 764   y = fabs (x);
 765   u.x = big + y;
 766   y = y - (u.x - big);
 767   res = do_sin_slow (u, y, 0, 0, &cor);
 768   if (res == res + cor)
 769     return (x > 0) ? res : -res;
 770   else
 771     {
 772       __dubsin (fabs (x), 0, w);
 773       if (w[0] == w[0] + 1.000000005 * w[1])
 774         return (x > 0) ? w[0] : -w[0];
 775       else
 776         return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false);
 777     }
 778 }
 779
 780 /**************************************************************************/
 781 /*  Routine compute sin(x) for   0.855469  <|x|<2.426265  by  __sincostab.tbl  */
 782 /* and if result still doesn't accurate enough by mpsin   or dubsin       */
 783 /**************************************************************************/
 784 static double
 785 SECTION
 786 slow2 (double x)
 787 {
 788   mynumber u;
 789   double w[2], y, y1, y2, cor, res, del;
 790
 791   y = fabs (x);
 792   y = hp0 - y;
 793   if (y >= 0)
 794     {
 795       u.x = big + y;
 796       y = y - (u.x - big);
 797       del = hp1;
 798     }
 799   else
 800     {
 801       u.x = big - y;
 802       y = -(y + (u.x - big));
 803       del = -hp1;
 804     }
 805   res = do_cos_slow (u, y, del, 0, &cor);
 806   if (res == res + cor)
 807     return (x > 0) ? res : -res;
 808   else
 809     {
 810       y = fabs (x) - hp0;
 811       y1 = y - hp1;
 812       y2 = (y - y1) - hp1;
 813       __docos (y1, y2, w);
 814       if (w[0] == w[0] + 1.000000005 * w[1])
 815         return (x > 0) ? w[0] : -w[0];
 816       else
 817         return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false);
 818     }
 819 }
 820
 821 /***************************************************************************/
 822 /*  Routine compute sin(x+dx) (Double-Length number) where x is small enough*/
 823 /* to use Taylor series around zero and   (x+dx)                            */
 824 /* in first or third quarter of unit circle.Routine receive also            */
 825 /* (right argument) the  original   value of x for computing error of      */
 826 /* result.And if result not accurate enough routine calls mpsin1 or dubsin */
 827 /***************************************************************************/
 828
 829 static double
 830 SECTION
 831 sloww (double x, double dx, double orig)
 832 {
 833   double y, t, res, cor, w[2], a, da, xn;
 834   mynumber v;
 835   int4 n;
 836   res = TAYLOR_SLOW (x, dx, cor);
 837   if (cor > 0)
 838     cor = 1.0005 * cor + fabs (orig) * 3.1e-30;
 839   else
 840     cor = 1.0005 * cor - fabs (orig) * 3.1e-30;
 841
 842   if (res == res + cor)
 843     return res;
 844   else
 845     {
 846       (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w);
 847       if (w[1] > 0)
 848         cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-30;
 849       else
 850         cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-30;
 851
 852       if (w[0] == w[0] + cor)
 853         return (x > 0) ? w[0] : -w[0];
 854       else
 855         {
 856           t = (orig * hpinv + toint);
 857           xn = t - toint;
 858           v.x = t;
 859           y = (orig - xn * mp1) - xn * mp2;
 860           n = v.i[LOW_HALF] & 3;
 861           da = xn * pp3;
 862           t = y - da;
 863           da = (y - t) - da;
 864           y = xn * pp4;
 865           a = t - y;
 866           da = ((t - a) - y) + da;
 867           if (n & 2)
 868             {
 869               a = -a;
 870               da = -da;
 871             }
 872           (a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w);
 873           if (w[1] > 0)
 874             cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-40;
 875           else
 876             cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-40;
 877
 878           if (w[0] == w[0] + cor)
 879             return (a > 0) ? w[0] : -w[0];
 880           else
 881             return __mpsin (orig, 0, true);
 882         }
 883     }
 884 }
 885
 886 /***************************************************************************/
 887 /*  Routine compute sin(x+dx)   (Double-Length number) where x in first or  */
 888 /*  third quarter of unit circle.Routine receive also (right argument) the  */
 889 /*  original   value of x for computing error of result.And if result not  */
 890 /* accurate enough routine calls  mpsin1   or dubsin                       */
 891 /***************************************************************************/
 892
 893 static double
 894 SECTION
 895 sloww1 (double x, double dx, double orig, int m)
 896 {
 897   mynumber u;
 898   double w[2], y, cor, res;
 899
 900   u.x = big + x;
 901   y = x - (u.x - big);
 902   res = do_sin_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor);
 903
 904   if (res == res + cor)
 905     return (m > 0) ? res : -res;
 906   else
 907     {
 908       __dubsin (x, dx, w);
 909
 910       if (w[1] > 0)
 911         cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig);
 912       else
 913         cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig);
 914
 915       if (w[0] == w[0] + cor)
 916         return (m > 0) ? w[0] : -w[0];
 917       else
 918         return __mpsin (orig, 0, true);
 919     }
 920 }
 921
 922 /***************************************************************************/
 923 /*  Routine compute sin(x+dx)   (Double-Length number) where x in second or */
 924 /*  fourth quarter of unit circle.Routine receive also  the  original value */
 925 /* and quarter(n= 1or 3)of x for computing error of result.And if result not*/
 926 /* accurate enough routine calls  mpsin1   or dubsin                       */
 927 /***************************************************************************/
 928
 929 static double
 930 SECTION
 931 sloww2 (double x, double dx, double orig, int n)
 932 {
 933   mynumber u;
 934   double w[2], y, cor, res;
 935
 936   u.x = big + x;
 937   y = x - (u.x - big);
 938   res = do_cos_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor);
 939
 940   if (res == res + cor)
 941     return (n & 2) ? -res : res;
 942   else
 943     {
 944       __docos (x, dx, w);
 945
 946       if (w[1] > 0)
 947         cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig);
 948       else
 949         cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig);
 950
 951       if (w[0] == w[0] + cor)
 952         return (n & 2) ? -w[0] : w[0];
 953       else
 954         return __mpsin (orig, 0, true);
 955     }
 956 }
 957
 958 /***************************************************************************/
 959 /*  Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x   */
 960 /* is small enough to use Taylor series around zero and   (x+dx)            */
 961 /* in first or third quarter of unit circle.Routine receive also            */
 962 /* (right argument) the  original   value of x for computing error of      */
 963 /* result.And if result not accurate enough routine calls other routines    */
 964 /***************************************************************************/
 965
 966 static double
 967 SECTION
 968 bsloww (double x, double dx, double orig, int n)
 969 {
 970   double res, cor, w[2];
 971
 972   res = TAYLOR_SLOW (x, dx, cor);
 973   cor = (cor > 0) ? 1.0005 * cor + 1.1e-24 : 1.0005 * cor - 1.1e-24;
 974   if (res == res + cor)
 975     return res;
 976   else
 977     {
 978       (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w);
 979       if (w[1] > 0)
 980         cor = 1.000000001 * w[1] + 1.1e-24;
 981       else
 982         cor = 1.000000001 * w[1] - 1.1e-24;
 983       if (w[0] == w[0] + cor)
 984         return (x > 0) ? w[0] : -w[0];
 985       else
 986         return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true);
 987     }
 988 }
 989
 990 /***************************************************************************/
 991 /*  Routine compute sin(x+dx)  or cos(x+dx) (Double-Length number) where x  */
 992 /* in first or third quarter of unit circle.Routine receive also            */
 993 /* (right argument) the original  value of x for computing error of result.*/
 994 /* And if result not  accurate enough routine calls  other routines         */
 995 /***************************************************************************/
 996
 997 static double
 998 SECTION
 999 bsloww1 (double x, double dx, double orig, int n)
1000 {
1001   mynumber u;
1002   double w[2], y, cor, res;
1003
1004   y = fabs (x);
1005   u.x = big + y;
1006   y = y - (u.x - big);
1007   dx = (x > 0) ? dx : -dx;
1008   res = do_sin_slow (u, y, dx, 1.1e-24, &cor);
1009   if (res == res + cor)
1010     return (x > 0) ? res : -res;
1011   else
1012     {
1013       __dubsin (fabs (x), dx, w);
1014
1015       if (w[1] > 0)
1016         cor = 1.000000005 * w[1] + 1.1e-24;
1017       else
1018         cor = 1.000000005 * w[1] - 1.1e-24;
1019
1020       if (w[0] == w[0] + cor)
1021         return (x > 0) ? w[0] : -w[0];
1022       else
1023         return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true);
1024     }
1025 }
1026
1027 /***************************************************************************/
1028 /*  Routine compute sin(x+dx)  or cos(x+dx) (Double-Length number) where x  */
1029 /* in second or fourth quarter of unit circle.Routine receive also  the     */
1030 /* original value and quarter(n= 1or 3)of x for computing error of result.  */
1031 /* And if result not accurate enough routine calls  other routines          */
1032 /***************************************************************************/
1033
1034 static double
1035 SECTION
1036 bsloww2 (double x, double dx, double orig, int n)
1037 {
1038   mynumber u;
1039   double w[2], y, cor, res;
1040
1041   y = fabs (x);
1042   u.x = big + y;
1043   y = y - (u.x - big);
1044   dx = (x > 0) ? dx : -dx;
1045   res = do_cos_slow (u, y, dx, 1.1e-24, &cor);
1046   if (res == res + cor)
1047     return (n & 2) ? -res : res;
1048   else
1049     {
1050       __docos (fabs (x), dx, w);
1051
1052       if (w[1] > 0)
1053         cor = 1.000000005 * w[1] + 1.1e-24;
1054       else
1055         cor = 1.000000005 * w[1] - 1.1e-24;
1056
1057       if (w[0] == w[0] + cor)
1058         return (n & 2) ? -w[0] : w[0];
1059       else
1060         return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true);
1061     }
1062 }
1063
1064 /************************************************************************/
1065 /*  Routine compute cos(x) for  2^-27 < |x|< 0.25 by  Taylor with more   */
1066 /* precision  and if still doesn't accurate enough by mpcos   or docos  */
1067 /************************************************************************/
1068
1069 static double
1070 SECTION
1071 cslow2 (double x)
1072 {
1073   mynumber u;
1074   double w[2], y, cor, res;
1075
1076   y = fabs (x);
1077   u.x = big + y;
1078   y = y - (u.x - big);
1079   res = do_cos_slow (u, y, 0, 0, &cor);
1080   if (res == res + cor)
1081     return res;
1082   else
1083     {
1084       y = fabs (x);
1085       __docos (y, 0, w);
1086       if (w[0] == w[0] + 1.000000005 * w[1])
1087         return w[0];
1088       else
1089         return __mpcos (x, 0, false);
1090     }
1091 }
1092
1093 /***************************************************************************/
1094 /*  Routine compute cos(x+dx) (Double-Length number) where x is small enough*/
1095 /* to use Taylor series around zero and   (x+dx) .Routine receive also      */
1096 /* (right argument) the  original   value of x for computing error of      */
1097 /* result.And if result not accurate enough routine calls other routines    */
1098 /***************************************************************************/
1099
1100
1101 static double
1102 SECTION
1103 csloww (double x, double dx, double orig)
1104 {
1105   double y, t, res, cor, w[2], a, da, xn;
1106   mynumber v;
1107   int4 n;
1108
1109   /* Taylor series */
1110   res = TAYLOR_SLOW (x, dx, cor);
1111
1112   if (cor > 0)
1113     cor = 1.0005 * cor + fabs (orig) * 3.1e-30;
1114   else
1115     cor = 1.0005 * cor - fabs (orig) * 3.1e-30;
1116
1117   if (res == res + cor)
1118     return res;
1119   else
1120     {
1121       (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w);
1122
1123       if (w[1] > 0)
1124         cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-30;
1125       else
1126         cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-30;
1127
1128       if (w[0] == w[0] + cor)
1129         return (x > 0) ? w[0] : -w[0];
1130       else
1131         {
1132           t = (orig * hpinv + toint);
1133           xn = t - toint;
1134           v.x = t;
1135           y = (orig - xn * mp1) - xn * mp2;
1136           n = v.i[LOW_HALF] & 3;
1137           da = xn * pp3;
1138           t = y - da;
1139           da = (y - t) - da;
1140           y = xn * pp4;
1141           a = t - y;
1142           da = ((t - a) - y) + da;
1143           if (n == 1)
1144             {
1145               a = -a;
1146               da = -da;
1147             }
1148           (a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w);
1149
1150           if (w[1] > 0)
1151             cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-40;
1152           else
1153             cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-40;
1154
1155           if (w[0] == w[0] + cor)
1156             return (a > 0) ? w[0] : -w[0];
1157           else
1158             return __mpcos (orig, 0, true);
1159         }
1160     }
1161 }
1162
1163 /***************************************************************************/
1164 /*  Routine compute sin(x+dx)   (Double-Length number) where x in first or  */
1165 /*  third quarter of unit circle.Routine receive also (right argument) the  */
1166 /*  original   value of x for computing error of result.And if result not  */
1167 /* accurate enough routine calls  other routines                            */
1168 /***************************************************************************/
1169
1170 static double
1171 SECTION
1172 csloww1 (double x, double dx, double orig, int m)
1173 {
1174   mynumber u;
1175   double w[2], y, cor, res;
1176
1177   u.x = big + x;
1178   y = x - (u.x - big);
1179   res = do_sin_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor);
1180
1181   if (res == res + cor)
1182     return (m > 0) ? res : -res;
1183   else
1184     {
1185       __dubsin (x, dx, w);
1186       if (w[1] > 0)
1187         cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig);
1188       else
1189         cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig);
1190       if (w[0] == w[0] + cor)
1191         return (m > 0) ? w[0] : -w[0];
1192       else
1193         return __mpcos (orig, 0, true);
1194     }
1195 }
1196
1197
1198 /***************************************************************************/
1199 /*  Routine compute sin(x+dx)   (Double-Length number) where x in second or */
1200 /*  fourth quarter of unit circle.Routine receive also  the  original value */
1201 /* and quarter(n= 1or 3)of x for computing error of result.And if result not*/
1202 /* accurate enough routine calls  other routines                            */
1203 /***************************************************************************/
1204
1205 static double
1206 SECTION
1207 csloww2 (double x, double dx, double orig, int n)
1208 {
1209   mynumber u;
1210   double w[2], y, cor, res;
1211
1212   u.x = big + x;
1213   y = x - (u.x - big);
1214   res = do_cos_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor);
1215
1216   if (res == res + cor)
1217     return (n) ? -res : res;
1218   else
1219     {
1220       __docos (x, dx, w);
1221       if (w[1] > 0)
1222         cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig);
1223       else
1224         cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig);
1225       if (w[0] == w[0] + cor)
1226         return (n) ? -w[0] : w[0];
1227       else
1228         return __mpcos (orig, 0, true);
1229     }
1230 }
1231
1232 #ifndef __cos
1233 weak_alias (__cos, cos)
1234 # ifdef NO_LONG_DOUBLE
1235 strong_alias (__cos, __cosl)
1236 weak_alias (__cos, cosl)
1237 # endif
1238 #endif
1239 #ifndef __sin
1240 weak_alias (__sin, sin)
1241 # ifdef NO_LONG_DOUBLE
1242 strong_alias (__sin, __sinl)
1243 weak_alias (__sin, sinl)
1244 # endif
1245 #endif