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