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 /************************************************************************/
46 void __mpatan(mp_no *,mp_no *,int); /* see definition in mpatan.c */
47 static double atanMp(double,const int[]);
48 double __signArctan(double,double);
49 /* An ultimate atan() routine. Given an IEEE double machine number x, */
50 /* routine computes the correctly rounded (to nearest) value of atan(x). */
51 double atan(double x) {
54 double cor,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,u,u2,u3,
63 static const int pr[M]={6,8,10,32};
66 mp_no mpt1,mpx,mpy,mpy1,mpy2,mperr;
69 num.d = x; ux = num.i[HIGH_HALF]; dx = num.i[LOW_HALF];
72 if (((ux&0x7ff00000)==0x7ff00000) && (((ux&0x000fffff)|dx)!=0x00000000))
75 /* Regular values of x, including denormals +-0 and +-INF */
76 u = (x<ZERO) ? -x : x;
79 if (u<A) { /* u < A */
81 else { /* A <= u < B */
82 v=x*x; yy=x*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
83 if ((y=x+(yy-U1*x)) == x+(yy+U1*x)) return y;
85 EMULV(x,x,v,vv,t1,t2,t3,t4,t5) /* v+vv=x^2 */
86 s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
87 ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
88 MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
89 ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
90 MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
91 ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
92 MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
93 ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
94 MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
95 MUL2(x,ZERO,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
96 ADD2(x,ZERO,s2,ss2,s1,ss1,t1,t2)
97 if ((y=s1+(ss1-U5*s1)) == s1+(ss1+U5*s1)) return y;
101 else { /* B <= u < C */
102 i=(TWO52+TWO8*u)-TWO52; i-=16;
104 yy=z*(cij[i][2].d+z*(cij[i][3].d+z*(cij[i][4].d+
105 z*(cij[i][5].d+z* cij[i][6].d))));
108 if (i<48) u2=U21; /* u < 1/4 */
109 else u2=U22; } /* 1/4 <= u < 1/2 */
111 if (i<176) u2=U23; /* 1/2 <= u < 3/4 */
112 else u2=U24; } /* 3/4 <= u <= 1 */
113 if ((y=t1+(yy-u2*t1)) == t1+(yy+u2*t1)) return __signArctan(x,y);
116 s1=z*(hij[i][11].d+z*(hij[i][12].d+z*(hij[i][13].d+
117 z*(hij[i][14].d+z* hij[i][15].d))));
118 ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
119 MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
120 ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
121 MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
122 ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
123 MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
124 ADD2(hij[i][3].d,hij[i][4].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][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
127 if ((y=s2+(ss2-U6*s2)) == s2+(ss2+U6*s2)) return __signArctan(x,y);
133 if (u<D) { /* C <= u < D */
135 EMULV(w,u,t1,t2,t3,t4,t5,t6,t7)
137 i=(TWO52+TWO8*w)-TWO52; i-=16;
138 z=(w-cij[i][0].d)+ww;
139 yy=HPI1-z*(cij[i][2].d+z*(cij[i][3].d+z*(cij[i][4].d+
140 z*(cij[i][5].d+z* cij[i][6].d))));
142 if (i<112) u3=U31; /* w < 1/2 */
143 else u3=U32; /* w >= 1/2 */
144 if ((y=t1+(yy-u3)) == t1+(yy+u3)) return __signArctan(x,y);
146 DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
149 s1=z*(hij[i][11].d+z*(hij[i][12].d+z*(hij[i][13].d+
150 z*(hij[i][14].d+z* hij[i][15].d))));
151 ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
152 MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
153 ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
154 MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
155 ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
156 MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
157 ADD2(hij[i][3].d,hij[i][4].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][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
160 SUB2(HPI,HPI1,s2,ss2,s1,ss1,t1,t2)
161 if ((y=s1+(ss1-U7)) == s1+(ss1+U7)) return __signArctan(x,y);
166 if (u<E) { /* D <= u < E */
168 EMULV(w,u,t1,t2,t3,t4,t5,t6,t7)
169 yy=w*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
172 yy=((HPI1+cor)-ww)-yy;
173 if ((y=t3+(yy-U4)) == t3+(yy+U4)) return __signArctan(x,y);
175 DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
176 MUL2(w,ww,w,ww,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
177 s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
178 ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
179 MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
180 ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
181 MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
182 ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
183 MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
184 ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
185 MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
186 MUL2(w,ww,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
187 ADD2(w,ww,s2,ss2,s1,ss1,t1,t2)
188 SUB2(HPI,HPI1,s1,ss1,s2,ss2,t1,t2)
189 if ((y=s2+(ss2-U8)) == s2+(ss2+U8)) return __signArctan(x,y);
203 /* Fix the sign of y and return */
204 double __signArctan(double x,double y){
206 if (x<ZERO) return -y;
210 /* Final stages. Compute atan(x) by multiple precision arithmetic */
211 static double atanMp(double x,const int pr[]){
212 mp_no mpx,mpy,mpy2,mperr,mpt1,mpy1;
216 for (i=0; i<M; i++) {
218 __dbl_mp(x,&mpx,p); __mpatan(&mpx,&mpy,p);
219 __dbl_mp(u9[i].d,&mpt1,p); __mul(&mpy,&mpt1,&mperr,p);
220 __add(&mpy,&mperr,&mpy1,p); __sub(&mpy,&mperr,&mpy2,p);
221 __mp_dbl(&mpy1,&y1,p); __mp_dbl(&mpy2,&y2,p);
222 if (y1==y2) return y1;
224 return y1; /*if unpossible to do exact computing */
227 #ifdef NO_LONG_DOUBLE
228 weak_alias (atan, atanl)