]>
Commit | Line | Data |
---|---|---|
f7eac6eb | 1 | /* |
e4d82761 | 2 | * IBM Accurate Mathematical Library |
aeb25823 | 3 | * written by International Business Machines Corp. |
f7a9f785 | 4 | * Copyright (C) 2001-2016 Free Software Foundation, Inc. |
f7eac6eb | 5 | * |
e4d82761 UD |
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 | |
cc7375ce | 8 | * the Free Software Foundation; either version 2.1 of the License, or |
e4d82761 | 9 | * (at your option) any later version. |
f7eac6eb | 10 | * |
e4d82761 UD |
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 | |
c6c6dd48 | 14 | * GNU Lesser General Public License for more details. |
f7eac6eb | 15 | * |
e4d82761 | 16 | * You should have received a copy of the GNU Lesser General Public License |
59ba27a6 | 17 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
f7eac6eb | 18 | */ |
e4d82761 UD |
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 | ||
c8b3296b | 39 | #include <dla.h> |
e4d82761 UD |
40 | #include "mpa.h" |
41 | #include "MathLib.h" | |
42 | #include "uatan.tbl" | |
43 | #include "atnat2.h" | |
8431838d | 44 | #include <fenv.h> |
4629c866 JM |
45 | #include <float.h> |
46 | #include <math.h> | |
1ed0291c | 47 | #include <math_private.h> |
10e1cf6b | 48 | #include <stap-probe.h> |
e859d1d9 | 49 | |
31d3cc00 UD |
50 | #ifndef SECTION |
51 | # define SECTION | |
52 | #endif | |
53 | ||
e4d82761 UD |
54 | /************************************************************************/ |
55 | /* An ultimate atan2 routine. Given two IEEE double machine numbers y,x */ | |
56 | /* it computes the correctly rounded (to nearest) value of atan2(y,x). */ | |
57 | /* Assumption: Machine arithmetic operations are performed in */ | |
58 | /* round to nearest mode of IEEE 754 standard. */ | |
59 | /************************************************************************/ | |
1728ab37 | 60 | static double atan2Mp (double, double, const int[]); |
af968f62 | 61 | /* Fix the sign and return after stage 1 or stage 2 */ |
1728ab37 SP |
62 | static double |
63 | signArctan2 (double y, double z) | |
af968f62 | 64 | { |
1728ab37 | 65 | return __copysign (z, y); |
af968f62 | 66 | } |
1728ab37 SP |
67 | |
68 | static double normalized (double, double, double, double); | |
69 | void __mpatan2 (mp_no *, mp_no *, mp_no *, int); | |
e4d82761 | 70 | |
31d3cc00 UD |
71 | double |
72 | SECTION | |
1728ab37 SP |
73 | __ieee754_atan2 (double y, double x) |
74 | { | |
75 | int i, de, ux, dx, uy, dy; | |
76 | static const int pr[MM] = { 6, 8, 10, 20, 32 }; | |
77 | double ax, ay, u, du, u9, ua, v, vv, dv, t1, t2, t3, t7, t8, | |
c5d5d574 | 78 | z, zz, cor, s1, ss1, s2, ss2; |
58985aa9 | 79 | #ifndef DLA_FMS |
1728ab37 | 80 | double t4, t5, t6; |
50944bca | 81 | #endif |
e4d82761 | 82 | number num; |
e4d82761 | 83 | |
c5d5d574 OB |
84 | static const int ep = 59768832, /* 57*16**5 */ |
85 | em = -59768832; /* -57*16**5 */ | |
e4d82761 UD |
86 | |
87 | /* x=NaN or y=NaN */ | |
1728ab37 SP |
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 + x; | |
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 | } | |
e4d82761 UD |
104 | |
105 | /* y=+-0 */ | |
1728ab37 SP |
106 | if (uy == 0x00000000) |
107 | { | |
108 | if (dy == 0x00000000) | |
109 | { | |
110 | if ((ux & 0x80000000) == 0x00000000) | |
a64d7e0e | 111 | return 0; |
1728ab37 SP |
112 | else |
113 | return opi.d; | |
114 | } | |
115 | } | |
116 | else if (uy == 0x80000000) | |
117 | { | |
118 | if (dy == 0x00000000) | |
119 | { | |
120 | if ((ux & 0x80000000) == 0x00000000) | |
a64d7e0e | 121 | return -0.0; |
1728ab37 SP |
122 | else |
123 | return mopi.d; | |
124 | } | |
125 | } | |
e4d82761 UD |
126 | |
127 | /* x=+-0 */ | |
a64d7e0e | 128 | if (x == 0) |
1728ab37 SP |
129 | { |
130 | if ((uy & 0x80000000) == 0x00000000) | |
131 | return hpi.d; | |
132 | else | |
133 | return mhpi.d; | |
134 | } | |
e4d82761 UD |
135 | |
136 | /* x=+-INF */ | |
1728ab37 SP |
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) | |
a64d7e0e | 154 | return 0; |
1728ab37 | 155 | else |
a64d7e0e | 156 | return -0.0; |
1728ab37 SP |
157 | } |
158 | } | |
e4d82761 | 159 | } |
1728ab37 SP |
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 | } | |
e4d82761 | 182 | } |
e4d82761 UD |
183 | |
184 | /* y=+-INF */ | |
1728ab37 SP |
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 | } | |
e4d82761 | 195 | |
8431838d | 196 | SET_RESTORE_ROUND (FE_TONEAREST); |
e4d82761 | 197 | /* either x/y or y/x is very close to zero */ |
a64d7e0e SP |
198 | ax = (x < 0) ? -x : x; |
199 | ay = (y < 0) ? -y : y; | |
e4d82761 | 200 | de = (uy & 0x7ff00000) - (ux & 0x7ff00000); |
1728ab37 SP |
201 | if (de >= ep) |
202 | { | |
a64d7e0e | 203 | return ((y > 0) ? hpi.d : mhpi.d); |
1728ab37 SP |
204 | } |
205 | else if (de <= em) | |
206 | { | |
a64d7e0e | 207 | if (x > 0) |
1728ab37 | 208 | { |
4629c866 | 209 | double ret; |
1728ab37 | 210 | if ((z = ay / ax) < TWOM1022) |
4629c866 | 211 | ret = normalized (ax, ay, y, z); |
1728ab37 | 212 | else |
4629c866 JM |
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; | |
1728ab37 SP |
221 | } |
222 | else | |
223 | { | |
a64d7e0e | 224 | return ((y > 0) ? opi.d : mopi.d); |
1728ab37 SP |
225 | } |
226 | } | |
e4d82761 UD |
227 | |
228 | /* if either x or y is extremely close to zero, scale abs(x), abs(y). */ | |
1728ab37 SP |
229 | if (ax < twom500.d || ay < twom500.d) |
230 | { | |
231 | ax *= two500.d; | |
232 | ay *= two500.d; | |
233 | } | |
e4d82761 | 234 | |
7726d6a9 JM |
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 | ||
e4d82761 | 242 | /* x,y which are neither special nor extreme */ |
1728ab37 SP |
243 | if (ay < ax) |
244 | { | |
245 | u = ay / ax; | |
246 | EMULV (ax, u, v, vv, t1, t2, t3, t4, t5); | |
247 | du = ((ay - v) - vv) / ax; | |
248 | } | |
249 | else | |
250 | { | |
251 | u = ax / ay; | |
252 | EMULV (ay, u, v, vv, t1, t2, t3, t4, t5); | |
253 | du = ((ax - v) - vv) / ay; | |
e4d82761 UD |
254 | } |
255 | ||
a64d7e0e | 256 | if (x > 0) |
1728ab37 SP |
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, t3, t4, t5, t6, t7, t8); | |
276 | s1 = v * (f11.d + v * (f13.d | |
277 | + v * (f15.d + v * (f17.d + v * f19.d)))); | |
a64d7e0e | 278 | ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2); |
1728ab37 SP |
279 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
280 | ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2); | |
281 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); | |
282 | ADD2 (f5.d, ff5.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 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2); | |
285 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); | |
286 | MUL2 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8); | |
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) | |
c5d5d574 | 309 | u9 = u91.d; /* u < 1/4 */ |
1728ab37 SP |
310 | else |
311 | u9 = u92.d; | |
c5d5d574 | 312 | } /* 1/4 <= u < 1/2 */ |
1728ab37 SP |
313 | else |
314 | { | |
315 | if (i < 176) | |
c5d5d574 | 316 | u9 = u93.d; /* 1/2 <= u < 3/4 */ |
1728ab37 SP |
317 | else |
318 | u9 = u94.d; | |
c5d5d574 | 319 | } /* 3/4 <= u <= 1 */ |
1728ab37 SP |
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 | |
c5d5d574 OB |
327 | + v * (hij[i][13].d |
328 | + v * (hij[i][14].d | |
329 | + v * hij[i][15].d)))); | |
a64d7e0e | 330 | ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2); |
1728ab37 SP |
331 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
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, t3, t4, t5, t6, t7, t8); | |
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, t3, t4, t5, t6, t7, t8); | |
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, t3, t4, t5, t6, t7, t8); | |
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, t3, t4, t5, t6, t7, t8); | |
361 | s1 = v * (f11.d | |
362 | + v * (f13.d | |
363 | + v * (f15.d + v * (f17.d + v * f19.d)))); | |
a64d7e0e | 364 | ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2); |
1728ab37 SP |
365 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
366 | ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2); | |
367 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); | |
368 | ADD2 (f5.d, ff5.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 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2); | |
371 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); | |
372 | MUL2 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8); | |
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 | ||
a64d7e0e | 407 | ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2); |
1728ab37 SP |
408 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
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, t3, t4, t5, t6, t7, t8); | |
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, t3, t4, t5, t6, t7, t8); | |
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, t3, t4, t5, t6, t7, t8); | |
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); | |
e4d82761 | 421 | } |
1728ab37 SP |
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, t3, t4, t5, t6, t7, t8); | |
440 | s1 = v * (f11.d | |
441 | + v * (f13.d + v * (f15.d + v * (f17.d + v * f19.d)))); | |
a64d7e0e | 442 | ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2); |
1728ab37 SP |
443 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
444 | ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2); | |
445 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); | |
446 | ADD2 (f5.d, ff5.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 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2); | |
449 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); | |
450 | MUL2 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8); | |
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)))); | |
a64d7e0e | 482 | ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2); |
1728ab37 SP |
483 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
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, t3, t4, t5, t6, t7, t8); | |
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, t3, t4, t5, t6, t7, t8); | |
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, t3, t4, t5, t6, t7, t8); | |
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); | |
e4d82761 UD |
496 | } |
497 | ||
1728ab37 SP |
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, t3, t4, t5, t6, t7, t8); | |
512 | s1 = v * (f11.d + v * (f13.d + v * (f15.d + v * (f17.d + v * f19.d)))); | |
a64d7e0e | 513 | ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2); |
1728ab37 SP |
514 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
515 | ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2); | |
516 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); | |
517 | ADD2 (f5.d, ff5.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 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2); | |
520 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); | |
521 | MUL2 (u, du, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8); | |
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); | |
e4d82761 | 528 | } |
1728ab37 SP |
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 | ||
a64d7e0e | 554 | ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2); |
1728ab37 SP |
555 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
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, t3, t4, t5, t6, t7, t8); | |
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, t3, t4, t5, t6, t7, t8); | |
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, t3, t4, t5, t6, t7, t8); | |
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); | |
e4d82761 | 568 | } |
1728ab37 | 569 | |
af968f62 | 570 | #ifndef __ieee754_atan2 |
0ac5ae23 | 571 | strong_alias (__ieee754_atan2, __atan2_finite) |
af968f62 | 572 | #endif |
0ac5ae23 | 573 | |
1728ab37 | 574 | /* Treat the Denormalized case */ |
31d3cc00 UD |
575 | static double |
576 | SECTION | |
1728ab37 SP |
577 | normalized (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); | |
e4d82761 | 590 | } |
1728ab37 SP |
591 | |
592 | /* Stage 3: Perform a multi-Precision computation */ | |
31d3cc00 UD |
593 | static double |
594 | SECTION | |
1728ab37 | 595 | atan2Mp (double x, double y, const int pr[]) |
e4d82761 | 596 | { |
1728ab37 SP |
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) | |
10e1cf6b SP |
613 | { |
614 | LIBC_PROBE (slowatan2, 4, &p, &x, &y, &z1); | |
615 | return z1; | |
616 | } | |
1728ab37 | 617 | } |
10e1cf6b | 618 | LIBC_PROBE (slowatan2_inexact, 4, &p, &x, &y, &z1); |
1728ab37 | 619 | return z1; /*if impossible to do exact computing */ |
f7eac6eb | 620 | } |