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