2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2014 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 /************************************************************************/
45 #include <stap-probe.h>
47 void __mpatan (mp_no
*, mp_no
*, int); /* see definition in mpatan.c */
48 static double atanMp (double, const int[]);
50 /* Fix the sign of y and return */
52 __signArctan (double x
, double y
)
54 return __copysign (y
, x
);
58 /* An ultimate atan() routine. Given an IEEE double machine number x, */
59 /* routine computes the correctly rounded (to nearest) value of atan(x). */
63 double cor
, s1
, ss1
, s2
, ss2
, t1
, t2
, t3
, t7
, t8
, t9
, t10
, u
, u2
, u3
,
64 v
, vv
, w
, ww
, y
, yy
, z
, zz
;
69 static const int pr
[M
] = { 6, 8, 10, 32 };
73 ux
= num
.i
[HIGH_HALF
];
77 if (((ux
& 0x7ff00000) == 0x7ff00000)
78 && (((ux
& 0x000fffff) | dx
) != 0x00000000))
81 /* Regular values of x, including denormals +-0 and +-INF */
92 yy
= d11
.d
+ v
* d13
.d
;
99 if ((y
= x
+ (yy
- U1
* x
)) == x
+ (yy
+ U1
* x
))
102 EMULV (x
, x
, v
, vv
, t1
, t2
, t3
, t4
, t5
); /* v+vv=x^2 */
104 s1
= f17
.d
+ v
* f19
.d
;
110 ADD2 (f9
.d
, ff9
.d
, s1
, 0, s2
, ss2
, t1
, t2
);
111 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
112 ADD2 (f7
.d
, ff7
.d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
113 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
114 ADD2 (f5
.d
, ff5
.d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
115 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
116 ADD2 (f3
.d
, ff3
.d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
117 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
118 MUL2 (x
, 0, s1
, ss1
, s2
, ss2
, t1
, t2
, t3
, t4
, t5
, t6
, t7
,
120 ADD2 (x
, 0, s2
, ss2
, s1
, ss1
, t1
, t2
);
121 if ((y
= s1
+ (ss1
- U5
* s1
)) == s1
+ (ss1
+ U5
* s1
))
124 return atanMp (x
, pr
);
129 i
= (TWO52
+ TWO8
* u
) - TWO52
;
132 yy
= cij
[i
][5].d
+ z
* cij
[i
][6].d
;
133 yy
= cij
[i
][4].d
+ z
* yy
;
134 yy
= cij
[i
][3].d
+ z
* yy
;
135 yy
= cij
[i
][2].d
+ z
* yy
;
142 u2
= U21
; /* u < 1/4 */
145 } /* 1/4 <= u < 1/2 */
149 u2
= U23
; /* 1/2 <= u < 3/4 */
152 } /* 3/4 <= u <= 1 */
153 if ((y
= t1
+ (yy
- u2
* t1
)) == t1
+ (yy
+ u2
* t1
))
154 return __signArctan (x
, y
);
158 s1
= hij
[i
][14].d
+ z
* hij
[i
][15].d
;
159 s1
= hij
[i
][13].d
+ z
* s1
;
160 s1
= hij
[i
][12].d
+ z
* s1
;
161 s1
= hij
[i
][11].d
+ z
* s1
;
164 ADD2 (hij
[i
][9].d
, hij
[i
][10].d
, s1
, 0, s2
, ss2
, t1
, t2
);
165 MUL2 (z
, 0, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
166 ADD2 (hij
[i
][7].d
, hij
[i
][8].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
167 MUL2 (z
, 0, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
168 ADD2 (hij
[i
][5].d
, hij
[i
][6].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
169 MUL2 (z
, 0, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
170 ADD2 (hij
[i
][3].d
, hij
[i
][4].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
171 MUL2 (z
, 0, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
172 ADD2 (hij
[i
][1].d
, hij
[i
][2].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
173 if ((y
= s2
+ (ss2
- U6
* s2
)) == s2
+ (ss2
+ U6
* s2
))
174 return __signArctan (x
, y
);
176 return atanMp (x
, pr
);
184 EMULV (w
, u
, t1
, t2
, t3
, t4
, t5
, t6
, t7
);
185 ww
= w
* ((1 - t1
) - t2
);
186 i
= (TWO52
+ TWO8
* w
) - TWO52
;
188 z
= (w
- cij
[i
][0].d
) + ww
;
190 yy
= cij
[i
][5].d
+ z
* cij
[i
][6].d
;
191 yy
= cij
[i
][4].d
+ z
* yy
;
192 yy
= cij
[i
][3].d
+ z
* yy
;
193 yy
= cij
[i
][2].d
+ z
* yy
;
196 t1
= HPI
- cij
[i
][1].d
;
198 u3
= U31
; /* w < 1/2 */
200 u3
= U32
; /* w >= 1/2 */
201 if ((y
= t1
+ (yy
- u3
)) == t1
+ (yy
+ u3
))
202 return __signArctan (x
, y
);
204 DIV2 (1, 0, u
, 0, w
, ww
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
, t9
,
206 t1
= w
- hij
[i
][0].d
;
207 EADD (t1
, ww
, z
, zz
);
209 s1
= hij
[i
][14].d
+ z
* hij
[i
][15].d
;
210 s1
= hij
[i
][13].d
+ z
* s1
;
211 s1
= hij
[i
][12].d
+ z
* s1
;
212 s1
= hij
[i
][11].d
+ z
* s1
;
215 ADD2 (hij
[i
][9].d
, hij
[i
][10].d
, s1
, 0, s2
, ss2
, t1
, t2
);
216 MUL2 (z
, zz
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
217 ADD2 (hij
[i
][7].d
, hij
[i
][8].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
218 MUL2 (z
, zz
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
219 ADD2 (hij
[i
][5].d
, hij
[i
][6].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
220 MUL2 (z
, zz
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
221 ADD2 (hij
[i
][3].d
, hij
[i
][4].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
222 MUL2 (z
, zz
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
223 ADD2 (hij
[i
][1].d
, hij
[i
][2].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
224 SUB2 (HPI
, HPI1
, s2
, ss2
, s1
, ss1
, t1
, t2
);
225 if ((y
= s1
+ (ss1
- U7
)) == s1
+ (ss1
+ U7
))
226 return __signArctan (x
, y
);
228 return atanMp (x
, pr
);
236 EMULV (w
, u
, t1
, t2
, t3
, t4
, t5
, t6
, t7
);
238 yy
= d11
.d
+ v
* d13
.d
;
245 ww
= w
* ((1 - t1
) - t2
);
246 ESUB (HPI
, w
, t3
, cor
);
247 yy
= ((HPI1
+ cor
) - ww
) - yy
;
248 if ((y
= t3
+ (yy
- U4
)) == t3
+ (yy
+ U4
))
249 return __signArctan (x
, y
);
251 DIV2 (1, 0, u
, 0, w
, ww
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
,
253 MUL2 (w
, ww
, w
, ww
, v
, vv
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
255 s1
= f17
.d
+ v
* f19
.d
;
261 ADD2 (f9
.d
, ff9
.d
, s1
, 0, s2
, ss2
, t1
, t2
);
262 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
263 ADD2 (f7
.d
, ff7
.d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
264 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
265 ADD2 (f5
.d
, ff5
.d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
266 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
267 ADD2 (f3
.d
, ff3
.d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
268 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
269 MUL2 (w
, ww
, s1
, ss1
, s2
, ss2
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
270 ADD2 (w
, ww
, s2
, ss2
, s1
, ss1
, t1
, t2
);
271 SUB2 (HPI
, HPI1
, s1
, ss1
, s2
, ss2
, t1
, t2
);
273 if ((y
= s2
+ (ss2
- U8
)) == s2
+ (ss2
+ U8
))
274 return __signArctan (x
, y
);
276 return atanMp (x
, pr
);
290 /* Final stages. Compute atan(x) by multiple precision arithmetic */
292 atanMp (double x
, const int pr
[])
294 mp_no mpx
, mpy
, mpy2
, mperr
, mpt1
, mpy1
;
298 for (i
= 0; i
< M
; i
++)
301 __dbl_mp (x
, &mpx
, p
);
302 __mpatan (&mpx
, &mpy
, p
);
303 __dbl_mp (u9
[i
].d
, &mpt1
, p
);
304 __mul (&mpy
, &mpt1
, &mperr
, p
);
305 __add (&mpy
, &mperr
, &mpy1
, p
);
306 __sub (&mpy
, &mperr
, &mpy2
, p
);
307 __mp_dbl (&mpy1
, &y1
, p
);
308 __mp_dbl (&mpy2
, &y2
, p
);
311 LIBC_PROBE (slowatan
, 3, &p
, &x
, &y1
);
315 LIBC_PROBE (slowatan_inexact
, 3, &p
, &x
, &y1
);
316 return y1
; /*if impossible to do exact computing */
319 #ifdef NO_LONG_DOUBLE
320 weak_alias (atan
, atanl
)