1 /*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2023 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 /* signArctan2 */
24 /* */
25 /* FILES NEEDED: dla.h endian.h mydefs.h atnat2.h */
26 /* uatan.tbl */
27 /* */
28 /************************************************************************/
29
30 #include <dla.h>
31 #include "mydefs.h"
32 #include "uatan.tbl"
33 #include "atnat2.h"
34 #include <fenv.h>
35 #include <float.h>
36 #include <math.h>
37 #include <math-barriers.h>
38 #include <math_private.h>
39 #include <fenv_private.h>
40 #include <libm-alias-finite.h>
41
42 #ifndef SECTION
43 # define SECTION
44 #endif
45
46 #define TWO52 0x1.0p52
47 #define TWOM1022 0x1.0p-1022
48
49 /* Fix the sign and return after stage 1 or stage 2 */
50 static double
51 signArctan2 (double y, double z)
52 {
53 return copysign (z, y);
54 }
55
56 /* atan2 with max ULP of ~0.524 based on random sampling. */
57 double
58 SECTION
59 __ieee754_atan2 (double y, double x)
60 {
61 int i, de, ux, dx, uy, dy;
62 double ax, ay, u, du, v, vv, dv, t1, t2, t3,
63 z, zz, cor;
64 mynumber num;
65
66 static const int ep = 59768832, /* 57*16**5 */
67 em = -59768832; /* -57*16**5 */
68
69 /* x=NaN or y=NaN */
70 num.d = x;
71 ux = num.i[HIGH_HALF];
72 dx = num.i[LOW_HALF];
73 if ((ux & 0x7ff00000) == 0x7ff00000)
74 {
75 if (((ux & 0x000fffff) | dx) != 0x00000000)
76 return x + y;
77 }
78 num.d = y;
79 uy = num.i[HIGH_HALF];
80 dy = num.i[LOW_HALF];
81 if ((uy & 0x7ff00000) == 0x7ff00000)
82 {
83 if (((uy & 0x000fffff) | dy) != 0x00000000)
84 return y + y;
85 }
86
87 /* y=+-0 */
88 if (uy == 0x00000000)
89 {
90 if (dy == 0x00000000)
91 {
92 if ((ux & 0x80000000) == 0x00000000)
93 return 0;
94 else
95 return opi.d;
96 }
97 }
98 else if (uy == 0x80000000)
99 {
100 if (dy == 0x00000000)
101 {
102 if ((ux & 0x80000000) == 0x00000000)
103 return -0.0;
104 else
105 return mopi.d;
106 }
107 }
108
109 /* x=+-0 */
110 if (x == 0)
111 {
112 if ((uy & 0x80000000) == 0x00000000)
113 return hpi.d;
114 else
115 return mhpi.d;
116 }
117
118 /* x=+-INF */
119 if (ux == 0x7ff00000)
120 {
121 if (dx == 0x00000000)
122 {
123 if (uy == 0x7ff00000)
124 {
125 if (dy == 0x00000000)
126 return qpi.d;
127 }
128 else if (uy == 0xfff00000)
129 {
130 if (dy == 0x00000000)
131 return mqpi.d;
132 }
133 else
134 {
135 if ((uy & 0x80000000) == 0x00000000)
136 return 0;
137 else
138 return -0.0;
139 }
140 }
141 }
142 else if (ux == 0xfff00000)
143 {
144 if (dx == 0x00000000)
145 {
146 if (uy == 0x7ff00000)
147 {
148 if (dy == 0x00000000)
149 return tqpi.d;
150 }
151 else if (uy == 0xfff00000)
152 {
153 if (dy == 0x00000000)
154 return mtqpi.d;
155 }
156 else
157 {
158 if ((uy & 0x80000000) == 0x00000000)
159 return opi.d;
160 else
161 return mopi.d;
162 }
163 }
164 }
165
166 /* y=+-INF */
167 if (uy == 0x7ff00000)
168 {
169 if (dy == 0x00000000)
170 return hpi.d;
171 }
172 else if (uy == 0xfff00000)
173 {
174 if (dy == 0x00000000)
175 return mhpi.d;
176 }
177
178 SET_RESTORE_ROUND (FE_TONEAREST);
179 /* either x/y or y/x is very close to zero */
180 ax = (x < 0) ? -x : x;
181 ay = (y < 0) ? -y : y;
182 de = (uy & 0x7ff00000) - (ux & 0x7ff00000);
183 if (de >= ep)
184 {
185 return ((y > 0) ? hpi.d : mhpi.d);
186 }
187 else if (de <= em)
188 {
189 if (x > 0)
190 {
191 double ret;
192 z = ay / ax;
193 ret = signArctan2 (y, z);
194 if (fabs (ret) < DBL_MIN)
195 {
196 double vret = ret ? ret : DBL_MIN;
197 double force_underflow = vret * vret;
198 math_force_eval (force_underflow);
199 }
200 return ret;
201 }
202 else
203 {
204 return ((y > 0) ? opi.d : mopi.d);
205 }
206 }
207
208 /* if either x or y is extremely close to zero, scale abs(x), abs(y). */
209 if (ax < twom500.d || ay < twom500.d)
210 {
211 ax *= two500.d;
212 ay *= two500.d;
213 }
214
215 /* Likewise for large x and y. */
216 if (ax > two500.d || ay > two500.d)
217 {
218 ax *= twom500.d;
219 ay *= twom500.d;
220 }
221
222 /* x,y which are neither special nor extreme */
223 if (ay < ax)
224 {
225 u = ay / ax;
226 EMULV (ax, u, v, vv);
227 du = ((ay - v) - vv) / ax;
228 }
229 else
230 {
231 u = ax / ay;
232 EMULV (ay, u, v, vv);
233 du = ((ax - v) - vv) / ay;
234 }
235
236 if (x > 0)
237 {
238 /* (i) x>0, abs(y)< abs(x): atan(ay/ax) */
239 if (ay < ax)
240 {
241 if (u < inv16.d)
242 {
243 v = u * u;
244
245 zz = du + u * v * (d3.d
246 + v * (d5.d
247 + v * (d7.d
248 + v * (d9.d
249 + v * (d11.d
250 + v * d13.d)))));
251
252 z = u + zz;
253 /* Max ULP is 0.504. */
254 return signArctan2 (y, z);
255 }
256
257 i = (TWO52 + 256 * u) - TWO52;
258 i -= 16;
259 t3 = u - cij[i][0].d;
260 EADD (t3, du, v, dv);
261 t1 = cij[i][1].d;
262 t2 = cij[i][2].d;
263 zz = v * t2 + (dv * t2
264 + v * v * (cij[i][3].d
265 + v * (cij[i][4].d
266 + v * (cij[i][5].d
267 + v * cij[i][6].d))));
268 z = t1 + zz;
269 /* Max ULP is 0.56. */
270 return signArctan2 (y, z);
271 }
272
273 /* (ii) x>0, abs(x)<=abs(y): pi/2-atan(ax/ay) */
274 if (u < inv16.d)
275 {
276 v = u * u;
277 zz = u * v * (d3.d
278 + v * (d5.d
279 + v * (d7.d
280 + v * (d9.d
281 + v * (d11.d
282 + v * d13.d)))));
283 ESUB (hpi.d, u, t2, cor);
284 t3 = ((hpi1.d + cor) - du) - zz;
285 z = t2 + t3;
286 /* Max ULP is 0.501. */
287 return signArctan2 (y, z);
288 }
289
290 i = (TWO52 + 256 * u) - TWO52;
291 i -= 16;
292 v = (u - cij[i][0].d) + du;
293
294 zz = hpi1.d - v * (cij[i][2].d
295 + v * (cij[i][3].d
296 + v * (cij[i][4].d
297 + v * (cij[i][5].d
298 + v * cij[i][6].d))));
299 t1 = hpi.d - cij[i][1].d;
300 z = t1 + zz;
301 /* Max ULP is 0.503. */
302 return signArctan2 (y, z);
303 }
304
305 /* (iii) x<0, abs(x)< abs(y): pi/2+atan(ax/ay) */
306 if (ax < ay)
307 {
308 if (u < inv16.d)
309 {
310 v = u * u;
311 zz = u * v * (d3.d
312 + v * (d5.d
313 + v * (d7.d
314 + v * (d9.d
315 + v * (d11.d + v * d13.d)))));
316 EADD (hpi.d, u, t2, cor);
317 t3 = ((hpi1.d + cor) + du) + zz;
318 z = t2 + t3;
319 /* Max ULP is 0.501. */
320 return signArctan2 (y, z);
321 }
322
323 i = (TWO52 + 256 * u) - TWO52;
324 i -= 16;
325 v = (u - cij[i][0].d) + du;
326 zz = hpi1.d + v * (cij[i][2].d
327 + v * (cij[i][3].d
328 + v * (cij[i][4].d
329 + v * (cij[i][5].d
330 + v * cij[i][6].d))));
331 t1 = hpi.d + cij[i][1].d;
332 z = t1 + zz;
333 /* Max ULP is 0.503. */
334 return signArctan2 (y, z);
335 }
336
337 /* (iv) x<0, abs(y)<=abs(x): pi-atan(ax/ay) */
338 if (u < inv16.d)
339 {
340 v = u * u;
341 zz = u * v * (d3.d
342 + v * (d5.d
343 + v * (d7.d
344 + v * (d9.d + v * (d11.d + v * d13.d)))));
345 ESUB (opi.d, u, t2, cor);
346 t3 = ((opi1.d + cor) - du) - zz;
347 z = t2 + t3;
348 /* Max ULP is 0.501. */
349 return signArctan2 (y, z);
350 }
351
352 i = (TWO52 + 256 * u) - TWO52;
353 i -= 16;
354 v = (u - cij[i][0].d) + du;
355 zz = opi1.d - v * (cij[i][2].d
356 + v * (cij[i][3].d
357 + v * (cij[i][4].d
358 + v * (cij[i][5].d + v * cij[i][6].d))));
359 t1 = opi.d - cij[i][1].d;
360 z = t1 + zz;
361 /* Max ULP is 0.502. */
362 return signArctan2 (y, z);
363 }
364
365 #ifndef __ieee754_atan2
366 libm_alias_finite (__ieee754_atan2, __atan2)
367 #endif