sysdeps/ieee754/dbl-64/s_sin.c

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