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