1/* 2 * IBM Accurate Mathematical Library 3 * written by International Business Machines Corp. 4 * Copyright (C) 2001 Free Software Foundation 5 * 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. 10 * 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. 15 * 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. 19 */ 20/************************************************************************/ 21/* MODULE_NAME: atnat.c */ 22/* */ 23/* FUNCTIONS: uatan */ 24/* atanMp */ 25/* signArctan */ 26/* */ 27/* */ 28/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h atnat.h */ 29/* mpatan.c mpatan2.c mpsqrt.c */ 30/* uatan.tbl */ 31/* */ 32/* An ultimate atan() routine. Given an IEEE double machine number x */ 33/* it computes the correctly rounded (to nearest) value of atan(x). */ 34/* */ 35/* Assumption: Machine arithmetic operations are performed in */ 36/* round to nearest mode of IEEE 754 standard. */ 37/* */ 38/************************************************************************/ 39 40#include "dla.h" 41#include "mpa.h" 42#include "MathLib.h" 43#include "uatan.tbl" 44#include "atnat.h" 45#include "math.h" 46 47void __mpatan(mp_no *,mp_no *,int); /* see definition in mpatan.c */ 48static double atanMp(double,const int[]); 49double __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). */ 52double atan(double x) { 53 54 55 double cor,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,u,u2,u3, 56 v,vv,w,ww,y,yy,z,zz; 57#if 0 58 double y1,y2; 59#endif 60 int i,ux,dx; 61#if 0 62 int p; 63#endif 64 static const int pr[M]={6,8,10,32}; 65 number num; 66#if 0 67 mp_no mpt1,mpx,mpy,mpy1,mpy2,mperr; 68#endif 69 70 num.d = x; ux = num.i[HIGH_HALF]; dx = num.i[LOW_HALF]; 71 72 /* x=NaN */ 73 if (((ux&0x7ff00000)==0x7ff00000) && (((ux&0x000fffff)|dx)!=0x00000000)) 74 return x+x; 75 76 /* Regular values of x, including denormals +-0 and +-INF */ 77 u = (x<ZERO) ? -x : x; 78 if (u<C) { 79 if (u<B) { 80 if (u<A) { /* u < A */ 81 return x; } 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; 85 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; 99 100 return atanMp(x,pr); 101 } } 102 else { /* B <= u < C */ 103 i=(TWO52+TWO8*u)-TWO52; i-=16; 104 z=u-cij[i][0].d; 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)))); 107 t1=cij[i][1].d; 108 if (i<112) { 109 if (i<48) u2=U21; /* u < 1/4 */ 110 else u2=U22; } /* 1/4 <= u < 1/2 */ 111 else { 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); 115 116 z=u-hij[i][0].d; 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); 129 130 return atanMp(x,pr); 131 } 132 } 133 else { 134 if (u<D) { /* C <= u < D */ 135 w=ONE/u; 136 EMULV(w,u,t1,t2,t3,t4,t5,t6,t7) 137 ww=w*((ONE-t1)-t2); 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)))); 142 t1=HPI-cij[i][1].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); 146 147 DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10) 148 t1=w-hij[i][0].d; 149 EADD(t1,ww,z,zz) 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); 163 164 return atanMp(x,pr); 165 } 166 else { 167 if (u<E) { /* D <= u < E */ 168 w=ONE/u; v=w*w; 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))))); 171 ww=w*((ONE-t1)-t2); 172 ESUB(HPI,w,t3,cor) 173 yy=((HPI1+cor)-ww)-yy; 174 if ((y=t3+(yy-U4)) == t3+(yy+U4)) return __signArctan(x,y); 175 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); 191 192 return atanMp(x,pr); 193 } 194 else { 195 /* u >= E */ 196 if (x>0) return HPI; 197 else return MHPI; } 198 } 199 } 200 201} 202 203 204 /* Fix the sign of y and return */ 205double __signArctan(double x,double y){ 206 207 if (x<ZERO) return -y; 208 else return y; 209} 210 211 /* Final stages. Compute atan(x) by multiple precision arithmetic */ 212static double atanMp(double x,const int pr[]){ 213 mp_no mpx,mpy,mpy2,mperr,mpt1,mpy1; 214 double y1,y2; 215 int i,p; 216 217for (i=0; i<M; i++) { 218 p = pr[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; 224 } 225 return y1; /*if unpossible to do exact computing */ 226} 227 228#ifdef NO_LONG_DOUBLE 229weak_alias (atan, atanl) 230#endif 231