1 | /* |
2 | * IBM Accurate Mathematical Library |
3 | * written by International Business Machines Corp. |
4 | * Copyright (C) 2001-2018 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: atnat.c */ |
21 | /* */ |
22 | /* FUNCTIONS: uatan */ |
23 | /* atanMp */ |
24 | /* signArctan */ |
25 | /* */ |
26 | /* */ |
27 | /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h atnat.h */ |
28 | /* mpatan.c mpatan2.c mpsqrt.c */ |
29 | /* uatan.tbl */ |
30 | /* */ |
31 | /* An ultimate atan() routine. Given an IEEE double machine number x */ |
32 | /* it computes the correctly rounded (to nearest) value of atan(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 "atnat.h" |
44 | #include <fenv.h> |
45 | #include <float.h> |
46 | #include <libm-alias-double.h> |
47 | #include <math.h> |
48 | #include <math_private.h> |
49 | #include <stap-probe.h> |
50 | |
51 | void __mpatan (mp_no *, mp_no *, int); /* see definition in mpatan.c */ |
52 | static double atanMp (double, const int[]); |
53 | |
54 | /* Fix the sign of y and return */ |
55 | static double |
56 | __signArctan (double x, double y) |
57 | { |
58 | return __copysign (y, x); |
59 | } |
60 | |
61 | |
62 | /* An ultimate atan() routine. Given an IEEE double machine number x, */ |
63 | /* routine computes the correctly rounded (to nearest) value of atan(x). */ |
64 | double |
65 | __atan (double x) |
66 | { |
67 | double cor, s1, ss1, s2, ss2, t1, t2, t3, t7, t8, t9, t10, u, u2, u3, |
68 | v, vv, w, ww, y, yy, z, zz; |
69 | #ifndef DLA_FMS |
70 | double t4, t5, t6; |
71 | #endif |
72 | int i, ux, dx; |
73 | static const int pr[M] = { 6, 8, 10, 32 }; |
74 | number num; |
75 | |
76 | num.d = x; |
77 | ux = num.i[HIGH_HALF]; |
78 | dx = num.i[LOW_HALF]; |
79 | |
80 | /* x=NaN */ |
81 | if (((ux & 0x7ff00000) == 0x7ff00000) |
82 | && (((ux & 0x000fffff) | dx) != 0x00000000)) |
83 | return x + x; |
84 | |
85 | /* Regular values of x, including denormals +-0 and +-INF */ |
86 | SET_RESTORE_ROUND (FE_TONEAREST); |
87 | u = (x < 0) ? -x : x; |
88 | if (u < C) |
89 | { |
90 | if (u < B) |
91 | { |
92 | if (u < A) |
93 | { |
94 | math_check_force_underflow_nonneg (u); |
95 | return x; |
96 | } |
97 | else |
98 | { /* A <= u < B */ |
99 | v = x * x; |
100 | yy = d11.d + v * d13.d; |
101 | yy = d9.d + v * yy; |
102 | yy = d7.d + v * yy; |
103 | yy = d5.d + v * yy; |
104 | yy = d3.d + v * yy; |
105 | yy *= x * v; |
106 | |
107 | if ((y = x + (yy - U1 * x)) == x + (yy + U1 * x)) |
108 | return y; |
109 | |
110 | EMULV (x, x, v, vv, t1, t2, t3, t4, t5); /* v+vv=x^2 */ |
111 | |
112 | s1 = f17.d + v * f19.d; |
113 | s1 = f15.d + v * s1; |
114 | s1 = f13.d + v * s1; |
115 | s1 = f11.d + v * s1; |
116 | s1 *= v; |
117 | |
118 | ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2); |
119 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
120 | ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2); |
121 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
122 | ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2); |
123 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
124 | ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2); |
125 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
126 | MUL2 (x, 0, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, |
127 | t8); |
128 | ADD2 (x, 0, s2, ss2, s1, ss1, t1, t2); |
129 | if ((y = s1 + (ss1 - U5 * s1)) == s1 + (ss1 + U5 * s1)) |
130 | return y; |
131 | |
132 | return atanMp (x, pr); |
133 | } |
134 | } |
135 | else |
136 | { /* B <= u < C */ |
137 | i = (TWO52 + TWO8 * u) - TWO52; |
138 | i -= 16; |
139 | z = u - cij[i][0].d; |
140 | yy = cij[i][5].d + z * cij[i][6].d; |
141 | yy = cij[i][4].d + z * yy; |
142 | yy = cij[i][3].d + z * yy; |
143 | yy = cij[i][2].d + z * yy; |
144 | yy *= z; |
145 | |
146 | t1 = cij[i][1].d; |
147 | if (i < 112) |
148 | { |
149 | if (i < 48) |
150 | u2 = U21; /* u < 1/4 */ |
151 | else |
152 | u2 = U22; |
153 | } /* 1/4 <= u < 1/2 */ |
154 | else |
155 | { |
156 | if (i < 176) |
157 | u2 = U23; /* 1/2 <= u < 3/4 */ |
158 | else |
159 | u2 = U24; |
160 | } /* 3/4 <= u <= 1 */ |
161 | if ((y = t1 + (yy - u2 * t1)) == t1 + (yy + u2 * t1)) |
162 | return __signArctan (x, y); |
163 | |
164 | z = u - hij[i][0].d; |
165 | |
166 | s1 = hij[i][14].d + z * hij[i][15].d; |
167 | s1 = hij[i][13].d + z * s1; |
168 | s1 = hij[i][12].d + z * s1; |
169 | s1 = hij[i][11].d + z * s1; |
170 | s1 *= z; |
171 | |
172 | ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2); |
173 | MUL2 (z, 0, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
174 | ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2); |
175 | MUL2 (z, 0, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
176 | ADD2 (hij[i][5].d, hij[i][6].d, s1, ss1, s2, ss2, t1, t2); |
177 | MUL2 (z, 0, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
178 | ADD2 (hij[i][3].d, hij[i][4].d, s1, ss1, s2, ss2, t1, t2); |
179 | MUL2 (z, 0, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
180 | ADD2 (hij[i][1].d, hij[i][2].d, s1, ss1, s2, ss2, t1, t2); |
181 | if ((y = s2 + (ss2 - U6 * s2)) == s2 + (ss2 + U6 * s2)) |
182 | return __signArctan (x, y); |
183 | |
184 | return atanMp (x, pr); |
185 | } |
186 | } |
187 | else |
188 | { |
189 | if (u < D) |
190 | { /* C <= u < D */ |
191 | w = 1 / u; |
192 | EMULV (w, u, t1, t2, t3, t4, t5, t6, t7); |
193 | ww = w * ((1 - t1) - t2); |
194 | i = (TWO52 + TWO8 * w) - TWO52; |
195 | i -= 16; |
196 | z = (w - cij[i][0].d) + ww; |
197 | |
198 | yy = cij[i][5].d + z * cij[i][6].d; |
199 | yy = cij[i][4].d + z * yy; |
200 | yy = cij[i][3].d + z * yy; |
201 | yy = cij[i][2].d + z * yy; |
202 | yy = HPI1 - z * yy; |
203 | |
204 | t1 = HPI - cij[i][1].d; |
205 | if (i < 112) |
206 | u3 = U31; /* w < 1/2 */ |
207 | else |
208 | u3 = U32; /* w >= 1/2 */ |
209 | if ((y = t1 + (yy - u3)) == t1 + (yy + u3)) |
210 | return __signArctan (x, y); |
211 | |
212 | DIV2 (1, 0, u, 0, w, ww, t1, t2, t3, t4, t5, t6, t7, t8, t9, |
213 | t10); |
214 | t1 = w - hij[i][0].d; |
215 | EADD (t1, ww, z, zz); |
216 | |
217 | s1 = hij[i][14].d + z * hij[i][15].d; |
218 | s1 = hij[i][13].d + z * s1; |
219 | s1 = hij[i][12].d + z * s1; |
220 | s1 = hij[i][11].d + z * s1; |
221 | s1 *= z; |
222 | |
223 | ADD2 (hij[i][9].d, hij[i][10].d, s1, 0, s2, ss2, t1, t2); |
224 | MUL2 (z, zz, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
225 | ADD2 (hij[i][7].d, hij[i][8].d, s1, ss1, s2, ss2, t1, t2); |
226 | MUL2 (z, zz, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
227 | ADD2 (hij[i][5].d, hij[i][6].d, s1, ss1, s2, ss2, t1, t2); |
228 | MUL2 (z, zz, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
229 | ADD2 (hij[i][3].d, hij[i][4].d, s1, ss1, s2, ss2, t1, t2); |
230 | MUL2 (z, zz, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
231 | ADD2 (hij[i][1].d, hij[i][2].d, s1, ss1, s2, ss2, t1, t2); |
232 | SUB2 (HPI, HPI1, s2, ss2, s1, ss1, t1, t2); |
233 | if ((y = s1 + (ss1 - U7)) == s1 + (ss1 + U7)) |
234 | return __signArctan (x, y); |
235 | |
236 | return atanMp (x, pr); |
237 | } |
238 | else |
239 | { |
240 | if (u < E) |
241 | { /* D <= u < E */ |
242 | w = 1 / u; |
243 | v = w * w; |
244 | EMULV (w, u, t1, t2, t3, t4, t5, t6, t7); |
245 | |
246 | yy = d11.d + v * d13.d; |
247 | yy = d9.d + v * yy; |
248 | yy = d7.d + v * yy; |
249 | yy = d5.d + v * yy; |
250 | yy = d3.d + v * yy; |
251 | yy *= w * v; |
252 | |
253 | ww = w * ((1 - t1) - t2); |
254 | ESUB (HPI, w, t3, cor); |
255 | yy = ((HPI1 + cor) - ww) - yy; |
256 | if ((y = t3 + (yy - U4)) == t3 + (yy + U4)) |
257 | return __signArctan (x, y); |
258 | |
259 | DIV2 (1, 0, u, 0, w, ww, t1, t2, t3, t4, t5, t6, t7, t8, |
260 | t9, t10); |
261 | MUL2 (w, ww, w, ww, v, vv, t1, t2, t3, t4, t5, t6, t7, t8); |
262 | |
263 | s1 = f17.d + v * f19.d; |
264 | s1 = f15.d + v * s1; |
265 | s1 = f13.d + v * s1; |
266 | s1 = f11.d + v * s1; |
267 | s1 *= v; |
268 | |
269 | ADD2 (f9.d, ff9.d, s1, 0, s2, ss2, t1, t2); |
270 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
271 | ADD2 (f7.d, ff7.d, s1, ss1, s2, ss2, t1, t2); |
272 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
273 | ADD2 (f5.d, ff5.d, s1, ss1, s2, ss2, t1, t2); |
274 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
275 | ADD2 (f3.d, ff3.d, s1, ss1, s2, ss2, t1, t2); |
276 | MUL2 (v, vv, s2, ss2, s1, ss1, t1, t2, t3, t4, t5, t6, t7, t8); |
277 | MUL2 (w, ww, s1, ss1, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8); |
278 | ADD2 (w, ww, s2, ss2, s1, ss1, t1, t2); |
279 | SUB2 (HPI, HPI1, s1, ss1, s2, ss2, t1, t2); |
280 | |
281 | if ((y = s2 + (ss2 - U8)) == s2 + (ss2 + U8)) |
282 | return __signArctan (x, y); |
283 | |
284 | return atanMp (x, pr); |
285 | } |
286 | else |
287 | { |
288 | /* u >= E */ |
289 | if (x > 0) |
290 | return HPI; |
291 | else |
292 | return MHPI; |
293 | } |
294 | } |
295 | } |
296 | } |
297 | |
298 | /* Final stages. Compute atan(x) by multiple precision arithmetic */ |
299 | static double |
300 | atanMp (double x, const int pr[]) |
301 | { |
302 | mp_no mpx, mpy, mpy2, mperr, mpt1, mpy1; |
303 | double y1, y2; |
304 | int i, p; |
305 | |
306 | for (i = 0; i < M; i++) |
307 | { |
308 | p = pr[i]; |
309 | __dbl_mp (x, &mpx, p); |
310 | __mpatan (&mpx, &mpy, p); |
311 | __dbl_mp (u9[i].d, &mpt1, p); |
312 | __mul (&mpy, &mpt1, &mperr, p); |
313 | __add (&mpy, &mperr, &mpy1, p); |
314 | __sub (&mpy, &mperr, &mpy2, p); |
315 | __mp_dbl (&mpy1, &y1, p); |
316 | __mp_dbl (&mpy2, &y2, p); |
317 | if (y1 == y2) |
318 | { |
319 | LIBC_PROBE (slowatan, 3, &p, &x, &y1); |
320 | return y1; |
321 | } |
322 | } |
323 | LIBC_PROBE (slowatan_inexact, 3, &p, &x, &y1); |
324 | return y1; /*if impossible to do exact computing */ |
325 | } |
326 | |
327 | #ifndef __atan |
328 | libm_alias_double (__atan, atan) |
329 | #endif |
330 | |