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