]>
git.ipfire.org Git - thirdparty/glibc.git/blob - sysdeps/ieee754/dbl-64/s_atan.c
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2013 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 /************************************************************************/
46 void __mpatan (mp_no
*, mp_no
*, int); /* see definition in mpatan.c */
47 static double atanMp (double, const int[]);
49 /* Fix the sign of y and return */
51 __signArctan (double x
, double y
)
53 return __copysign (y
, x
);
57 /* An ultimate atan() routine. Given an IEEE double machine number x, */
58 /* routine computes the correctly rounded (to nearest) value of atan(x). */
62 double cor
, s1
, ss1
, s2
, ss2
, t1
, t2
, t3
, t7
, t8
, t9
, t10
, u
, u2
, u3
,
63 v
, vv
, w
, ww
, y
, yy
, z
, zz
;
68 static const int pr
[M
] = { 6, 8, 10, 32 };
72 ux
= num
.i
[HIGH_HALF
];
76 if (((ux
& 0x7ff00000) == 0x7ff00000)
77 && (((ux
& 0x000fffff) | dx
) != 0x00000000))
80 /* Regular values of x, including denormals +-0 and +-INF */
91 yy
= d11
.d
+ v
* d13
.d
;
98 if ((y
= x
+ (yy
- U1
* x
)) == x
+ (yy
+ U1
* x
))
101 EMULV (x
, x
, v
, vv
, t1
, t2
, t3
, t4
, t5
); /* v+vv=x^2 */
103 s1
= f17
.d
+ v
* f19
.d
;
109 ADD2 (f9
.d
, ff9
.d
, s1
, 0, s2
, ss2
, t1
, t2
);
110 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
111 ADD2 (f7
.d
, ff7
.d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
112 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
113 ADD2 (f5
.d
, ff5
.d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
114 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
115 ADD2 (f3
.d
, ff3
.d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
116 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
117 MUL2 (x
, 0, s1
, ss1
, s2
, ss2
, t1
, t2
, t3
, t4
, t5
, t6
, t7
,
119 ADD2 (x
, 0, s2
, ss2
, s1
, ss1
, t1
, t2
);
120 if ((y
= s1
+ (ss1
- U5
* s1
)) == s1
+ (ss1
+ U5
* s1
))
123 return atanMp (x
, pr
);
128 i
= (TWO52
+ TWO8
* u
) - TWO52
;
131 yy
= cij
[i
][5].d
+ z
* cij
[i
][6].d
;
132 yy
= cij
[i
][4].d
+ z
* yy
;
133 yy
= cij
[i
][3].d
+ z
* yy
;
134 yy
= cij
[i
][2].d
+ z
* yy
;
141 u2
= U21
; /* u < 1/4 */
144 } /* 1/4 <= u < 1/2 */
148 u2
= U23
; /* 1/2 <= u < 3/4 */
151 } /* 3/4 <= u <= 1 */
152 if ((y
= t1
+ (yy
- u2
* t1
)) == t1
+ (yy
+ u2
* t1
))
153 return __signArctan (x
, y
);
157 s1
= hij
[i
][14].d
+ z
* hij
[i
][15].d
;
158 s1
= hij
[i
][13].d
+ z
* s1
;
159 s1
= hij
[i
][12].d
+ z
* s1
;
160 s1
= hij
[i
][11].d
+ z
* s1
;
163 ADD2 (hij
[i
][9].d
, hij
[i
][10].d
, s1
, 0, s2
, ss2
, t1
, t2
);
164 MUL2 (z
, 0, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
165 ADD2 (hij
[i
][7].d
, hij
[i
][8].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
166 MUL2 (z
, 0, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
167 ADD2 (hij
[i
][5].d
, hij
[i
][6].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
168 MUL2 (z
, 0, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
169 ADD2 (hij
[i
][3].d
, hij
[i
][4].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
170 MUL2 (z
, 0, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
171 ADD2 (hij
[i
][1].d
, hij
[i
][2].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
172 if ((y
= s2
+ (ss2
- U6
* s2
)) == s2
+ (ss2
+ U6
* s2
))
173 return __signArctan (x
, y
);
175 return atanMp (x
, pr
);
183 EMULV (w
, u
, t1
, t2
, t3
, t4
, t5
, t6
, t7
);
184 ww
= w
* ((ONE
- t1
) - t2
);
185 i
= (TWO52
+ TWO8
* w
) - TWO52
;
187 z
= (w
- cij
[i
][0].d
) + ww
;
189 yy
= cij
[i
][5].d
+ z
* cij
[i
][6].d
;
190 yy
= cij
[i
][4].d
+ z
* yy
;
191 yy
= cij
[i
][3].d
+ z
* yy
;
192 yy
= cij
[i
][2].d
+ z
* yy
;
195 t1
= HPI
- cij
[i
][1].d
;
197 u3
= U31
; /* w < 1/2 */
199 u3
= U32
; /* w >= 1/2 */
200 if ((y
= t1
+ (yy
- u3
)) == t1
+ (yy
+ u3
))
201 return __signArctan (x
, y
);
203 DIV2 (ONE
, 0, u
, 0, w
, ww
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
, t9
,
205 t1
= w
- hij
[i
][0].d
;
206 EADD (t1
, ww
, z
, zz
);
208 s1
= hij
[i
][14].d
+ z
* hij
[i
][15].d
;
209 s1
= hij
[i
][13].d
+ z
* s1
;
210 s1
= hij
[i
][12].d
+ z
* s1
;
211 s1
= hij
[i
][11].d
+ z
* s1
;
214 ADD2 (hij
[i
][9].d
, hij
[i
][10].d
, s1
, 0, s2
, ss2
, t1
, t2
);
215 MUL2 (z
, zz
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
216 ADD2 (hij
[i
][7].d
, hij
[i
][8].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
217 MUL2 (z
, zz
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
218 ADD2 (hij
[i
][5].d
, hij
[i
][6].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
219 MUL2 (z
, zz
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
220 ADD2 (hij
[i
][3].d
, hij
[i
][4].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
221 MUL2 (z
, zz
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
222 ADD2 (hij
[i
][1].d
, hij
[i
][2].d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
223 SUB2 (HPI
, HPI1
, s2
, ss2
, s1
, ss1
, t1
, t2
);
224 if ((y
= s1
+ (ss1
- U7
)) == s1
+ (ss1
+ U7
))
225 return __signArctan (x
, y
);
227 return atanMp (x
, pr
);
235 EMULV (w
, u
, t1
, t2
, t3
, t4
, t5
, t6
, t7
);
237 yy
= d11
.d
+ v
* d13
.d
;
244 ww
= w
* ((ONE
- t1
) - t2
);
245 ESUB (HPI
, w
, t3
, cor
);
246 yy
= ((HPI1
+ cor
) - ww
) - yy
;
247 if ((y
= t3
+ (yy
- U4
)) == t3
+ (yy
+ U4
))
248 return __signArctan (x
, y
);
250 DIV2 (ONE
, 0, u
, 0, w
, ww
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
,
252 MUL2 (w
, ww
, w
, ww
, v
, vv
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
254 s1
= f17
.d
+ v
* f19
.d
;
260 ADD2 (f9
.d
, ff9
.d
, s1
, 0, s2
, ss2
, t1
, t2
);
261 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
262 ADD2 (f7
.d
, ff7
.d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
263 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
264 ADD2 (f5
.d
, ff5
.d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
265 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
266 ADD2 (f3
.d
, ff3
.d
, s1
, ss1
, s2
, ss2
, t1
, t2
);
267 MUL2 (v
, vv
, s2
, ss2
, s1
, ss1
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
268 MUL2 (w
, ww
, s1
, ss1
, s2
, ss2
, t1
, t2
, t3
, t4
, t5
, t6
, t7
, t8
);
269 ADD2 (w
, ww
, s2
, ss2
, s1
, ss1
, t1
, t2
);
270 SUB2 (HPI
, HPI1
, s1
, ss1
, s2
, ss2
, t1
, t2
);
272 if ((y
= s2
+ (ss2
- U8
)) == s2
+ (ss2
+ U8
))
273 return __signArctan (x
, y
);
275 return atanMp (x
, pr
);
289 /* Final stages. Compute atan(x) by multiple precision arithmetic */
291 atanMp (double x
, const int pr
[])
293 mp_no mpx
, mpy
, mpy2
, mperr
, mpt1
, mpy1
;
297 for (i
= 0; i
< M
; i
++)
300 __dbl_mp (x
, &mpx
, p
);
301 __mpatan (&mpx
, &mpy
, p
);
302 __dbl_mp (u9
[i
].d
, &mpt1
, p
);
303 __mul (&mpy
, &mpt1
, &mperr
, p
);
304 __add (&mpy
, &mperr
, &mpy1
, p
);
305 __sub (&mpy
, &mperr
, &mpy2
, p
);
306 __mp_dbl (&mpy1
, &y1
, p
);
307 __mp_dbl (&mpy2
, &y2
, p
);
311 return y1
; /*if impossible to do exact computing */
314 #ifdef NO_LONG_DOUBLE
315 weak_alias (atan
, atanl
)