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 | /************************************************************************/ |
62 | static double atan2Mp (double, double, const int[]); |
63 | /* Fix the sign and return after stage 1 or stage 2 */ |
64 | static double |
65 | signArctan2 (double y, double z) |
66 | { |
67 | return copysign (z, y); |
68 | } |
69 | |
70 | static double normalized (double, double, double, double); |
71 | void __mpatan2 (mp_no *, mp_no *, mp_no *, int); |
72 | |
73 | double |
74 | SECTION |
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 |
573 | strong_alias (__ieee754_atan2, __atan2_finite) |
574 | #endif |
575 | |
576 | /* Treat the Denormalized case */ |
577 | static double |
578 | SECTION |
579 | normalized (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 */ |
595 | static double |
596 | SECTION |
597 | atan2Mp (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 | |