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