📄 e_atan2.cpp
字号:
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 */
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -