Back to index

glibc  2.9
Functions
e_atan2.c File Reference
#include "dla.h"
#include "mpa.h"
#include "MathLib.h"
#include "uatan.tbl"
#include "atnat2.h"
#include "math_private.h"

Go to the source code of this file.

Functions

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)

Function Documentation

double __ieee754_atan2 ( double  y,
double  x 
)

Definition at line 58 of file e_atan2.c.

                                          {

  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);
      }
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void __mpatan2 ( mp_no ,
mp_no ,
mp_no ,
int   
)

Definition at line 46 of file mpatan2.c.

                                                    {

  static const double ZERO = 0.0, ONE = 1.0;

  mp_no mpone = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
                    0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
                    0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
  mp_no mpt1,mpt2,mpt3;


  if (X[0] <= ZERO) {
    mpone.e = 1;                 mpone.d[0] = mpone.d[1] = ONE;
    __dvd(x,y,&mpt1,p);          __mul(&mpt1,&mpt1,&mpt2,p);
    if (mpt1.d[0] != ZERO)       mpt1.d[0] = ONE;
    __add(&mpt2,&mpone,&mpt3,p); __mpsqrt(&mpt3,&mpt2,p);
    __add(&mpt1,&mpt2,&mpt3,p);  mpt3.d[0]=Y[0];
    __mpatan(&mpt3,&mpt1,p);     __add(&mpt1,&mpt1,z,p);
  }
  else
  { __dvd(y,x,&mpt1,p);
    __mpatan(&mpt1,z,p);
  }

  return;
}

Here is the caller graph for this function:

static double atan2Mp ( double  x,
double  y,
const int  pr[] 
) [static]

Definition at line 391 of file e_atan2.c.

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

Here is the call graph for this function:

Here is the caller graph for this function:

static double normalized ( double  ax,
double  ay,
double  y,
double  z 
) [static]

Definition at line 376 of file e_atan2.c.

    { 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);
}

Here is the call graph for this function:

Here is the caller graph for this function:

static double signArctan2 ( double  y,
double  z 
) [static]

Definition at line 386 of file e_atan2.c.

{
  return ((y<ZERO) ? -z : z);
}

Here is the caller graph for this function: