sysdeps/ieee754/dbl-64/s_atan.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001, 2011 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, see <http://www.gnu.org/licenses/>.
  18  */
  19 /************************************************************************/
  20 /*  MODULE_NAME: atnat.c                                                */
  21 /*                                                                      */
  22 /*  FUNCTIONS:  uatan                                                   */
  23 /*              atanMp                                                  */
  24 /*              signArctan                                              */
  25 /*                                                                      */
  26 /*                                                                      */
  27 /*  FILES NEEDED: dla.h endian.h mpa.h mydefs.h atnat.h                 */
  28 /*                mpatan.c mpatan2.c mpsqrt.c                           */
  29 /*                uatan.tbl                                             */
  30 /*                                                                      */
  31 /* An ultimate atan() routine. Given an IEEE double machine number x    */
  32 /* it computes the correctly rounded (to nearest) value of atan(x).     */
  33 /*                                                                      */
  34 /* Assumption: Machine arithmetic operations are performed in           */
  35 /* round to nearest mode of IEEE 754 standard.                          */
  36 /*                                                                      */
  37 /************************************************************************/
  38
  39 #include <dla.h>
  40 #include "mpa.h"
  41 #include "MathLib.h"
  42 #include "uatan.tbl"
  43 #include "atnat.h"
  44 #include "math.h"
  45
  46 void __mpatan(mp_no *,mp_no *,int);          /* see definition in mpatan.c */
  47 static double atanMp(double,const int[]);
  48
  49   /* Fix the sign of y and return */
  50 static double  __signArctan(double x,double y){
  51   return __copysign(y, x);
  52 }
  53
  54
  55 /* An ultimate atan() routine. Given an IEEE double machine number x,    */
  56 /* routine computes the correctly rounded (to nearest) value of atan(x). */
  57 double atan(double x) {
  58
  59
  60   double cor,s1,ss1,s2,ss2,t1,t2,t3,t7,t8,t9,t10,u,u2,u3,
  61          v,vv,w,ww,y,yy,z,zz;
  62 #ifndef DLA_FMS
  63   double t4,t5,t6;
  64 #endif
  65 #if 0
  66   double y1,y2;
  67 #endif
  68   int i,ux,dx;
  69 #if 0
  70   int p;
  71 #endif
  72   static const int pr[M]={6,8,10,32};
  73   number num;
  74 #if 0
  75   mp_no mpt1,mpx,mpy,mpy1,mpy2,mperr;
  76 #endif
  77
  78   num.d = x;  ux = num.i[HIGH_HALF];  dx = num.i[LOW_HALF];
  79
  80   /* x=NaN */
  81   if (((ux&0x7ff00000)==0x7ff00000) && (((ux&0x000fffff)|dx)!=0x00000000))
  82     return x+x;
  83
  84   /* Regular values of x, including denormals +-0 and +-INF */
  85   u = (x<ZERO) ? -x : x;
  86   if (u<C) {
  87     if (u<B) {
  88       if (u<A) {                                           /* u < A */
  89          return x; }
  90       else {                                               /* A <= u < B */
  91         v=x*x;  yy=x*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
  92         if ((y=x+(yy-U1*x)) == x+(yy+U1*x))  return y;
  93
  94         EMULV(x,x,v,vv,t1,t2,t3,t4,t5)                       /* v+vv=x^2 */
  95         s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
  96         ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
  97         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
  98         ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
  99         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 100         ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
 101         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 102         ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
 103         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 104         MUL2(x,ZERO,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
 105         ADD2(x,ZERO,s2,ss2,s1,ss1,t1,t2)
 106         if ((y=s1+(ss1-U5*s1)) == s1+(ss1+U5*s1))  return y;
 107
 108         return atanMp(x,pr);
 109       } }
 110     else {  /* B <= u < C */
 111       i=(TWO52+TWO8*u)-TWO52;  i-=16;
 112       z=u-cij[i][0].d;
 113       yy=z*(cij[i][2].d+z*(cij[i][3].d+z*(cij[i][4].d+
 114                         z*(cij[i][5].d+z* cij[i][6].d))));
 115       t1=cij[i][1].d;
 116       if (i<112) {
 117         if (i<48)  u2=U21;    /* u < 1/4        */
 118         else       u2=U22; }  /* 1/4 <= u < 1/2 */
 119       else {
 120         if (i<176) u2=U23;    /* 1/2 <= u < 3/4 */
 121         else       u2=U24; }  /* 3/4 <= u <= 1  */
 122       if ((y=t1+(yy-u2*t1)) == t1+(yy+u2*t1))  return __signArctan(x,y);
 123
 124       z=u-hij[i][0].d;
 125       s1=z*(hij[i][11].d+z*(hij[i][12].d+z*(hij[i][13].d+
 126          z*(hij[i][14].d+z* hij[i][15].d))));
 127       ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
 128       MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 129       ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
 130       MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 131       ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
 132       MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 133       ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
 134       MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 135       ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
 136       if ((y=s2+(ss2-U6*s2)) == s2+(ss2+U6*s2))  return __signArctan(x,y);
 137
 138       return atanMp(x,pr);
 139     }
 140   }
 141   else {
 142     if (u<D) { /* C <= u < D */
 143       w=ONE/u;
 144       EMULV(w,u,t1,t2,t3,t4,t5,t6,t7)
 145       ww=w*((ONE-t1)-t2);
 146       i=(TWO52+TWO8*w)-TWO52;  i-=16;
 147       z=(w-cij[i][0].d)+ww;
 148       yy=HPI1-z*(cij[i][2].d+z*(cij[i][3].d+z*(cij[i][4].d+
 149                              z*(cij[i][5].d+z* cij[i][6].d))));
 150       t1=HPI-cij[i][1].d;
 151       if (i<112)  u3=U31;  /* w <  1/2 */
 152       else        u3=U32;  /* w >= 1/2 */
 153       if ((y=t1+(yy-u3)) == t1+(yy+u3))  return __signArctan(x,y);
 154
 155       DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
 156       t1=w-hij[i][0].d;
 157       EADD(t1,ww,z,zz)
 158       s1=z*(hij[i][11].d+z*(hij[i][12].d+z*(hij[i][13].d+
 159          z*(hij[i][14].d+z* hij[i][15].d))));
 160       ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
 161       MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 162       ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
 163       MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 164       ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
 165       MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 166       ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
 167       MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 168       ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
 169       SUB2(HPI,HPI1,s2,ss2,s1,ss1,t1,t2)
 170       if ((y=s1+(ss1-U7)) == s1+(ss1+U7))  return __signArctan(x,y);
 171
 172     return atanMp(x,pr);
 173     }
 174     else {
 175       if (u<E) { /* D <= u < E */
 176         w=ONE/u;   v=w*w;
 177         EMULV(w,u,t1,t2,t3,t4,t5,t6,t7)
 178         yy=w*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
 179         ww=w*((ONE-t1)-t2);
 180         ESUB(HPI,w,t3,cor)
 181         yy=((HPI1+cor)-ww)-yy;
 182         if ((y=t3+(yy-U4)) == t3+(yy+U4))  return __signArctan(x,y);
 183
 184         DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
 185         MUL2(w,ww,w,ww,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
 186         s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
 187         ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
 188         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 189         ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
 190         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 191         ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
 192         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 193         ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
 194         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
 195         MUL2(w,ww,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
 196         ADD2(w,ww,s2,ss2,s1,ss1,t1,t2)
 197         SUB2(HPI,HPI1,s1,ss1,s2,ss2,t1,t2)
 198         if ((y=s2+(ss2-U8)) == s2+(ss2+U8))  return __signArctan(x,y);
 199
 200       return atanMp(x,pr);
 201       }
 202       else {
 203         /* u >= E */
 204         if (x>0) return  HPI;
 205         else     return MHPI; }
 206     }
 207   }
 208
 209 }
 210
 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;
 214   double y1,y2;
 215   int i,p;
 216
 217 for (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
 229 weak_alias (atan, atanl)
 230 #endif