mirror of
https://sourceware.org/git/glibc.git
synced 2024-11-23 05:20:06 +00:00
c6c6dd4803
* manual/contrib.texi: Removed licenses, added acknowledgements
for contributions by Intel, IBM, Craig Metz.
* LICENSES: New file, contains the text of all non-FSF licenses in the
distribution that require putting the notice in the accompanying
documentation.
* README.template, README: Mention LICENSES.
* sysdeps/mach/hurd/net/if_ppp.h: Replaced CMU license with a
new one modelled on the modern BSD license, per recent letter
of permission from CMU.
* sysdeps/unix/sysv/linux/net/if_ppp.h: Likewise.
* sysdeps/ieee754/dbl-64/MathLib.h: Changed the copyright holder
from IBM to FSF, per the recent Software Letter. Changed the
distribution terms from GPL to LGPL.
* sysdeps/ieee754/dbl-64/asincos.tbl: Added FSF copyright and
copying permission notice (Lesser GPL), per recent IBM Software Letter.
* sysdeps/ieee754/dbl-64/powtwo.tbl: Likewise.
* sysdeps/ieee754/dbl-64/root.tbl: Likewise.
* sysdeps/ieee754/dbl-64/sincos.tbl: Likewise.
* sysdeps/ieee754/dbl-64/uatan.tbl: Likewise.
* sysdeps/ieee754/dbl-64/uexp.tbl: Likewise.
* sysdeps/ieee754/dbl-64/ulog.tbl: Likewise.
* sysdeps/ieee754/dbl-64/upow.tbl: Likewise.
* sysdeps/ieee754/dbl-64/utan.tbl: Likewise.
* sysdeps/ieee754/dbl-64/atnat.h: Changed the copyright holder
from IBM to FSF, per the recent Software Letter. Corrected the
text of the copying permission notice to say Lesser GPL instead
of GPL in warranty disclaimer paragraph.
* sysdeps/ieee754/dbl-64/atnat2.h: Likewise.
* sysdeps/ieee754/dbl-64/branred.h: Likewise.
* sysdeps/ieee754/dbl-64/dla.h: Likewise.
* sysdeps/ieee754/dbl-64/doasin.h: Likewise.
* sysdeps/ieee754/dbl-64/dosincos.h: Likewise.
* sysdeps/ieee754/dbl-64/mpa.h: Likewise.
* sysdeps/ieee754/dbl-64/mpa2.h: Likewise.
* sysdeps/ieee754/dbl-64/mpatan.h: Likewise.
* sysdeps/ieee754/dbl-64/mpexp.h: Likewise.
* sysdeps/ieee754/dbl-64/mplog.h: Likewise.
* sysdeps/ieee754/dbl-64/mpsqrt.h: Likewise.
* sysdeps/ieee754/dbl-64/mydefs.h: Likewise.
* sysdeps/ieee754/dbl-64/sincos32.h: Likewise.
* sysdeps/ieee754/dbl-64/uasncs.h: Likewise.
* sysdeps/ieee754/dbl-64/uexp.h: Likewise.
* sysdeps/ieee754/dbl-64/ulog.h: Likewise.
* sysdeps/ieee754/dbl-64/upow.h: Likewise.
* sysdeps/ieee754/dbl-64/urem.h: Likewise.
* sysdeps/ieee754/dbl-64/uroot.h: Likewise.
* sysdeps/ieee754/dbl-64/usncs.h: Likewise.
* sysdeps/ieee754/dbl-64/utan.h: Likewise.
* sysdeps/ieee754/dbl-64/branred.c: Corrected the text of the copying
permission notice to say Lesser GPL instead of GPL in warranty
disclaimer paragraph.
* sysdeps/ieee754/dbl-64/doasin.c: Likewise.
* sysdeps/ieee754/dbl-64/dosincos.c: Likewise.
* sysdeps/ieee754/dbl-64/e_asin.c: Likewise.
* sysdeps/ieee754/dbl-64/e_atan2.c: Likewise.
* sysdeps/ieee754/dbl-64/e_exp.c: Likewise.
* sysdeps/ieee754/dbl-64/e_log.c: Likewise.
* sysdeps/ieee754/dbl-64/e_pow.c: Likewise.
* sysdeps/ieee754/dbl-64/e_remainder.c: Likewise.
* sysdeps/ieee754/dbl-64/e_sqrt.c: Likewise.
* sysdeps/ieee754/dbl-64/halfulp.c: Likewise.
* sysdeps/ieee754/dbl-64/mpa.c: Likewise.
* sysdeps/ieee754/dbl-64/mpatan.c: Likewise.
* sysdeps/ieee754/dbl-64/mpatan2.c: Likewise.
* sysdeps/ieee754/dbl-64/mpexp.c: Likewise.
* sysdeps/ieee754/dbl-64/mplog.c: Likewise.
* sysdeps/ieee754/dbl-64/mpsqrt.c: Likewise.
* sysdeps/ieee754/dbl-64/mptan.c: Likewise.
* sysdeps/ieee754/dbl-64/s_atan.c: Likewise.
* sysdeps/ieee754/dbl-64/s_sin.c: Likewise.
* sysdeps/ieee754/dbl-64/s_tan.c: Likewise.
* sysdeps/ieee754/dbl-64/sincos32.c: Likewise.
* sysdeps/ieee754/dbl-64/slowexp.c: Likewise.
* sysdeps/ieee754/dbl-64/slowpow.c: Likewise.
407 lines
16 KiB
C
407 lines
16 KiB
C
/*
|
|
* IBM Accurate Mathematical Library
|
|
* written by International Business Machines Corp.
|
|
* Copyright (C) 2001 Free Software Foundation
|
|
*
|
|
* 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.
|
|
*
|
|
* 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 Lesser General Public License for more details.
|
|
*
|
|
* 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"
|
|
#include "math_private.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;
|
|
#if 0
|
|
int p;
|
|
#endif
|
|
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,cor,s1,ss1,s2,ss2;
|
|
#if 0
|
|
double z1,z2;
|
|
#endif
|
|
number num;
|
|
#if 0
|
|
mp_no mperr,mpt1,mpx,mpy,mpz,mpz1,mpz2;
|
|
#endif
|
|
|
|
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; }
|
|
|
|
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 */
|
|
}
|