Back to index

glibc  2.9
s_atan.c
Go to the documentation of this file.
00001 /*
00002  * IBM Accurate Mathematical Library
00003  * written by International Business Machines Corp.
00004  * Copyright (C) 2001 Free Software Foundation
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU Lesser General Public License as published by
00008  * the Free Software Foundation; either version 2.1 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU Lesser General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU Lesser General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00019  */
00020 /************************************************************************/
00021 /*  MODULE_NAME: atnat.c                                                */
00022 /*                                                                      */
00023 /*  FUNCTIONS:  uatan                                                   */
00024 /*              atanMp                                                  */
00025 /*              signArctan                                              */
00026 /*                                                                      */
00027 /*                                                                      */
00028 /*  FILES NEEDED: dla.h endian.h mpa.h mydefs.h atnat.h                 */
00029 /*                mpatan.c mpatan2.c mpsqrt.c                           */
00030 /*                uatan.tbl                                             */
00031 /*                                                                      */
00032 /* An ultimate atan() routine. Given an IEEE double machine number x    */
00033 /* it computes the correctly rounded (to nearest) value of atan(x).     */
00034 /*                                                                      */
00035 /* Assumption: Machine arithmetic operations are performed in           */
00036 /* round to nearest mode of IEEE 754 standard.                          */
00037 /*                                                                      */
00038 /************************************************************************/
00039 
00040 #include "dla.h"
00041 #include "mpa.h"
00042 #include "MathLib.h"
00043 #include "uatan.tbl"
00044 #include "atnat.h"
00045 #include "math.h"
00046 
00047 void __mpatan(mp_no *,mp_no *,int);          /* see definition in mpatan.c */
00048 static double atanMp(double,const int[]);
00049 double __signArctan(double,double);
00050 /* An ultimate atan() routine. Given an IEEE double machine number x,    */
00051 /* routine computes the correctly rounded (to nearest) value of atan(x). */
00052 double atan(double x) {
00053 
00054 
00055   double cor,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,u,u2,u3,
00056          v,vv,w,ww,y,yy,z,zz;
00057 #if 0
00058   double y1,y2;
00059 #endif
00060   int i,ux,dx;
00061 #if 0
00062   int p;
00063 #endif
00064   static const int pr[M]={6,8,10,32};
00065   number num;
00066 #if 0
00067   mp_no mpt1,mpx,mpy,mpy1,mpy2,mperr;
00068 #endif
00069 
00070   num.d = x;  ux = num.i[HIGH_HALF];  dx = num.i[LOW_HALF];
00071 
00072   /* x=NaN */
00073   if (((ux&0x7ff00000)==0x7ff00000) && (((ux&0x000fffff)|dx)!=0x00000000))
00074     return x+x;
00075 
00076   /* Regular values of x, including denormals +-0 and +-INF */
00077   u = (x<ZERO) ? -x : x;
00078   if (u<C) {
00079     if (u<B) {
00080       if (u<A) {                                           /* u < A */
00081          return x; }
00082       else {                                               /* A <= u < B */
00083         v=x*x;  yy=x*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
00084         if ((y=x+(yy-U1*x)) == x+(yy+U1*x))  return y;
00085 
00086         EMULV(x,x,v,vv,t1,t2,t3,t4,t5)                       /* v+vv=x^2 */
00087         s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
00088         ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
00089         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00090         ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
00091         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00092         ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
00093         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00094         ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
00095         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00096         MUL2(x,ZERO,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
00097         ADD2(x,ZERO,s2,ss2,s1,ss1,t1,t2)
00098         if ((y=s1+(ss1-U5*s1)) == s1+(ss1+U5*s1))  return y;
00099 
00100         return atanMp(x,pr);
00101       } }
00102     else {  /* B <= u < C */
00103       i=(TWO52+TWO8*u)-TWO52;  i-=16;
00104       z=u-cij[i][0].d;
00105       yy=z*(cij[i][2].d+z*(cij[i][3].d+z*(cij[i][4].d+
00106                         z*(cij[i][5].d+z* cij[i][6].d))));
00107       t1=cij[i][1].d;
00108       if (i<112) {
00109         if (i<48)  u2=U21;    /* u < 1/4        */
00110         else       u2=U22; }  /* 1/4 <= u < 1/2 */
00111       else {
00112         if (i<176) u2=U23;    /* 1/2 <= u < 3/4 */
00113         else       u2=U24; }  /* 3/4 <= u <= 1  */
00114       if ((y=t1+(yy-u2*t1)) == t1+(yy+u2*t1))  return __signArctan(x,y);
00115 
00116       z=u-hij[i][0].d;
00117       s1=z*(hij[i][11].d+z*(hij[i][12].d+z*(hij[i][13].d+
00118          z*(hij[i][14].d+z* hij[i][15].d))));
00119       ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
00120       MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00121       ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
00122       MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00123       ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
00124       MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00125       ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
00126       MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00127       ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
00128       if ((y=s2+(ss2-U6*s2)) == s2+(ss2+U6*s2))  return __signArctan(x,y);
00129 
00130       return atanMp(x,pr);
00131     }
00132   }
00133   else {
00134     if (u<D) { /* C <= u < D */
00135       w=ONE/u;
00136       EMULV(w,u,t1,t2,t3,t4,t5,t6,t7)
00137       ww=w*((ONE-t1)-t2);
00138       i=(TWO52+TWO8*w)-TWO52;  i-=16;
00139       z=(w-cij[i][0].d)+ww;
00140       yy=HPI1-z*(cij[i][2].d+z*(cij[i][3].d+z*(cij[i][4].d+
00141                              z*(cij[i][5].d+z* cij[i][6].d))));
00142       t1=HPI-cij[i][1].d;
00143       if (i<112)  u3=U31;  /* w <  1/2 */
00144       else        u3=U32;  /* w >= 1/2 */
00145       if ((y=t1+(yy-u3)) == t1+(yy+u3))  return __signArctan(x,y);
00146 
00147       DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
00148       t1=w-hij[i][0].d;
00149       EADD(t1,ww,z,zz)
00150       s1=z*(hij[i][11].d+z*(hij[i][12].d+z*(hij[i][13].d+
00151          z*(hij[i][14].d+z* hij[i][15].d))));
00152       ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
00153       MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00154       ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
00155       MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00156       ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
00157       MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00158       ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
00159       MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00160       ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
00161       SUB2(HPI,HPI1,s2,ss2,s1,ss1,t1,t2)
00162       if ((y=s1+(ss1-U7)) == s1+(ss1+U7))  return __signArctan(x,y);
00163 
00164     return atanMp(x,pr);
00165     }
00166     else {
00167       if (u<E) { /* D <= u < E */
00168         w=ONE/u;   v=w*w;
00169         EMULV(w,u,t1,t2,t3,t4,t5,t6,t7)
00170         yy=w*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
00171         ww=w*((ONE-t1)-t2);
00172         ESUB(HPI,w,t3,cor)
00173         yy=((HPI1+cor)-ww)-yy;
00174         if ((y=t3+(yy-U4)) == t3+(yy+U4))  return __signArctan(x,y);
00175 
00176         DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
00177         MUL2(w,ww,w,ww,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
00178         s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
00179         ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
00180         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00181         ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
00182         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00183         ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
00184         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00185         ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
00186         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00187         MUL2(w,ww,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
00188         ADD2(w,ww,s2,ss2,s1,ss1,t1,t2)
00189         SUB2(HPI,HPI1,s1,ss1,s2,ss2,t1,t2)
00190         if ((y=s2+(ss2-U8)) == s2+(ss2+U8))  return __signArctan(x,y);
00191 
00192       return atanMp(x,pr);
00193       }
00194       else {
00195         /* u >= E */
00196         if (x>0) return  HPI;
00197         else     return MHPI; }
00198     }
00199   }
00200 
00201 }
00202 
00203 
00204   /* Fix the sign of y and return */
00205 double  __signArctan(double x,double y){
00206 
00207     if (x<ZERO) return -y;
00208     else        return  y;
00209 }
00210 
00211  /* Final stages. Compute atan(x) by multiple precision arithmetic */
00212 static double atanMp(double x,const int pr[]){
00213   mp_no mpx,mpy,mpy2,mperr,mpt1,mpy1;
00214   double y1,y2;
00215   int i,p;
00216 
00217 for (i=0; i<M; i++) {
00218     p = pr[i];
00219     __dbl_mp(x,&mpx,p);          __mpatan(&mpx,&mpy,p);
00220     __dbl_mp(u9[i].d,&mpt1,p);   __mul(&mpy,&mpt1,&mperr,p);
00221     __add(&mpy,&mperr,&mpy1,p);  __sub(&mpy,&mperr,&mpy2,p);
00222     __mp_dbl(&mpy1,&y1,p);       __mp_dbl(&mpy2,&y2,p);
00223     if (y1==y2)   return y1;
00224   }
00225   return y1; /*if unpossible to do exact computing */
00226 }
00227 
00228 #ifdef NO_LONG_DOUBLE
00229 weak_alias (atan, atanl)
00230 #endif