]>
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, 2011 Free Software Foundation
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, write to the Free Software
18 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
20 /************************************************************************/
21 /* MODULE_NAME: atnat.c */
23 /* FUNCTIONS: uatan */
28 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h atnat.h */
29 /* mpatan.c mpatan2.c mpsqrt.c */
32 /* An ultimate atan() routine. Given an IEEE double machine number x */
33 /* it computes the correctly rounded (to nearest) value of atan(x). */
35 /* Assumption: Machine arithmetic operations are performed in */
36 /* round to nearest mode of IEEE 754 standard. */
38 /************************************************************************/
47 void __mpatan(mp_no
*,mp_no
*,int); /* see definition in mpatan.c */
48 static double atanMp(double,const int[]);
49 double __signArctan(double,double);
50 /* An ultimate atan() routine. Given an IEEE double machine number x, */
51 /* routine computes the correctly rounded (to nearest) value of atan(x). */
52 double atan(double x
) {
55 double cor
,s1
,ss1
,s2
,ss2
,t1
,t2
,t3
,t7
,t8
,t9
,t10
,u
,u2
,u3
,
67 static const int pr
[M
]={6,8,10,32};
70 mp_no mpt1
,mpx
,mpy
,mpy1
,mpy2
,mperr
;
73 num
.d
= x
; ux
= num
.i
[HIGH_HALF
]; dx
= num
.i
[LOW_HALF
];
76 if (((ux
&0x7ff00000)==0x7ff00000) && (((ux
&0x000fffff)|dx
)!=0x00000000))
79 /* Regular values of x, including denormals +-0 and +-INF */
80 u
= (x
<ZERO
) ? -x
: x
;
83 if (u
<A
) { /* u < A */
85 else { /* A <= u < B */
86 v
=x
*x
; yy
=x
*v
*(d3
.d
+v
*(d5
.d
+v
*(d7
.d
+v
*(d9
.d
+v
*(d11
.d
+v
*d13
.d
)))));
87 if ((y
=x
+(yy
-U1
*x
)) == x
+(yy
+U1
*x
)) return y
;
89 EMULV(x
,x
,v
,vv
,t1
,t2
,t3
,t4
,t5
) /* v+vv=x^2 */
90 s1
=v
*(f11
.d
+v
*(f13
.d
+v
*(f15
.d
+v
*(f17
.d
+v
*f19
.d
))));
91 ADD2(f9
.d
,ff9
.d
,s1
,ZERO
,s2
,ss2
,t1
,t2
)
92 MUL2(v
,vv
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
93 ADD2(f7
.d
,ff7
.d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
94 MUL2(v
,vv
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
95 ADD2(f5
.d
,ff5
.d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
96 MUL2(v
,vv
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
97 ADD2(f3
.d
,ff3
.d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
98 MUL2(v
,vv
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
99 MUL2(x
,ZERO
,s1
,ss1
,s2
,ss2
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
100 ADD2(x
,ZERO
,s2
,ss2
,s1
,ss1
,t1
,t2
)
101 if ((y
=s1
+(ss1
-U5
*s1
)) == s1
+(ss1
+U5
*s1
)) return y
;
105 else { /* B <= u < C */
106 i
=(TWO52
+TWO8
*u
)-TWO52
; i
-=16;
108 yy
=z
*(cij
[i
][2].d
+z
*(cij
[i
][3].d
+z
*(cij
[i
][4].d
+
109 z
*(cij
[i
][5].d
+z
* cij
[i
][6].d
))));
112 if (i
<48) u2
=U21
; /* u < 1/4 */
113 else u2
=U22
; } /* 1/4 <= u < 1/2 */
115 if (i
<176) u2
=U23
; /* 1/2 <= u < 3/4 */
116 else u2
=U24
; } /* 3/4 <= u <= 1 */
117 if ((y
=t1
+(yy
-u2
*t1
)) == t1
+(yy
+u2
*t1
)) return __signArctan(x
,y
);
120 s1
=z
*(hij
[i
][11].d
+z
*(hij
[i
][12].d
+z
*(hij
[i
][13].d
+
121 z
*(hij
[i
][14].d
+z
* hij
[i
][15].d
))));
122 ADD2(hij
[i
][9].d
,hij
[i
][10].d
,s1
,ZERO
,s2
,ss2
,t1
,t2
)
123 MUL2(z
,ZERO
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
124 ADD2(hij
[i
][7].d
,hij
[i
][8].d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
125 MUL2(z
,ZERO
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
126 ADD2(hij
[i
][5].d
,hij
[i
][6].d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
127 MUL2(z
,ZERO
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
128 ADD2(hij
[i
][3].d
,hij
[i
][4].d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
129 MUL2(z
,ZERO
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
130 ADD2(hij
[i
][1].d
,hij
[i
][2].d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
131 if ((y
=s2
+(ss2
-U6
*s2
)) == s2
+(ss2
+U6
*s2
)) return __signArctan(x
,y
);
137 if (u
<D
) { /* C <= u < D */
139 EMULV(w
,u
,t1
,t2
,t3
,t4
,t5
,t6
,t7
)
141 i
=(TWO52
+TWO8
*w
)-TWO52
; i
-=16;
142 z
=(w
-cij
[i
][0].d
)+ww
;
143 yy
=HPI1
-z
*(cij
[i
][2].d
+z
*(cij
[i
][3].d
+z
*(cij
[i
][4].d
+
144 z
*(cij
[i
][5].d
+z
* cij
[i
][6].d
))));
146 if (i
<112) u3
=U31
; /* w < 1/2 */
147 else u3
=U32
; /* w >= 1/2 */
148 if ((y
=t1
+(yy
-u3
)) == t1
+(yy
+u3
)) return __signArctan(x
,y
);
150 DIV2(ONE
,ZERO
,u
,ZERO
,w
,ww
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
,t9
,t10
)
153 s1
=z
*(hij
[i
][11].d
+z
*(hij
[i
][12].d
+z
*(hij
[i
][13].d
+
154 z
*(hij
[i
][14].d
+z
* hij
[i
][15].d
))));
155 ADD2(hij
[i
][9].d
,hij
[i
][10].d
,s1
,ZERO
,s2
,ss2
,t1
,t2
)
156 MUL2(z
,zz
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
157 ADD2(hij
[i
][7].d
,hij
[i
][8].d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
158 MUL2(z
,zz
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
159 ADD2(hij
[i
][5].d
,hij
[i
][6].d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
160 MUL2(z
,zz
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
161 ADD2(hij
[i
][3].d
,hij
[i
][4].d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
162 MUL2(z
,zz
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
163 ADD2(hij
[i
][1].d
,hij
[i
][2].d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
164 SUB2(HPI
,HPI1
,s2
,ss2
,s1
,ss1
,t1
,t2
)
165 if ((y
=s1
+(ss1
-U7
)) == s1
+(ss1
+U7
)) return __signArctan(x
,y
);
170 if (u
<E
) { /* D <= u < E */
172 EMULV(w
,u
,t1
,t2
,t3
,t4
,t5
,t6
,t7
)
173 yy
=w
*v
*(d3
.d
+v
*(d5
.d
+v
*(d7
.d
+v
*(d9
.d
+v
*(d11
.d
+v
*d13
.d
)))));
176 yy
=((HPI1
+cor
)-ww
)-yy
;
177 if ((y
=t3
+(yy
-U4
)) == t3
+(yy
+U4
)) return __signArctan(x
,y
);
179 DIV2(ONE
,ZERO
,u
,ZERO
,w
,ww
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
,t9
,t10
)
180 MUL2(w
,ww
,w
,ww
,v
,vv
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
181 s1
=v
*(f11
.d
+v
*(f13
.d
+v
*(f15
.d
+v
*(f17
.d
+v
*f19
.d
))));
182 ADD2(f9
.d
,ff9
.d
,s1
,ZERO
,s2
,ss2
,t1
,t2
)
183 MUL2(v
,vv
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
184 ADD2(f7
.d
,ff7
.d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
185 MUL2(v
,vv
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
186 ADD2(f5
.d
,ff5
.d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
187 MUL2(v
,vv
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
188 ADD2(f3
.d
,ff3
.d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
189 MUL2(v
,vv
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
190 MUL2(w
,ww
,s1
,ss1
,s2
,ss2
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
191 ADD2(w
,ww
,s2
,ss2
,s1
,ss1
,t1
,t2
)
192 SUB2(HPI
,HPI1
,s1
,ss1
,s2
,ss2
,t1
,t2
)
193 if ((y
=s2
+(ss2
-U8
)) == s2
+(ss2
+U8
)) return __signArctan(x
,y
);
207 /* Fix the sign of y and return */
208 double __signArctan(double x
,double y
){
210 if (x
<ZERO
) return -y
;
214 /* Final stages. Compute atan(x) by multiple precision arithmetic */
215 static double atanMp(double x
,const int pr
[]){
216 mp_no mpx
,mpy
,mpy2
,mperr
,mpt1
,mpy1
;
220 for (i
=0; i
<M
; i
++) {
222 __dbl_mp(x
,&mpx
,p
); __mpatan(&mpx
,&mpy
,p
);
223 __dbl_mp(u9
[i
].d
,&mpt1
,p
); __mul(&mpy
,&mpt1
,&mperr
,p
);
224 __add(&mpy
,&mperr
,&mpy1
,p
); __sub(&mpy
,&mperr
,&mpy2
,p
);
225 __mp_dbl(&mpy1
,&y1
,p
); __mp_dbl(&mpy2
,&y2
,p
);
226 if (y1
==y2
) return y1
;
228 return y1
; /*if unpossible to do exact computing */
231 #ifdef NO_LONG_DOUBLE
232 weak_alias (atan
, atanl
)