aboutsummaryrefslogtreecommitdiff
path: root/sysdeps/ieee754/dbl-64/e_atan2.c
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps/ieee754/dbl-64/e_atan2.c')
-rw-r--r--sysdeps/ieee754/dbl-64/e_atan2.c511
1 files changed, 388 insertions, 123 deletions
diff --git a/sysdeps/ieee754/dbl-64/e_atan2.c b/sysdeps/ieee754/dbl-64/e_atan2.c
index ae7d759a9f..11342d87d3 100644
--- a/sysdeps/ieee754/dbl-64/e_atan2.c
+++ b/sysdeps/ieee754/dbl-64/e_atan2.c
@@ -1,130 +1,395 @@
-/* @(#)e_atan2.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.
- * ====================================================
- */
-
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: e_atan2.c,v 1.8 1995/05/10 20:44:51 jtc Exp $";
-#endif
-
-/* __ieee754_atan2(y,x)
- * Method :
- * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
- * 2. Reduce x to positive by (if x and y are unexceptional):
- * ARG (x+iy) = arctan(y/x) ... if x > 0,
- * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0,
- *
- * Special cases:
+ * 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.
*
- * ATAN2((anything), NaN ) is NaN;
- * ATAN2(NAN , (anything) ) is NaN;
- * ATAN2(+-0, +(anything but NaN)) is +-0 ;
- * ATAN2(+-0, -(anything but NaN)) is +-pi ;
- * ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
- * ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
- * ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
- * ATAN2(+-INF,+INF ) is +-pi/4 ;
- * ATAN2(+-INF,-INF ) is +-3pi/4;
- * ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
+ * 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: atnat2.c */
+/* */
+/* FUNCTIONS: uatan2 */
+/* atan2Mp */
+/* signArctan2 */
+/* normalized */
+/* */
+/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h atnat2.h */
+/* mpatan.c mpatan2.c mpsqrt.c */
+/* uatan.tbl */
+/* */
+/* An ultimate atan2() routine. Given two IEEE double machine numbers y,*/
+/* x it computes the correctly rounded (to nearest) value of atan2(y,x).*/
+/* */
+/* Assumption: Machine arithmetic operations are performed in */
+/* round to nearest mode of IEEE 754 standard. */
+/* */
+/************************************************************************/
+
+#include "dla.h"
+#include "mpa.h"
+#include "MathLib.h"
+#include "uatan.tbl"
+#include "atnat2.h"
+/************************************************************************/
+/* An ultimate atan2 routine. Given two IEEE double machine numbers y,x */
+/* it computes the correctly rounded (to nearest) value of atan2(y,x). */
+/* Assumption: Machine arithmetic operations are performed in */
+/* round to nearest mode of IEEE 754 standard. */
+/************************************************************************/
+static double atan2Mp(double ,double ,const int[]);
+static double signArctan2(double ,double);
+static double normalized(double ,double,double ,double);
+void __mpatan2(mp_no *,mp_no *,mp_no *,int);
+
+double __ieee754_atan2(double y,double x) {
+
+ int i,de,ux,dx,uy,dy,p;
+ static const int pr[MM]={6,8,10,20,32};
+ double ax,ay,u,du,u9,ua,v,vv,dv,t1,t2,t3,t4,t5,t6,t7,t8,
+ z,zz,z1,z2,cor,s1,ss1,s2,ss2;
+ number num;
+ mp_no mperr,mpt1,mpx,mpy,mpz,mpz1,mpz2;
+
+ static const int ep= 59768832, /* 57*16**5 */
+ em=-59768832; /* -57*16**5 */
+
+ /* x=NaN or y=NaN */
+ num.d = x; ux = num.i[HIGH_HALF]; dx = num.i[LOW_HALF];
+ if ((ux&0x7ff00000) ==0x7ff00000) {
+ if (((ux&0x000fffff)|dx)!=0x00000000) return x+x; }
+ num.d = y; uy = num.i[HIGH_HALF]; dy = num.i[LOW_HALF];
+ if ((uy&0x7ff00000) ==0x7ff00000) {
+ if (((uy&0x000fffff)|dy)!=0x00000000) return y+y; }
+
+ /* y=+-0 */
+ if (uy==0x00000000) {
+ if (dy==0x00000000) {
+ if ((ux&0x80000000)==0x00000000) return ZERO;
+ else return opi.d; } }
+ else if (uy==0x80000000) {
+ if (dy==0x00000000) {
+ if ((ux&0x80000000)==0x00000000) return MZERO;
+ else return mopi.d;} }
+
+ /* x=+-0 */
+ if (x==ZERO) {
+ if ((uy&0x80000000)==0x00000000) return hpi.d;
+ else return mhpi.d; }
+
+ /* x=+-INF */
+ if (ux==0x7ff00000) {
+ if (dx==0x00000000) {
+ if (uy==0x7ff00000) {
+ if (dy==0x00000000) return qpi.d; }
+ else if (uy==0xfff00000) {
+ if (dy==0x00000000) return mqpi.d; }
+ else {
+ if ((uy&0x80000000)==0x00000000) return ZERO;
+ else return MZERO; }
+ }
+ }
+ else if (ux==0xfff00000) {
+ if (dx==0x00000000) {
+ if (uy==0x7ff00000) {
+ if (dy==0x00000000) return tqpi.d; }
+ else if (uy==0xfff00000) {
+ if (dy==0x00000000) return mtqpi.d; }
+ else {
+ if ((uy&0x80000000)==0x00000000) return opi.d;
+ else return mopi.d; }
+ }
+ }
+
+ /* y=+-INF */
+ if (uy==0x7ff00000) {
+ if (dy==0x00000000) return hpi.d; }
+ else if (uy==0xfff00000) {
+ if (dy==0x00000000) return mhpi.d; }
+
+ /* either x/y or y/x is very close to zero */
+ ax = (x<ZERO) ? -x : x; ay = (y<ZERO) ? -y : y;
+ de = (uy & 0x7ff00000) - (ux & 0x7ff00000);
+ if (de>=ep) { return ((y>ZERO) ? hpi.d : mhpi.d); }
+ else if (de<=em) {
+ if (x>ZERO) {
+ if ((z=ay/ax)<TWOM1022) return normalized(ax,ay,y,z);
+ else return signArctan2(y,z); }
+ else { return ((y>ZERO) ? opi.d : mopi.d); } }
+
+ /* if either x or y is extremely close to zero, scale abs(x), abs(y). */
+ if (ax<twom500.d || ay<twom500.d) { ax*=two500.d; ay*=two500.d; }
+
+ /* x,y which are neither special nor extreme */
+ if (ay<ax) {
+ u=ay/ax;
+ EMULV(ax,u,v,vv,t1,t2,t3,t4,t5)
+ du=((ay-v)-vv)/ax; }
+ else {
+ u=ax/ay;
+ EMULV(ay,u,v,vv,t1,t2,t3,t4,t5)
+ du=((ax-v)-vv)/ay; }
-#include "math.h"
-#include "math_private.h"
-
-#ifdef __STDC__
-static const double
-#else
-static double
-#endif
-tiny = 1.0e-300,
-zero = 0.0,
-pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */
-pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */
-pi = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */
-pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
-
-#ifdef __STDC__
- double __ieee754_atan2(double y, double x)
-#else
- double __ieee754_atan2(y,x)
- double y,x;
-#endif
-{
- double z;
- int32_t k,m,hx,hy,ix,iy;
- u_int32_t lx,ly;
-
- EXTRACT_WORDS(hx,lx,x);
- ix = hx&0x7fffffff;
- EXTRACT_WORDS(hy,ly,y);
- iy = hy&0x7fffffff;
- if(((ix|((lx|-lx)>>31))>0x7ff00000)||
- ((iy|((ly|-ly)>>31))>0x7ff00000)) /* x or y is NaN */
- return x+y;
- if(((hx-0x3ff00000)|lx)==0) return __atan(y); /* x=1.0 */
- m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */
-
- /* when y = 0 */
- if((iy|ly)==0) {
- switch(m) {
- case 0:
- case 1: return y; /* atan(+-0,+anything)=+-0 */
- case 2: return pi+tiny;/* atan(+0,-anything) = pi */
- case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
- }
- }
- /* when x = 0 */
- if((ix|lx)==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
-
- /* when x is INF */
- if(ix==0x7ff00000) {
- if(iy==0x7ff00000) {
- switch(m) {
- case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */
- case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
- case 2: return 3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/
- case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/
- }
- } else {
- switch(m) {
- case 0: return zero ; /* atan(+...,+INF) */
- case 1: return -zero ; /* atan(-...,+INF) */
- case 2: return pi+tiny ; /* atan(+...,-INF) */
- case 3: return -pi-tiny ; /* atan(-...,-INF) */
- }
- }
- }
- /* when y is INF */
- if(iy==0x7ff00000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
-
- /* compute y/x */
- k = (iy-ix)>>20;
- if(k > 60) z=pi_o_2+0.5*pi_lo; /* |y/x| > 2**60 */
- else if(hx<0&&k<-60) z=0.0; /* |y|/x < -2**60 */
- else z=__atan(fabs(y/x)); /* safe to do y/x */
- switch (m) {
- case 0: return z ; /* atan(+,+) */
- case 1: {
- u_int32_t zh;
- GET_HIGH_WORD(zh,z);
- SET_HIGH_WORD(z,zh ^ 0x80000000);
- }
- return z ; /* atan(-,+) */
- case 2: return pi-(z-pi_lo);/* atan(+,-) */
- default: /* case 3 */
- return (z-pi_lo)-pi;/* atan(-,-) */
- }
+ if (x>ZERO) {
+
+ /* (i) x>0, abs(y)< abs(x): atan(ay/ax) */
+ if (ay<ax) {
+ if (u<inv16.d) {
+ v=u*u; zz=du+u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
+ if ((z=u+(zz-u1.d*u)) == u+(zz+u1.d*u)) return signArctan2(y,z);
+
+ MUL2(u,du,u,du,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(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
+ if ((z=s1+(ss1-u5.d*s1)) == s1+(ss1+u5.d*s1)) return signArctan2(y,z);
+ return atan2Mp(x,y,pr);
+ }
+ else {
+ i=(TWO52+TWO8*u)-TWO52; i-=16;
+ t3=u-cij[i][0].d;
+ EADD(t3,du,v,dv)
+ t1=cij[i][1].d; t2=cij[i][2].d;
+ zz=v*t2+(dv*t2+v*v*(cij[i][3].d+v*(cij[i][4].d+
+ v*(cij[i][5].d+v* cij[i][6].d))));
+ if (i<112) {
+ if (i<48) u9=u91.d; /* u < 1/4 */
+ else u9=u92.d; } /* 1/4 <= u < 1/2 */
+ else {
+ if (i<176) u9=u93.d; /* 1/2 <= u < 3/4 */
+ else u9=u94.d; } /* 3/4 <= u <= 1 */
+ if ((z=t1+(zz-u9*t1)) == t1+(zz+u9*t1)) return signArctan2(y,z);
+
+ t1=u-hij[i][0].d;
+ EADD(t1,du,v,vv)
+ s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
+ v*(hij[i][14].d+v* hij[i][15].d))));
+ ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
+ MUL2(v,vv,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(v,vv,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(v,vv,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(v,vv,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 ((z=s2+(ss2-ub.d*s2)) == s2+(ss2+ub.d*s2)) return signArctan2(y,z);
+ return atan2Mp(x,y,pr);
+ }
+ }
+
+ /* (ii) x>0, abs(x)<=abs(y): pi/2-atan(ax/ay) */
+ else {
+ if (u<inv16.d) {
+ v=u*u;
+ zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
+ ESUB(hpi.d,u,t2,cor)
+ t3=((hpi1.d+cor)-du)-zz;
+ if ((z=t2+(t3-u2.d)) == t2+(t3+u2.d)) return signArctan2(y,z);
+
+ MUL2(u,du,u,du,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(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
+ SUB2(hpi.d,hpi1.d,s1,ss1,s2,ss2,t1,t2)
+ if ((z=s2+(ss2-u6.d)) == s2+(ss2+u6.d)) return signArctan2(y,z);
+ return atan2Mp(x,y,pr);
+ }
+ else {
+ i=(TWO52+TWO8*u)-TWO52; i-=16;
+ v=(u-cij[i][0].d)+du;
+ zz=hpi1.d-v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
+ v*(cij[i][5].d+v* cij[i][6].d))));
+ t1=hpi.d-cij[i][1].d;
+ if (i<112) ua=ua1.d; /* w < 1/2 */
+ else ua=ua2.d; /* w >= 1/2 */
+ if ((z=t1+(zz-ua)) == t1+(zz+ua)) return signArctan2(y,z);
+
+ t1=u-hij[i][0].d;
+ EADD(t1,du,v,vv)
+ s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
+ v*(hij[i][14].d+v* hij[i][15].d))));
+ ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
+ MUL2(v,vv,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(v,vv,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(v,vv,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(v,vv,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.d,hpi1.d,s2,ss2,s1,ss1,t1,t2)
+ if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d)) return signArctan2(y,z);
+ return atan2Mp(x,y,pr);
+ }
+ }
+ }
+ else {
+
+ /* (iii) x<0, abs(x)< abs(y): pi/2+atan(ax/ay) */
+ if (ax<ay) {
+ if (u<inv16.d) {
+ v=u*u;
+ zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
+ EADD(hpi.d,u,t2,cor)
+ t3=((hpi1.d+cor)+du)+zz;
+ if ((z=t2+(t3-u3.d)) == t2+(t3+u3.d)) return signArctan2(y,z);
+
+ MUL2(u,du,u,du,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(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
+ ADD2(hpi.d,hpi1.d,s1,ss1,s2,ss2,t1,t2)
+ if ((z=s2+(ss2-u7.d)) == s2+(ss2+u7.d)) return signArctan2(y,z);
+ return atan2Mp(x,y,pr);
+ }
+ else {
+ i=(TWO52+TWO8*u)-TWO52; i-=16;
+ v=(u-cij[i][0].d)+du;
+ zz=hpi1.d+v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
+ v*(cij[i][5].d+v* cij[i][6].d))));
+ t1=hpi.d+cij[i][1].d;
+ if (i<112) ua=ua1.d; /* w < 1/2 */
+ else ua=ua2.d; /* w >= 1/2 */
+ if ((z=t1+(zz-ua)) == t1+(zz+ua)) return signArctan2(y,z);
+
+ t1=u-hij[i][0].d;
+ EADD(t1,du,v,vv)
+ s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
+ v*(hij[i][14].d+v* hij[i][15].d))));
+ ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
+ MUL2(v,vv,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(v,vv,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(v,vv,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(v,vv,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)
+ ADD2(hpi.d,hpi1.d,s2,ss2,s1,ss1,t1,t2)
+ if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d)) return signArctan2(y,z);
+ return atan2Mp(x,y,pr);
+ }
+ }
+
+ /* (iv) x<0, abs(y)<=abs(x): pi-atan(ax/ay) */
+ else {
+ if (u<inv16.d) {
+ v=u*u;
+ zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
+ ESUB(opi.d,u,t2,cor)
+ t3=((opi1.d+cor)-du)-zz;
+ if ((z=t2+(t3-u4.d)) == t2+(t3+u4.d)) return signArctan2(y,z);
+
+ MUL2(u,du,u,du,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(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
+ SUB2(opi.d,opi1.d,s1,ss1,s2,ss2,t1,t2)
+ if ((z=s2+(ss2-u8.d)) == s2+(ss2+u8.d)) return signArctan2(y,z);
+ return atan2Mp(x,y,pr);
+ }
+ else {
+ i=(TWO52+TWO8*u)-TWO52; i-=16;
+ v=(u-cij[i][0].d)+du;
+ zz=opi1.d-v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
+ v*(cij[i][5].d+v* cij[i][6].d))));
+ t1=opi.d-cij[i][1].d;
+ if (i<112) ua=ua1.d; /* w < 1/2 */
+ else ua=ua2.d; /* w >= 1/2 */
+ if ((z=t1+(zz-ua)) == t1+(zz+ua)) return signArctan2(y,z);
+
+ t1=u-hij[i][0].d;
+ EADD(t1,du,v,vv)
+ s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
+ v*(hij[i][14].d+v* hij[i][15].d))));
+ ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
+ MUL2(v,vv,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(v,vv,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(v,vv,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(v,vv,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(opi.d,opi1.d,s2,ss2,s1,ss1,t1,t2)
+ if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d)) return signArctan2(y,z);
+ return atan2Mp(x,y,pr);
+ }
+ }
+ }
+}
+ /* Treat the Denormalized case */
+static double normalized(double ax,double ay,double y, double z)
+ { int p;
+ mp_no mpx,mpy,mpz,mperr,mpz2,mpt1;
+ p=6;
+ dbl_mp(ax,&mpx,p); dbl_mp(ay,&mpy,p); dvd(&mpy,&mpx,&mpz,p);
+ dbl_mp(ue.d,&mpt1,p); mul(&mpz,&mpt1,&mperr,p);
+ sub(&mpz,&mperr,&mpz2,p); mp_dbl(&mpz2,&z,p);
+ return signArctan2(y,z);
+}
+ /* Fix the sign and return after stage 1 or stage 2 */
+static double signArctan2(double y,double z)
+{
+ return ((y<ZERO) ? -z : z);
+}
+ /* Stage 3: Perform a multi-Precision computation */
+static double atan2Mp(double x,double y,const int pr[])
+{
+ double z1,z2;
+ int i,p;
+ mp_no mpx,mpy,mpz,mpz1,mpz2,mperr,mpt1;
+ for (i=0; i<MM; i++) {
+ p = pr[i];
+ dbl_mp(x,&mpx,p); dbl_mp(y,&mpy,p);
+ __mpatan2(&mpy,&mpx,&mpz,p);
+ dbl_mp(ud[i].d,&mpt1,p); mul(&mpz,&mpt1,&mperr,p);
+ add(&mpz,&mperr,&mpz1,p); sub(&mpz,&mperr,&mpz2,p);
+ mp_dbl(&mpz1,&z1,p); mp_dbl(&mpz2,&z2,p);
+ if (z1==z2) return z1;
+ }
+ return z1; /*if unpossible to do exact computing */
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-07 02:53:29 +0000