sysdeps/ieee754/dbl-64/s_sin.c

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