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