e_atan2.c

Glibc 2.3.2源代码(解压后有100多M)

        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 */
}