2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2015 Free Software Foundation, Inc.
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.
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.
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/>.
19 /************************************************************************/
20 /* MODULE_NAME: atnat.c */
22 /* FUNCTIONS: uatan */
27 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h atnat.h */
28 /* mpatan.c mpatan2.c mpsqrt.c */
31 /* An ultimate atan() routine. Given an IEEE double machine number x */
32 /* it computes the correctly rounded (to nearest) value of atan(x). */
34 /* Assumption: Machine arithmetic operations are performed in */
35 /* round to nearest mode of IEEE 754 standard. */
37 /************************************************************************/
47 #include <math_private.h>
48 #include <stap-probe.h>
50 void __mpatan (mp_no
*, mp_no
*, int); /* see definition in mpatan.c */
51 static double atanMp (double, const int[]);
53 /* Fix the sign of y and return */
55 __signArctan (double x
, double y
)
57 return __copysign (y
, x
);
61 /* An ultimate atan() routine. Given an IEEE double machine number x, */
62 /* routine computes the correctly rounded (to nearest) value of atan(x). */
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
;
72 static const int pr
[M
] = { 6, 8, 10, 32 };
76 ux
= num
.i
[HIGH_HALF
];
80 if (((ux
& 0x7ff00000) == 0x7ff00000)
81 && (((ux
& 0x000fffff) | dx
) != 0x00000000))
84 /* Regular values of x, including denormals +-0 and +-INF */
85 SET_RESTORE_ROUND (FE_TONEAREST
);
95 double force_underflow
= x
* x
;
96 math_force_eval (force_underflow
);
103 yy
= d11
.d
+ v
* d13
.d
;
110 if ((y
= x
+ (yy
- U1
* x
)) == x
+ (yy
+ U1
* x
))
113 EMULV (x
, x
, v
, vv
, t1
, t2
, t3
, t4
, t5
); /* v+vv=x^2 */
115 s1
= f17
.d
+ v
* f19
.d
;
121 ADD2 (f9
.d
, ff9
.d
, s1
, 0, s2
, ss2
, t1
, t2
);
122 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
123 ADD2 (f7
.d
, ff7
.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 (f5
.d
, ff5
.d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
126 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
127 ADD2 (f3
.d
, ff3
.d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
128 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
129 MUL2 (x
, 0, s1
, ss1
, s2
, ss2
, t1
, t2
, t3
, t4
, t5
, t6
, t7
,
131 ADD2 (x
, 0, s2
, ss2
, s1
, ss1
, t1
, t2
);
132 if ((y
= s1
+ (ss1
- U5
* s1
)) == s1
+ (ss1
+ U5
* s1
))
135 return atanMp (x
, pr
);
140 i
= (TWO52
+ TWO8
* u
) - TWO52
;
143 yy
= cij
[i
][5].d
+ z
* cij
[i
][6].d
;
144 yy
= cij
[i
][4].d
+ z
* yy
;
145 yy
= cij
[i
][3].d
+ z
* yy
;
146 yy
= cij
[i
][2].d
+ z
* yy
;
153 u2
= U21
; /* u < 1/4 */
156 } /* 1/4 <= u < 1/2 */
160 u2
= U23
; /* 1/2 <= u < 3/4 */
163 } /* 3/4 <= u <= 1 */
164 if ((y
= t1
+ (yy
- u2
* t1
)) == t1
+ (yy
+ u2
* t1
))
165 return __signArctan (x
, y
);
169 s1
= hij
[i
][14].d
+ z
* hij
[i
][15].d
;
170 s1
= hij
[i
][13].d
+ z
* s1
;
171 s1
= hij
[i
][12].d
+ z
* s1
;
172 s1
= hij
[i
][11].d
+ z
* s1
;
175 ADD2 (hij
[i
][9].d
, hij
[i
][10].d
, s1
, 0, s2
, ss2
, t1
, t2
);
176 MUL2 (z
, 0, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
177 ADD2 (hij
[i
][7].d
, hij
[i
][8].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
][5].d
, hij
[i
][6].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
180 MUL2 (z
, 0, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
181 ADD2 (hij
[i
][3].d
, hij
[i
][4].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
182 MUL2 (z
, 0, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
183 ADD2 (hij
[i
][1].d
, hij
[i
][2].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
184 if ((y
= s2
+ (ss2
- U6
* s2
)) == s2
+ (ss2
+ U6
* s2
))
185 return __signArctan (x
, y
);
187 return atanMp (x
, pr
);
195 EMULV (w
, u
, t1
, t2
, t3
, t4
, t5
, t6
, t7
);
196 ww
= w
* ((1 - t1
) - t2
);
197 i
= (TWO52
+ TWO8
* w
) - TWO52
;
199 z
= (w
- cij
[i
][0].d
) + ww
;
201 yy
= cij
[i
][5].d
+ z
* cij
[i
][6].d
;
202 yy
= cij
[i
][4].d
+ z
* yy
;
203 yy
= cij
[i
][3].d
+ z
* yy
;
204 yy
= cij
[i
][2].d
+ z
* yy
;
207 t1
= HPI
- cij
[i
][1].d
;
209 u3
= U31
; /* w < 1/2 */
211 u3
= U32
; /* w >= 1/2 */
212 if ((y
= t1
+ (yy
- u3
)) == t1
+ (yy
+ u3
))
213 return __signArctan (x
, y
);
215 DIV2 (1, 0, u
, 0, w
, ww
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
, t9
,
217 t1
= w
- hij
[i
][0].d
;
218 EADD (t1
, ww
, z
, zz
);
220 s1
= hij
[i
][14].d
+ z
* hij
[i
][15].d
;
221 s1
= hij
[i
][13].d
+ z
* s1
;
222 s1
= hij
[i
][12].d
+ z
* s1
;
223 s1
= hij
[i
][11].d
+ z
* s1
;
226 ADD2 (hij
[i
][9].d
, hij
[i
][10].d
, s1
, 0, s2
, ss2
, t1
, t2
);
227 MUL2 (z
, zz
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
228 ADD2 (hij
[i
][7].d
, hij
[i
][8].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
][5].d
, hij
[i
][6].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
][3].d
, hij
[i
][4].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
233 MUL2 (z
, zz
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
234 ADD2 (hij
[i
][1].d
, hij
[i
][2].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
235 SUB2 (HPI
, HPI1
, s2
, ss2
, s1
, ss1
, t1
, t2
);
236 if ((y
= s1
+ (ss1
- U7
)) == s1
+ (ss1
+ U7
))
237 return __signArctan (x
, y
);
239 return atanMp (x
, pr
);
247 EMULV (w
, u
, t1
, t2
, t3
, t4
, t5
, t6
, t7
);
249 yy
= d11
.d
+ v
* d13
.d
;
256 ww
= w
* ((1 - t1
) - t2
);
257 ESUB (HPI
, w
, t3
, cor
);
258 yy
= ((HPI1
+ cor
) - ww
) - yy
;
259 if ((y
= t3
+ (yy
- U4
)) == t3
+ (yy
+ U4
))
260 return __signArctan (x
, y
);
262 DIV2 (1, 0, u
, 0, w
, ww
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
,
264 MUL2 (w
, ww
, w
, ww
, v
, vv
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
266 s1
= f17
.d
+ v
* f19
.d
;
272 ADD2 (f9
.d
, ff9
.d
, s1
, 0, s2
, ss2
, t1
, t2
);
273 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
274 ADD2 (f7
.d
, ff7
.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 (f5
.d
, ff5
.d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
277 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
278 ADD2 (f3
.d
, ff3
.d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
279 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
280 MUL2 (w
, ww
, s1
, ss1
, s2
, ss2
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
281 ADD2 (w
, ww
, s2
, ss2
, s1
, ss1
, t1
, t2
);
282 SUB2 (HPI
, HPI1
, s1
, ss1
, s2
, ss2
, t1
, t2
);
284 if ((y
= s2
+ (ss2
- U8
)) == s2
+ (ss2
+ U8
))
285 return __signArctan (x
, y
);
287 return atanMp (x
, pr
);
301 /* Final stages. Compute atan(x) by multiple precision arithmetic */
303 atanMp (double x
, const int pr
[])
305 mp_no mpx
, mpy
, mpy2
, mperr
, mpt1
, mpy1
;
309 for (i
= 0; i
< M
; i
++)
312 __dbl_mp (x
, &mpx
, p
);
313 __mpatan (&mpx
, &mpy
, p
);
314 __dbl_mp (u9
[i
].d
, &mpt1
, p
);
315 __mul (&mpy
, &mpt1
, &mperr
, p
);
316 __add (&mpy
, &mperr
, &mpy1
, p
);
317 __sub (&mpy
, &mperr
, &mpy2
, p
);
318 __mp_dbl (&mpy1
, &y1
, p
);
319 __mp_dbl (&mpy2
, &y2
, p
);
322 LIBC_PROBE (slowatan
, 3, &p
, &x
, &y1
);
326 LIBC_PROBE (slowatan_inexact
, 3, &p
, &x
, &y1
);
327 return y1
; /*if impossible to do exact computing */
330 #ifdef NO_LONG_DOUBLE
331 weak_alias (atan
, atanl
)