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