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