summaryrefslogtreecommitdiff
path: root/sysdeps/ieee754/dbl-64/s_atan.c
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps/ieee754/dbl-64/s_atan.c')
-rw-r--r--sysdeps/ieee754/dbl-64/s_atan.c364
1 files changed, 211 insertions, 153 deletions
diff --git a/sysdeps/ieee754/dbl-64/s_atan.c b/sysdeps/ieee754/dbl-64/s_atan.c
index cad3ba12a8..997e29b5e4 100644
--- a/sysdeps/ieee754/dbl-64/s_atan.c
+++ b/sysdeps/ieee754/dbl-64/s_atan.c
@@ -1,163 +1,221 @@
-/* @(#)s_atan.c 5.1 93/09/24 */
/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * IBM Accurate Mathematical Library
+ * Copyright (c) International Business Machines Corp., 2001
*
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-/* Modified by Naohiko Shimizu/Tokai University, Japan 1997/08/25,
- for performance improvement on pipelined processors.
-*/
-
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: s_atan.c,v 1.8 1995/05/10 20:46:45 jtc Exp $";
-#endif
-
-/* atan(x)
- * Method
- * 1. Reduce x to positive by atan(x) = -atan(-x).
- * 2. According to the integer k=4t+0.25 chopped, t=x, the argument
- * is further reduced to one of the following intervals and the
- * arctangent of t is evaluated by the corresponding formula:
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
*
- * [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
- * [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
- * [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
- * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
- * [39/16,INF] atan(x) = atan(INF) + atan( -1/t )
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
*
- * Constants:
- * The hexadecimal values are the intended ones for the following
- * constants. The decimal values may be used, provided that the
- * compiler will convert from decimal to binary accurately enough
- * to produce the hexadecimal values shown.
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
+/************************************************************************/
+/* MODULE_NAME: atnat.c */
+/* */
+/* FUNCTIONS: uatan */
+/* atanMp */
+/* signArctan */
+/* */
+/* */
+/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h atnat.h */
+/* mpatan.c mpatan2.c mpsqrt.c */
+/* uatan.tbl */
+/* */
+/* An ultimate atan() routine. Given an IEEE double machine number x */
+/* it computes the correctly rounded (to nearest) value of atan(x). */
+/* */
+/* Assumption: Machine arithmetic operations are performed in */
+/* round to nearest mode of IEEE 754 standard. */
+/* */
+/************************************************************************/
-#include "math.h"
-#include "math_private.h"
+#include "dla.h"
+#include "mpa.h"
+#include "MathLib.h"
+#include "uatan.tbl"
+#include "atnat.h"
+
+void __mpatan(mp_no *,mp_no *,int); /* see definition in mpatan.c */
+static double atanMp(double,const int[]);
+double __signArctan(double,double);
+/* An ultimate atan() routine. Given an IEEE double machine number x, */
+/* routine computes the correctly rounded (to nearest) value of atan(x). */
+double atan(double x) {
+
+
+ double cor,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,u,u2,u3,
+ v,vv,w,ww,y,yy,y1,y2,z,zz;
+ int i,ux,dx,p;
+ static const int pr[M]={6,8,10,32};
+ number num;
+
+ mp_no mpt1,mpx,mpy,mpy1,mpy2,mperr;
+
+ num.d = x; ux = num.i[HIGH_HALF]; dx = num.i[LOW_HALF];
+
+ /* x=NaN */
+ if (((ux&0x7ff00000)==0x7ff00000) && (((ux&0x000fffff)|dx)!=0x00000000))
+ return x+x;
+
+ /* Regular values of x, including denormals +-0 and +-INF */
+ u = (x<ZERO) ? -x : x;
+ if (u<C) {
+ if (u<B) {
+ if (u<A) { /* u < A */
+ return x; }
+ else { /* A <= u < B */
+ v=x*x; yy=x*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
+ if ((y=x+(yy-U1*x)) == x+(yy+U1*x)) return y;
+
+ EMULV(x,x,v,vv,t1,t2,t3,t4,t5) /* v+vv=x^2 */
+ s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
+ ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
+ MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
+ MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
+ MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
+ MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ MUL2(x,ZERO,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(x,ZERO,s2,ss2,s1,ss1,t1,t2)
+ if ((y=s1+(ss1-U5*s1)) == s1+(ss1+U5*s1)) return y;
+
+ return atanMp(x,pr);
+ } }
+ else { /* B <= u < C */
+ i=(TWO52+TWO8*u)-TWO52; i-=16;
+ z=u-cij[i][0].d;
+ yy=z*(cij[i][2].d+z*(cij[i][3].d+z*(cij[i][4].d+
+ z*(cij[i][5].d+z* cij[i][6].d))));
+ t1=cij[i][1].d;
+ if (i<112) {
+ if (i<48) u2=U21; /* u < 1/4 */
+ else u2=U22; } /* 1/4 <= u < 1/2 */
+ else {
+ if (i<176) u2=U23; /* 1/2 <= u < 3/4 */
+ else u2=U24; } /* 3/4 <= u <= 1 */
+ if ((y=t1+(yy-u2*t1)) == t1+(yy+u2*t1)) return __signArctan(x,y);
+
+ z=u-hij[i][0].d;
+ s1=z*(hij[i][11].d+z*(hij[i][12].d+z*(hij[i][13].d+
+ z*(hij[i][14].d+z* hij[i][15].d))));
+ ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
+ MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
+ MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
+ MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
+ MUL2(z,ZERO,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
+ if ((y=s2+(ss2-U6*s2)) == s2+(ss2+U6*s2)) return __signArctan(x,y);
+
+ return atanMp(x,pr);
+ }
+ }
+ else {
+ if (u<D) { /* C <= u < D */
+ w=ONE/u;
+ EMULV(w,u,t1,t2,t3,t4,t5,t6,t7)
+ ww=w*((ONE-t1)-t2);
+ i=(TWO52+TWO8*w)-TWO52; i-=16;
+ z=(w-cij[i][0].d)+ww;
+ yy=HPI1-z*(cij[i][2].d+z*(cij[i][3].d+z*(cij[i][4].d+
+ z*(cij[i][5].d+z* cij[i][6].d))));
+ t1=HPI-cij[i][1].d;
+ if (i<112) u3=U31; /* w < 1/2 */
+ else u3=U32; /* w >= 1/2 */
+ if ((y=t1+(yy-u3)) == t1+(yy+u3)) return __signArctan(x,y);
+
+ DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+ t1=w-hij[i][0].d;
+ EADD(t1,ww,z,zz)
+ s1=z*(hij[i][11].d+z*(hij[i][12].d+z*(hij[i][13].d+
+ z*(hij[i][14].d+z* hij[i][15].d))));
+ ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
+ MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
+ MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
+ MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
+ MUL2(z,zz,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
+ SUB2(HPI,HPI1,s2,ss2,s1,ss1,t1,t2)
+ if ((y=s1+(ss1-U7)) == s1+(ss1+U7)) return __signArctan(x,y);
+
+ return atanMp(x,pr);
+ }
+ else {
+ if (u<E) { /* D <= u < E */
+ w=ONE/u; v=w*w;
+ EMULV(w,u,t1,t2,t3,t4,t5,t6,t7)
+ yy=w*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
+ ww=w*((ONE-t1)-t2);
+ ESUB(HPI,w,t3,cor)
+ yy=((HPI1+cor)-ww)-yy;
+ if ((y=t3+(yy-U4)) == t3+(yy+U4)) return __signArctan(x,y);
+
+ DIV2(ONE,ZERO,u,ZERO,w,ww,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+ MUL2(w,ww,w,ww,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
+ s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
+ ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
+ MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
+ MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
+ MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
+ MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
+ MUL2(w,ww,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(w,ww,s2,ss2,s1,ss1,t1,t2)
+ SUB2(HPI,HPI1,s1,ss1,s2,ss2,t1,t2)
+ if ((y=s2+(ss2-U8)) == s2+(ss2+U8)) return __signArctan(x,y);
+
+ return atanMp(x,pr);
+ }
+ else {
+ /* u >= E */
+ if (x>0) return HPI;
+ else return MHPI; }
+ }
+ }
-#ifdef __STDC__
-static const double atanhi[] = {
-#else
-static double atanhi[] = {
-#endif
- 4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */
- 7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */
- 9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */
- 1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */
-};
-
-#ifdef __STDC__
-static const double atanlo[] = {
-#else
-static double atanlo[] = {
-#endif
- 2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */
- 3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */
- 1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */
- 6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */
-};
-
-#ifdef __STDC__
-static const double aT[] = {
-#else
-static double aT[] = {
-#endif
- 3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */
- -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */
- 1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */
- -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */
- 9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */
- -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */
- 6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */
- -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */
- 4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */
- -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */
- 1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
-};
-
-#ifdef __STDC__
- static const double
-#else
- static double
-#endif
-one = 1.0,
-huge = 1.0e300;
-
-#ifdef __STDC__
- double __atan(double x)
-#else
- double __atan(x)
- double x;
-#endif
-{
- double w,s1,z,s,w2,w4,s11,s12,s13,s21,s22,s23;
- int32_t ix,hx,id;
-
- GET_HIGH_WORD(hx,x);
- ix = hx&0x7fffffff;
- if(ix>=0x44100000) { /* if |x| >= 2^66 */
- u_int32_t low;
- GET_LOW_WORD(low,x);
- if(ix>0x7ff00000||
- (ix==0x7ff00000&&(low!=0)))
- return x+x; /* NaN */
- if(hx>0) return atanhi[3]+atanlo[3];
- else return -atanhi[3]-atanlo[3];
- } if (ix < 0x3fdc0000) { /* |x| < 0.4375 */
- if (ix < 0x3e200000) { /* |x| < 2^-29 */
- if(huge+x>one) return x; /* raise inexact */
- }
- id = -1;
- } else {
- x = fabs(x);
- if (ix < 0x3ff30000) { /* |x| < 1.1875 */
- if (ix < 0x3fe60000) { /* 7/16 <=|x|<11/16 */
- id = 0; x = (2.0*x-one)/(2.0+x);
- } else { /* 11/16<=|x|< 19/16 */
- id = 1; x = (x-one)/(x+one);
- }
- } else {
- if (ix < 0x40038000) { /* |x| < 2.4375 */
- id = 2; x = (x-1.5)/(one+1.5*x);
- } else { /* 2.4375 <= |x| < 2^66 */
- id = 3; x = -1.0/x;
- }
- }}
- /* end of argument reduction */
- z = x*x;
- w = z*z;
- /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
-#ifdef DO_NOT_USE_THIS
- s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10])))));
- s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9]))));
- if (id<0) return x - x*(s1+s2);
- else {
- z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
- return (hx<0)? -z:z;
- }
-#else
- s11 = aT[8]+w*aT[10]; w2=w*w;
- s12 = aT[4]+w*aT[6]; w4=w2*w2;
- s13 = aT[0]+w*aT[2];
- s21 = aT[7]+w*aT[9];
- s22 = aT[3]+w*aT[5];
- s23 = w*aT[1];
- s1 = s13 + w2*s12 + w4*s11;
- s = s23 + w2*s22 + w4*s21 + z*s1;
- if (id<0) return x - x*(s);
- else {
- z = atanhi[id] - ((x*(s) - atanlo[id]) - x);
- return (hx<0)? -z:z;
- }
-#endif
}
-weak_alias (__atan, atan)
+
+
+ /* Fix the sign of y and return */
+double __signArctan(double x,double y){
+
+ if (x<ZERO) return -y;
+ else return y;
+}
+
+ /* Final stages. Compute atan(x) by multiple precision arithmetic */
+static double atanMp(double x,const int pr[]){
+ mp_no mpx,mpy,mpy2,mperr,mpt1,mpy1;
+ double y1,y2;
+ int i,p;
+
+for (i=0; i<M; i++) {
+ p = pr[i];
+ dbl_mp(x,&mpx,p); __mpatan(&mpx,&mpy,p);
+ dbl_mp(u9[i].d,&mpt1,p); mul(&mpy,&mpt1,&mperr,p);
+ add(&mpy,&mperr,&mpy1,p); sub(&mpy,&mperr,&mpy2,p);
+ mp_dbl(&mpy1,&y1,p); mp_dbl(&mpy2,&y2,p);
+ if (y1==y2) return y1;
+ }
+ return y1; /*if unpossible to do exact computing */
+}
+
#ifdef NO_LONG_DOUBLE
-strong_alias (__atan, __atanl)
-weak_alias (__atan, atanl)
+weak_alias (atan, atanl)
#endif
View Lorry Controller Status