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