sysdeps/ieee754/dbl-64/e_atan2.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2016 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 /*  MODULE_NAME: atnat2.c                                               */
  21 /*                                                                      */
  22 /*  FUNCTIONS: uatan2                                                   */
  23 /*             atan2Mp                                                  */
  24 /*             signArctan2                                              */
  25 /*             normalized                                               */
  26 /*                                                                      */
  27 /*  FILES NEEDED: dla.h endian.h mpa.h mydefs.h atnat2.h                */
  28 /*                mpatan.c mpatan2.c mpsqrt.c                           */
  29 /*                uatan.tbl                                             */
  30 /*                                                                      */
  31 /* An ultimate atan2() routine. Given two IEEE double machine numbers y,*/
  32 /* x it computes the correctly rounded (to nearest) value of atan2(y,x).*/
  33 /*                                                                      */
  34 /* Assumption: Machine arithmetic operations are performed in           */
  35 /* round to nearest mode of IEEE 754 standard.                          */
  36 /*                                                                      */
  37 /************************************************************************/
  38
  39 #include <dla.h>
  40 #include "mpa.h"
  41 #include "MathLib.h"
  42 #include "uatan.tbl"
  43 #include "atnat2.h"
  44 #include <fenv.h>
  45 #include <float.h>
  46 #include <math.h>
  47 #include <math_private.h>
  48 #include <stap-probe.h>
  49
  50 #ifndef SECTION
  51 # define SECTION
  52 #endif
  53
  54 /************************************************************************/
  55 /* An ultimate atan2 routine. Given two IEEE double machine numbers y,x */
  56 /* it computes the correctly rounded (to nearest) value of atan2(y,x).  */
  57 /* Assumption: Machine arithmetic operations are performed in           */
  58 /* round to nearest mode of IEEE 754 standard.                          */
  59 /************************************************************************/
  60 static double atan2Mp (double, double, const int[]);
  61   /* Fix the sign and return after stage 1 or stage 2 */
  62 static double
  63 signArctan2 (double y, double z)
  64 {
  65   return __copysign (z, y);
  66 }
  67
  68 static double normalized (double, double, double, double);
  69 void __mpatan2 (mp_no *, mp_no *, mp_no *, int);
  70
  71 double
  72 SECTION
  73 __ieee754_atan2 (double y, double x)
  74 {
  75   int i, de, ux, dx, uy, dy;
  76   static const int pr[MM] = { 6, 8, 10, 20, 32 };
  77   double ax, ay, u, du, u9, ua, v, vv, dv, t1, t2, t3, t7, t8,
  78          z, zz, cor, s1, ss1, s2, ss2;
  79 #ifndef DLA_FMS
  80   double t4, t5, t6;
  81 #endif
  82   number num;
  83
  84   static const int ep = 59768832,      /*  57*16**5   */
  85                    em = -59768832;      /* -57*16**5   */
  86
  87   /* x=NaN or y=NaN */
  88   num.d = x;
  89   ux = num.i[HIGH_HALF];
  90   dx = num.i[LOW_HALF];
  91   if ((ux & 0x7ff00000) == 0x7ff00000)
  92     {
  93       if (((ux & 0x000fffff) | dx) != 0x00000000)
  94         return x + x;
  95     }
  96   num.d = y;
  97   uy = num.i[HIGH_HALF];
  98   dy = num.i[LOW_HALF];
  99   if ((uy & 0x7ff00000) == 0x7ff00000)
 100     {
 101       if (((uy & 0x000fffff) | dy) != 0x00000000)
 102         return y + y;
 103     }
 104
 105   /* y=+-0 */
 106   if (uy == 0x00000000)
 107     {
 108       if (dy == 0x00000000)
 109         {
 110           if ((ux & 0x80000000) == 0x00000000)
 111             return 0;
 112           else
 113             return opi.d;
 114         }
 115     }
 116   else if (uy == 0x80000000)
 117     {
 118       if (dy == 0x00000000)
 119         {
 120           if ((ux & 0x80000000) == 0x00000000)
 121             return -0.0;
 122           else
 123             return mopi.d;
 124         }
 125     }
 126
 127   /* x=+-0 */
 128   if (x == 0)
 129     {
 130       if ((uy & 0x80000000) == 0x00000000)
 131         return hpi.d;
 132       else
 133         return mhpi.d;
 134     }
 135
 136   /* x=+-INF */
 137   if (ux == 0x7ff00000)
 138     {
 139       if (dx == 0x00000000)
 140         {
 141           if (uy == 0x7ff00000)
 142             {
 143               if (dy == 0x00000000)
 144                 return qpi.d;
 145             }
 146           else if (uy == 0xfff00000)
 147             {
 148               if (dy == 0x00000000)
 149                 return mqpi.d;
 150             }
 151           else
 152             {
 153               if ((uy & 0x80000000) == 0x00000000)
 154                 return 0;
 155               else
 156                 return -0.0;
 157             }
 158         }
 159     }
 160   else if (ux == 0xfff00000)
 161     {
 162       if (dx == 0x00000000)
 163         {
 164           if (uy == 0x7ff00000)
 165             {
 166               if (dy == 0x00000000)
 167                 return tqpi.d;
 168             }
 169           else if (uy == 0xfff00000)
 170             {
 171               if (dy == 0x00000000)
 172                 return mtqpi.d;
 173             }
 174           else
 175             {
 176               if ((uy & 0x80000000) == 0x00000000)
 177                 return opi.d;
 178               else
 179                 return mopi.d;
 180             }
 181         }
 182     }
 183
 184   /* y=+-INF */
 185   if (uy == 0x7ff00000)
 186     {
 187       if (dy == 0x00000000)
 188         return hpi.d;
 189     }
 190   else if (uy == 0xfff00000)
 191     {
 192       if (dy == 0x00000000)
 193         return mhpi.d;
 194     }
 195
 196   SET_RESTORE_ROUND (FE_TONEAREST);
 197   /* either x/y or y/x is very close to zero */
 198   ax = (x < 0) ? -x : x;
 199   ay = (y < 0) ? -y : y;
 200   de = (uy & 0x7ff00000) - (ux & 0x7ff00000);
 201   if (de >= ep)
 202     {
 203       return ((y > 0) ? hpi.d : mhpi.d);
 204     }
 205   else if (de <= em)
 206     {
 207       if (x > 0)
 208         {
 209           double ret;
 210           if ((z = ay / ax) < TWOM1022)
 211             ret = normalized (ax, ay, y, z);
 212           else
 213             ret = signArctan2 (y, z);
 214           if (fabs (ret) < DBL_MIN)
 215             {
 216               double vret = ret ? ret : DBL_MIN;
 217               double force_underflow = vret * vret;
 218               math_force_eval (force_underflow);
 219             }
 220           return ret;
 221         }
 222       else
 223         {
 224           return ((y > 0) ? opi.d : mopi.d);
 225         }
 226     }
 227
 228   /* if either x or y is extremely close to zero, scale abs(x), abs(y). */
 229   if (ax < twom500.d || ay < twom500.d)
 230     {
 231       ax *= two500.d;
 232       ay *= two500.d;
 233     }
 234
 235   /* Likewise for large x and y.  */
 236   if (ax > two500.d || ay > two500.d)
 237     {
 238       ax *= twom500.d;
 239       ay *= twom500.d;
 240     }
 241
 242   /* x,y which are neither special nor extreme */
 243   if (ay < ax)
 244     {
 245       u = ay / ax;
 246       EMULV (ax, u, v, vv, t1, t2, t3, t4, t5);
 247       du = ((ay - v) - vv) / ax;
 248     }
 249   else
 250     {
 251       u = ax / ay;
 252       EMULV (ay, u, v, vv, t1, t2, t3, t4, t5);
 253       du = ((ax - v) - vv) / ay;
 254     }
 255
 256   if (x > 0)
 257     {
 258       /* (i)   x>0, abs(y)< abs(x):  atan(ay/ax) */
 259       if (ay < ax)
 260         {
 261           if (u < inv16.d)
 262             {
 263               v = u * u;
 264
 265               zz = du + u * v * (d3.d
 266                                  + v * (d5.d
 267                                         + v * (d7.d
 268                                                + v * (d9.d
 269                                                       + v * (d11.d
 270                                                              + v * d13.d)))));
 271
 272               if ((z = u + (zz - u1.d * u)) == u + (zz + u1.d * u))
 273                 return signArctan2 (y, z);
 274
 275               MUL2 (u, du, u, du, v, vv, t1, t2, t3, t4, t5, t6, t7, t8);
 276               s1 = v * (f11.d + v * (f13.d
 277                                      + v * (f15.d + v * (f17.d + v * f19.d))));
 278               ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2);
 279               MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 280               ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2);
 281               MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 282               ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2);
 283               MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 284               ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2);
 285               MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 286               MUL2 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
 287               ADD2 (u, du, s2, ss2, s1, ss1, t1, t2);
 288
 289               if ((z = s1 + (ss1 - u5.d * s1)) == s1 + (ss1 + u5.d * s1))
 290                 return signArctan2 (y, z);
 291
 292               return atan2Mp (x, y, pr);
 293             }
 294
 295           i = (TWO52 + TWO8 * u) - TWO52;
 296           i -= 16;
 297           t3 = u - cij[i][0].d;
 298           EADD (t3, du, v, dv);
 299           t1 = cij[i][1].d;
 300           t2 = cij[i][2].d;
 301           zz = v * t2 + (dv * t2
 302                          + v * v * (cij[i][3].d
 303                                     + v * (cij[i][4].d
 304                                            + v * (cij[i][5].d
 305                                                   + v * cij[i][6].d))));
 306           if (i < 112)
 307             {
 308               if (i < 48)
 309                 u9 = u91.d;     /* u < 1/4      */
 310               else
 311                 u9 = u92.d;
 312             }           /* 1/4 <= u < 1/2 */
 313           else
 314             {
 315               if (i < 176)
 316                 u9 = u93.d;     /* 1/2 <= u < 3/4 */
 317               else
 318                 u9 = u94.d;
 319             }           /* 3/4 <= u <= 1  */
 320           if ((z = t1 + (zz - u9 * t1)) == t1 + (zz + u9 * t1))
 321             return signArctan2 (y, z);
 322
 323           t1 = u - hij[i][0].d;
 324           EADD (t1, du, v, vv);
 325           s1 = v * (hij[i][11].d
 326                     + v * (hij[i][12].d
 327                            + v * (hij[i][13].d
 328                                   + v * (hij[i][14].d
 329                                          + v * hij[i][15].d))));
 330           ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2);
 331           MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 332           ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2);
 333           MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 334           ADD2 (hij[i][5].d, hij[i][6].d, s1, ss1, s2, ss2, t1, t2);
 335           MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 336           ADD2 (hij[i][3].d, hij[i][4].d, s1, ss1, s2, ss2, t1, t2);
 337           MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 338           ADD2 (hij[i][1].d, hij[i][2].d, s1, ss1, s2, ss2, t1, t2);
 339
 340           if ((z = s2 + (ss2 - ub.d * s2)) == s2 + (ss2 + ub.d * s2))
 341             return signArctan2 (y, z);
 342           return atan2Mp (x, y, pr);
 343         }
 344
 345       /* (ii)  x>0, abs(x)<=abs(y):  pi/2-atan(ax/ay) */
 346       if (u < inv16.d)
 347         {
 348           v = u * u;
 349           zz = u * v * (d3.d
 350                         + v * (d5.d
 351                                + v * (d7.d
 352                                       + v * (d9.d
 353                                              + v * (d11.d
 354                                                     + v * d13.d)))));
 355           ESUB (hpi.d, u, t2, cor);
 356           t3 = ((hpi1.d + cor) - du) - zz;
 357           if ((z = t2 + (t3 - u2.d)) == t2 + (t3 + u2.d))
 358             return signArctan2 (y, z);
 359
 360           MUL2 (u, du, u, du, v, vv, t1, t2, t3, t4, t5, t6, t7, t8);
 361           s1 = v * (f11.d
 362                     + v * (f13.d
 363                            + v * (f15.d + v * (f17.d + v * f19.d))));
 364           ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2);
 365           MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 366           ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2);
 367           MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 368           ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2);
 369           MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 370           ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2);
 371           MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 372           MUL2 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
 373           ADD2 (u, du, s2, ss2, s1, ss1, t1, t2);
 374           SUB2 (hpi.d, hpi1.d, s1, ss1, s2, ss2, t1, t2);
 375
 376           if ((z = s2 + (ss2 - u6.d)) == s2 + (ss2 + u6.d))
 377             return signArctan2 (y, z);
 378           return atan2Mp (x, y, pr);
 379         }
 380
 381       i = (TWO52 + TWO8 * u) - TWO52;
 382       i -= 16;
 383       v = (u - cij[i][0].d) + du;
 384
 385       zz = hpi1.d - v * (cij[i][2].d
 386                          + v * (cij[i][3].d
 387                                 + v * (cij[i][4].d
 388                                        + v * (cij[i][5].d
 389                                               + v * cij[i][6].d))));
 390       t1 = hpi.d - cij[i][1].d;
 391       if (i < 112)
 392         ua = ua1.d;     /* w <  1/2 */
 393       else
 394         ua = ua2.d;     /* w >= 1/2 */
 395       if ((z = t1 + (zz - ua)) == t1 + (zz + ua))
 396         return signArctan2 (y, z);
 397
 398       t1 = u - hij[i][0].d;
 399       EADD (t1, du, v, vv);
 400
 401       s1 = v * (hij[i][11].d
 402                 + v * (hij[i][12].d
 403                        + v * (hij[i][13].d
 404                               + v * (hij[i][14].d
 405                                      + v * hij[i][15].d))));
 406
 407       ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2);
 408       MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 409       ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2);
 410       MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 411       ADD2 (hij[i][5].d, hij[i][6].d, s1, ss1, s2, ss2, t1, t2);
 412       MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 413       ADD2 (hij[i][3].d, hij[i][4].d, s1, ss1, s2, ss2, t1, t2);
 414       MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 415       ADD2 (hij[i][1].d, hij[i][2].d, s1, ss1, s2, ss2, t1, t2);
 416       SUB2 (hpi.d, hpi1.d, s2, ss2, s1, ss1, t1, t2);
 417
 418       if ((z = s1 + (ss1 - uc.d)) == s1 + (ss1 + uc.d))
 419         return signArctan2 (y, z);
 420       return atan2Mp (x, y, pr);
 421     }
 422
 423   /* (iii) x<0, abs(x)< abs(y):  pi/2+atan(ax/ay) */
 424   if (ax < ay)
 425     {
 426       if (u < inv16.d)
 427         {
 428           v = u * u;
 429           zz = u * v * (d3.d
 430                         + v * (d5.d
 431                                + v * (d7.d
 432                                       + v * (d9.d
 433                                              + v * (d11.d + v * d13.d)))));
 434           EADD (hpi.d, u, t2, cor);
 435           t3 = ((hpi1.d + cor) + du) + zz;
 436           if ((z = t2 + (t3 - u3.d)) == t2 + (t3 + u3.d))
 437             return signArctan2 (y, z);
 438
 439           MUL2 (u, du, u, du, v, vv, t1, t2, t3, t4, t5, t6, t7, t8);
 440           s1 = v * (f11.d
 441                     + v * (f13.d + v * (f15.d + v * (f17.d + v * f19.d))));
 442           ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2);
 443           MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 444           ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2);
 445           MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 446           ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2);
 447           MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 448           ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2);
 449           MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 450           MUL2 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
 451           ADD2 (u, du, s2, ss2, s1, ss1, t1, t2);
 452           ADD2 (hpi.d, hpi1.d, s1, ss1, s2, ss2, t1, t2);
 453
 454           if ((z = s2 + (ss2 - u7.d)) == s2 + (ss2 + u7.d))
 455             return signArctan2 (y, z);
 456           return atan2Mp (x, y, pr);
 457         }
 458
 459       i = (TWO52 + TWO8 * u) - TWO52;
 460       i -= 16;
 461       v = (u - cij[i][0].d) + du;
 462       zz = hpi1.d + v * (cij[i][2].d
 463                          + v * (cij[i][3].d
 464                                 + v * (cij[i][4].d
 465                                        + v * (cij[i][5].d
 466                                               + v * cij[i][6].d))));
 467       t1 = hpi.d + cij[i][1].d;
 468       if (i < 112)
 469         ua = ua1.d;     /* w <  1/2 */
 470       else
 471         ua = ua2.d;     /* w >= 1/2 */
 472       if ((z = t1 + (zz - ua)) == t1 + (zz + ua))
 473         return signArctan2 (y, z);
 474
 475       t1 = u - hij[i][0].d;
 476       EADD (t1, du, v, vv);
 477       s1 = v * (hij[i][11].d
 478                 + v * (hij[i][12].d
 479                        + v * (hij[i][13].d
 480                               + v * (hij[i][14].d
 481                                      + v * hij[i][15].d))));
 482       ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2);
 483       MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 484       ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2);
 485       MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 486       ADD2 (hij[i][5].d, hij[i][6].d, s1, ss1, s2, ss2, t1, t2);
 487       MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 488       ADD2 (hij[i][3].d, hij[i][4].d, s1, ss1, s2, ss2, t1, t2);
 489       MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 490       ADD2 (hij[i][1].d, hij[i][2].d, s1, ss1, s2, ss2, t1, t2);
 491       ADD2 (hpi.d, hpi1.d, s2, ss2, s1, ss1, t1, t2);
 492
 493       if ((z = s1 + (ss1 - uc.d)) == s1 + (ss1 + uc.d))
 494         return signArctan2 (y, z);
 495       return atan2Mp (x, y, pr);
 496     }
 497
 498   /* (iv)  x<0, abs(y)<=abs(x):  pi-atan(ax/ay) */
 499   if (u < inv16.d)
 500     {
 501       v = u * u;
 502       zz = u * v * (d3.d
 503                     + v * (d5.d
 504                            + v * (d7.d
 505                                   + v * (d9.d + v * (d11.d + v * d13.d)))));
 506       ESUB (opi.d, u, t2, cor);
 507       t3 = ((opi1.d + cor) - du) - zz;
 508       if ((z = t2 + (t3 - u4.d)) == t2 + (t3 + u4.d))
 509         return signArctan2 (y, z);
 510
 511       MUL2 (u, du, u, du, v, vv, t1, t2, t3, t4, t5, t6, t7, t8);
 512       s1 = v * (f11.d + v * (f13.d + v * (f15.d + v * (f17.d + v * f19.d))));
 513       ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2);
 514       MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 515       ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2);
 516       MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 517       ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2);
 518       MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 519       ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2);
 520       MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 521       MUL2 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
 522       ADD2 (u, du, s2, ss2, s1, ss1, t1, t2);
 523       SUB2 (opi.d, opi1.d, s1, ss1, s2, ss2, t1, t2);
 524
 525       if ((z = s2 + (ss2 - u8.d)) == s2 + (ss2 + u8.d))
 526         return signArctan2 (y, z);
 527       return atan2Mp (x, y, pr);
 528     }
 529
 530   i = (TWO52 + TWO8 * u) - TWO52;
 531   i -= 16;
 532   v = (u - cij[i][0].d) + du;
 533   zz = opi1.d - v * (cij[i][2].d
 534                      + v * (cij[i][3].d
 535                             + v * (cij[i][4].d
 536                                    + v * (cij[i][5].d + v * cij[i][6].d))));
 537   t1 = opi.d - cij[i][1].d;
 538   if (i < 112)
 539     ua = ua1.d; /* w <  1/2 */
 540   else
 541     ua = ua2.d; /* w >= 1/2 */
 542   if ((z = t1 + (zz - ua)) == t1 + (zz + ua))
 543     return signArctan2 (y, z);
 544
 545   t1 = u - hij[i][0].d;
 546
 547   EADD (t1, du, v, vv);
 548
 549   s1 = v * (hij[i][11].d
 550             + v * (hij[i][12].d
 551                    + v * (hij[i][13].d
 552                           + v * (hij[i][14].d + v * hij[i][15].d))));
 553
 554   ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2);
 555   MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 556   ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2);
 557   MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 558   ADD2 (hij[i][5].d, hij[i][6].d, s1, ss1, s2, ss2, t1, t2);
 559   MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 560   ADD2 (hij[i][3].d, hij[i][4].d, s1, ss1, s2, ss2, t1, t2);
 561   MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8);
 562   ADD2 (hij[i][1].d, hij[i][2].d, s1, ss1, s2, ss2, t1, t2);
 563   SUB2 (opi.d, opi1.d, s2, ss2, s1, ss1, t1, t2);
 564
 565   if ((z = s1 + (ss1 - uc.d)) == s1 + (ss1 + uc.d))
 566     return signArctan2 (y, z);
 567   return atan2Mp (x, y, pr);
 568 }
 569
 570 #ifndef __ieee754_atan2
 571 strong_alias (__ieee754_atan2, __atan2_finite)
 572 #endif
 573
 574 /* Treat the Denormalized case */
 575 static double
 576 SECTION
 577 normalized (double ax, double ay, double y, double z)
 578 {
 579   int p;
 580   mp_no mpx, mpy, mpz, mperr, mpz2, mpt1;
 581   p = 6;
 582   __dbl_mp (ax, &mpx, p);
 583   __dbl_mp (ay, &mpy, p);
 584   __dvd (&mpy, &mpx, &mpz, p);
 585   __dbl_mp (ue.d, &mpt1, p);
 586   __mul (&mpz, &mpt1, &mperr, p);
 587   __sub (&mpz, &mperr, &mpz2, p);
 588   __mp_dbl (&mpz2, &z, p);
 589   return signArctan2 (y, z);
 590 }
 591
 592 /* Stage 3: Perform a multi-Precision computation */
 593 static double
 594 SECTION
 595 atan2Mp (double x, double y, const int pr[])
 596 {
 597   double z1, z2;
 598   int i, p;
 599   mp_no mpx, mpy, mpz, mpz1, mpz2, mperr, mpt1;
 600   for (i = 0; i < MM; i++)
 601     {
 602       p = pr[i];
 603       __dbl_mp (x, &mpx, p);
 604       __dbl_mp (y, &mpy, p);
 605       __mpatan2 (&mpy, &mpx, &mpz, p);
 606       __dbl_mp (ud[i].d, &mpt1, p);
 607       __mul (&mpz, &mpt1, &mperr, p);
 608       __add (&mpz, &mperr, &mpz1, p);
 609       __sub (&mpz, &mperr, &mpz2, p);
 610       __mp_dbl (&mpz1, &z1, p);
 611       __mp_dbl (&mpz2, &z2, p);
 612       if (z1 == z2)
 613         {
 614           LIBC_PROBE (slowatan2, 4, &p, &x, &y, &z1);
 615           return z1;
 616         }
 617     }
 618   LIBC_PROBE (slowatan2_inexact, 4, &p, &x, &y, &z1);
 619   return z1;                    /*if impossible to do exact computing */
 620 }