2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001 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
,t4
,t5
,t6
,t7
,t8
,t9
,t10
,u
,u2
,u3
,
64 static const int pr
[M
]={6,8,10,32};
67 mp_no mpt1
,mpx
,mpy
,mpy1
,mpy2
,mperr
;
70 num
.d
= x
; ux
= num
.i
[HIGH_HALF
]; dx
= num
.i
[LOW_HALF
];
73 if (((ux
&0x7ff00000)==0x7ff00000) && (((ux
&0x000fffff)|dx
)!=0x00000000))
76 /* Regular values of x, including denormals +-0 and +-INF */
77 u
= (x
<ZERO
) ? -x
: x
;
80 if (u
<A
) { /* u < A */
82 else { /* A <= u < B */
83 v
=x
*x
; yy
=x
*v
*(d3
.d
+v
*(d5
.d
+v
*(d7
.d
+v
*(d9
.d
+v
*(d11
.d
+v
*d13
.d
)))));
84 if ((y
=x
+(yy
-U1
*x
)) == x
+(yy
+U1
*x
)) return y
;
86 EMULV(x
,x
,v
,vv
,t1
,t2
,t3
,t4
,t5
) /* v+vv=x^2 */
87 s1
=v
*(f11
.d
+v
*(f13
.d
+v
*(f15
.d
+v
*(f17
.d
+v
*f19
.d
))));
88 ADD2(f9
.d
,ff9
.d
,s1
,ZERO
,s2
,ss2
,t1
,t2
)
89 MUL2(v
,vv
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
90 ADD2(f7
.d
,ff7
.d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
91 MUL2(v
,vv
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
92 ADD2(f5
.d
,ff5
.d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
93 MUL2(v
,vv
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
94 ADD2(f3
.d
,ff3
.d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
95 MUL2(v
,vv
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
96 MUL2(x
,ZERO
,s1
,ss1
,s2
,ss2
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
97 ADD2(x
,ZERO
,s2
,ss2
,s1
,ss1
,t1
,t2
)
98 if ((y
=s1
+(ss1
-U5
*s1
)) == s1
+(ss1
+U5
*s1
)) return y
;
102 else { /* B <= u < C */
103 i
=(TWO52
+TWO8
*u
)-TWO52
; i
-=16;
105 yy
=z
*(cij
[i
][2].d
+z
*(cij
[i
][3].d
+z
*(cij
[i
][4].d
+
106 z
*(cij
[i
][5].d
+z
* cij
[i
][6].d
))));
109 if (i
<48) u2
=U21
; /* u < 1/4 */
110 else u2
=U22
; } /* 1/4 <= u < 1/2 */
112 if (i
<176) u2
=U23
; /* 1/2 <= u < 3/4 */
113 else u2
=U24
; } /* 3/4 <= u <= 1 */
114 if ((y
=t1
+(yy
-u2
*t1
)) == t1
+(yy
+u2
*t1
)) return __signArctan(x
,y
);
117 s1
=z
*(hij
[i
][11].d
+z
*(hij
[i
][12].d
+z
*(hij
[i
][13].d
+
118 z
*(hij
[i
][14].d
+z
* hij
[i
][15].d
))));
119 ADD2(hij
[i
][9].d
,hij
[i
][10].d
,s1
,ZERO
,s2
,ss2
,t1
,t2
)
120 MUL2(z
,ZERO
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
121 ADD2(hij
[i
][7].d
,hij
[i
][8].d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
122 MUL2(z
,ZERO
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
123 ADD2(hij
[i
][5].d
,hij
[i
][6].d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
124 MUL2(z
,ZERO
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
125 ADD2(hij
[i
][3].d
,hij
[i
][4].d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
126 MUL2(z
,ZERO
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
127 ADD2(hij
[i
][1].d
,hij
[i
][2].d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
128 if ((y
=s2
+(ss2
-U6
*s2
)) == s2
+(ss2
+U6
*s2
)) return __signArctan(x
,y
);
134 if (u
<D
) { /* C <= u < D */
136 EMULV(w
,u
,t1
,t2
,t3
,t4
,t5
,t6
,t7
)
138 i
=(TWO52
+TWO8
*w
)-TWO52
; i
-=16;
139 z
=(w
-cij
[i
][0].d
)+ww
;
140 yy
=HPI1
-z
*(cij
[i
][2].d
+z
*(cij
[i
][3].d
+z
*(cij
[i
][4].d
+
141 z
*(cij
[i
][5].d
+z
* cij
[i
][6].d
))));
143 if (i
<112) u3
=U31
; /* w < 1/2 */
144 else u3
=U32
; /* w >= 1/2 */
145 if ((y
=t1
+(yy
-u3
)) == t1
+(yy
+u3
)) return __signArctan(x
,y
);
147 DIV2(ONE
,ZERO
,u
,ZERO
,w
,ww
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
,t9
,t10
)
150 s1
=z
*(hij
[i
][11].d
+z
*(hij
[i
][12].d
+z
*(hij
[i
][13].d
+
151 z
*(hij
[i
][14].d
+z
* hij
[i
][15].d
))));
152 ADD2(hij
[i
][9].d
,hij
[i
][10].d
,s1
,ZERO
,s2
,ss2
,t1
,t2
)
153 MUL2(z
,zz
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
154 ADD2(hij
[i
][7].d
,hij
[i
][8].d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
155 MUL2(z
,zz
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
156 ADD2(hij
[i
][5].d
,hij
[i
][6].d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
157 MUL2(z
,zz
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
158 ADD2(hij
[i
][3].d
,hij
[i
][4].d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
159 MUL2(z
,zz
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
160 ADD2(hij
[i
][1].d
,hij
[i
][2].d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
161 SUB2(HPI
,HPI1
,s2
,ss2
,s1
,ss1
,t1
,t2
)
162 if ((y
=s1
+(ss1
-U7
)) == s1
+(ss1
+U7
)) return __signArctan(x
,y
);
167 if (u
<E
) { /* D <= u < E */
169 EMULV(w
,u
,t1
,t2
,t3
,t4
,t5
,t6
,t7
)
170 yy
=w
*v
*(d3
.d
+v
*(d5
.d
+v
*(d7
.d
+v
*(d9
.d
+v
*(d11
.d
+v
*d13
.d
)))));
173 yy
=((HPI1
+cor
)-ww
)-yy
;
174 if ((y
=t3
+(yy
-U4
)) == t3
+(yy
+U4
)) return __signArctan(x
,y
);
176 DIV2(ONE
,ZERO
,u
,ZERO
,w
,ww
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
,t9
,t10
)
177 MUL2(w
,ww
,w
,ww
,v
,vv
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
178 s1
=v
*(f11
.d
+v
*(f13
.d
+v
*(f15
.d
+v
*(f17
.d
+v
*f19
.d
))));
179 ADD2(f9
.d
,ff9
.d
,s1
,ZERO
,s2
,ss2
,t1
,t2
)
180 MUL2(v
,vv
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
181 ADD2(f7
.d
,ff7
.d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
182 MUL2(v
,vv
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
183 ADD2(f5
.d
,ff5
.d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
184 MUL2(v
,vv
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
185 ADD2(f3
.d
,ff3
.d
,s1
,ss1
,s2
,ss2
,t1
,t2
)
186 MUL2(v
,vv
,s2
,ss2
,s1
,ss1
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
187 MUL2(w
,ww
,s1
,ss1
,s2
,ss2
,t1
,t2
,t3
,t4
,t5
,t6
,t7
,t8
)
188 ADD2(w
,ww
,s2
,ss2
,s1
,ss1
,t1
,t2
)
189 SUB2(HPI
,HPI1
,s1
,ss1
,s2
,ss2
,t1
,t2
)
190 if ((y
=s2
+(ss2
-U8
)) == s2
+(ss2
+U8
)) return __signArctan(x
,y
);
204 /* Fix the sign of y and return */
205 double __signArctan(double x
,double y
){
207 if (x
<ZERO
) return -y
;
211 /* Final stages. Compute atan(x) by multiple precision arithmetic */
212 static double atanMp(double x
,const int pr
[]){
213 mp_no mpx
,mpy
,mpy2
,mperr
,mpt1
,mpy1
;
217 for (i
=0; i
<M
; i
++) {
219 __dbl_mp(x
,&mpx
,p
); __mpatan(&mpx
,&mpy
,p
);
220 __dbl_mp(u9
[i
].d
,&mpt1
,p
); __mul(&mpy
,&mpt1
,&mperr
,p
);
221 __add(&mpy
,&mperr
,&mpy1
,p
); __sub(&mpy
,&mperr
,&mpy2
,p
);
222 __mp_dbl(&mpy1
,&y1
,p
); __mp_dbl(&mpy2
,&y2
,p
);
223 if (y1
==y2
) return y1
;
225 return y1
; /*if unpossible to do exact computing */
228 #ifdef NO_LONG_DOUBLE
229 weak_alias (atan
, atanl
)