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