1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2020 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 <https://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#include <libm-alias-finite.h>
52
53#ifndef SECTION
54# define SECTION
55#endif
56
57/************************************************************************/
58/* An ultimate atan2 routine. Given two IEEE double machine numbers y,x */
59/* it computes the correctly rounded (to nearest) value of atan2(y,x). */
60/* Assumption: Machine arithmetic operations are performed in */
61/* round to nearest mode of IEEE 754 standard. */
62/************************************************************************/
63static double atan2Mp (double, double, const int[]);
64 /* Fix the sign and return after stage 1 or stage 2 */
65static double
66signArctan2 (double y, double z)
67{
68 return copysign (z, y);
69}
70
71static double normalized (double, double, double, double);
72void __mpatan2 (mp_no *, mp_no *, mp_no *, int);
73
74double
75SECTION
76__ieee754_atan2 (double y, double x)
77{
78 int i, de, ux, dx, uy, dy;
79 static const int pr[MM] = { 6, 8, 10, 20, 32 };
80 double ax, ay, u, du, u9, ua, v, vv, dv, t1, t2, t3,
81 z, zz, cor, s1, ss1, s2, ss2;
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 + y;
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);
247 du = ((ay - v) - vv) / ax;
248 }
249 else
250 {
251 u = ax / ay;
252 EMULV (ay, u, v, vv);
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);
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);
280 ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2);
281 MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2);
282 ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2);
283 MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2);
284 ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2);
285 MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2);
286 MUL2 (u, du, s1, ss1, s2, ss2, t1, t2);
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);
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);
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);
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);
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);
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);
366 ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2);
367 MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2);
368 ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2);
369 MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2);
370 ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2);
371 MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2);
372 MUL2 (u, du, s1, ss1, s2, ss2, t1, t2);
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);
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);
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);
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);
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);
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);
444 ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2);
445 MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2);
446 ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2);
447 MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2);
448 ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2);
449 MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2);
450 MUL2 (u, du, s1, ss1, s2, ss2, t1, t2);
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);
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);
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);
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);
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);
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);
515 ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2);
516 MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2);
517 ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2);
518 MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2);
519 ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2);
520 MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2);
521 MUL2 (u, du, s1, ss1, s2, ss2, t1, t2);
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);
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);
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);
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);
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
571libm_alias_finite (__ieee754_atan2, __atan2)
572#endif
573
574/* Treat the Denormalized case */
575static double
576SECTION
577normalized (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 */
593static double
594SECTION
595atan2Mp (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}
621