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
51void __mpatan (mp_no *, mp_no *, int); /* see definition in mpatan.c */
52static double atanMp (double, const int[]);
53
54 /* Fix the sign of y and return */
55static 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). */
64double
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 */
299static double
300atanMp (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
328libm_alias_double (__atan, atan)
329#endif
330