sysdeps/ieee754/dbl-64/s_sin.c

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