sysdeps/ieee754/dbl-64/s_sin.c

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