1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2016 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
50void __mpatan (mp_no *, mp_no *, int); /* see definition in mpatan.c */
51static double atanMp (double, const int[]);
52
53 /* Fix the sign of y and return */
54static 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). */
63double
64atan (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 */
298static double
299atanMp (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
327weak_alias (atan, atanl)
328#endif
329