Back to index

glibc  2.9
e_atan2.c
Go to the documentation of this file.
00001 /*
00002  * IBM Accurate Mathematical Library
00003  * written by International Business Machines Corp.
00004  * Copyright (C) 2001 Free Software Foundation
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU Lesser General Public License as published by
00008  * the Free Software Foundation; either version 2.1 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU Lesser General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU Lesser General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00019  */
00020 /************************************************************************/
00021 /*  MODULE_NAME: atnat2.c                                               */
00022 /*                                                                      */
00023 /*  FUNCTIONS: uatan2                                                   */
00024 /*             atan2Mp                                                  */
00025 /*             signArctan2                                              */
00026 /*             normalized                                               */
00027 /*                                                                      */
00028 /*  FILES NEEDED: dla.h endian.h mpa.h mydefs.h atnat2.h                */
00029 /*                mpatan.c mpatan2.c mpsqrt.c                           */
00030 /*                uatan.tbl                                             */
00031 /*                                                                      */
00032 /* An ultimate atan2() routine. Given two IEEE double machine numbers y,*/
00033 /* x it computes the correctly rounded (to nearest) value of atan2(y,x).*/
00034 /*                                                                      */
00035 /* Assumption: Machine arithmetic operations are performed in           */
00036 /* round to nearest mode of IEEE 754 standard.                          */
00037 /*                                                                      */
00038 /************************************************************************/
00039 
00040 #include "dla.h"
00041 #include "mpa.h"
00042 #include "MathLib.h"
00043 #include "uatan.tbl"
00044 #include "atnat2.h"
00045 #include "math_private.h"
00046 
00047 /************************************************************************/
00048 /* An ultimate atan2 routine. Given two IEEE double machine numbers y,x */
00049 /* it computes the correctly rounded (to nearest) value of atan2(y,x).  */
00050 /* Assumption: Machine arithmetic operations are performed in           */
00051 /* round to nearest mode of IEEE 754 standard.                          */
00052 /************************************************************************/
00053 static double atan2Mp(double ,double ,const int[]);
00054 static double signArctan2(double ,double);
00055 static double normalized(double ,double,double ,double);
00056 void __mpatan2(mp_no *,mp_no *,mp_no *,int);
00057 
00058 double __ieee754_atan2(double y,double x) {
00059 
00060   int i,de,ux,dx,uy,dy;
00061 #if 0
00062   int p;
00063 #endif
00064   static const int pr[MM]={6,8,10,20,32};
00065   double ax,ay,u,du,u9,ua,v,vv,dv,t1,t2,t3,t4,t5,t6,t7,t8,
00066          z,zz,cor,s1,ss1,s2,ss2;
00067 #if 0
00068   double z1,z2;
00069 #endif
00070   number num;
00071 #if 0
00072   mp_no mperr,mpt1,mpx,mpy,mpz,mpz1,mpz2;
00073 #endif
00074 
00075   static const int ep= 59768832,   /*  57*16**5   */
00076                    em=-59768832;   /* -57*16**5   */
00077 
00078   /* x=NaN or y=NaN */
00079   num.d = x;  ux = num.i[HIGH_HALF];  dx = num.i[LOW_HALF];
00080   if   ((ux&0x7ff00000)    ==0x7ff00000) {
00081     if (((ux&0x000fffff)|dx)!=0x00000000) return x+x; }
00082   num.d = y;  uy = num.i[HIGH_HALF];  dy = num.i[LOW_HALF];
00083   if   ((uy&0x7ff00000)    ==0x7ff00000) {
00084     if (((uy&0x000fffff)|dy)!=0x00000000) return y+y; }
00085 
00086   /* y=+-0 */
00087   if      (uy==0x00000000) {
00088     if    (dy==0x00000000) {
00089       if  ((ux&0x80000000)==0x00000000)  return ZERO;
00090       else                               return opi.d; } }
00091   else if (uy==0x80000000) {
00092     if    (dy==0x00000000) {
00093       if  ((ux&0x80000000)==0x00000000)  return MZERO;
00094       else                               return mopi.d;} }
00095 
00096   /* x=+-0 */
00097   if (x==ZERO) {
00098     if ((uy&0x80000000)==0x00000000)     return hpi.d;
00099     else                                 return mhpi.d; }
00100 
00101   /* x=+-INF */
00102   if          (ux==0x7ff00000) {
00103     if        (dx==0x00000000) {
00104       if      (uy==0x7ff00000) {
00105         if    (dy==0x00000000)  return qpi.d; }
00106       else if (uy==0xfff00000) {
00107         if    (dy==0x00000000)  return mqpi.d; }
00108       else {
00109         if    ((uy&0x80000000)==0x00000000)  return ZERO;
00110         else                                 return MZERO; }
00111     }
00112   }
00113   else if     (ux==0xfff00000) {
00114     if        (dx==0x00000000) {
00115       if      (uy==0x7ff00000) {
00116         if    (dy==0x00000000)  return tqpi.d; }
00117       else if (uy==0xfff00000) {
00118         if    (dy==0x00000000)  return mtqpi.d; }
00119       else                     {
00120         if    ((uy&0x80000000)==0x00000000)  return opi.d;
00121         else                                 return mopi.d; }
00122     }
00123   }
00124 
00125   /* y=+-INF */
00126   if      (uy==0x7ff00000) {
00127     if    (dy==0x00000000)  return hpi.d; }
00128   else if (uy==0xfff00000) {
00129     if    (dy==0x00000000)  return mhpi.d; }
00130 
00131   /* either x/y or y/x is very close to zero */
00132   ax = (x<ZERO) ? -x : x;    ay = (y<ZERO) ? -y : y;
00133   de = (uy & 0x7ff00000) - (ux & 0x7ff00000);
00134   if      (de>=ep)  { return ((y>ZERO) ? hpi.d : mhpi.d); }
00135   else if (de<=em)  {
00136     if    (x>ZERO)  {
00137       if  ((z=ay/ax)<TWOM1022)  return normalized(ax,ay,y,z);
00138       else                      return signArctan2(y,z); }
00139     else            { return ((y>ZERO) ? opi.d : mopi.d); } }
00140 
00141   /* if either x or y is extremely close to zero, scale abs(x), abs(y). */
00142   if (ax<twom500.d || ay<twom500.d) { ax*=two500.d;  ay*=two500.d; }
00143 
00144   /* x,y which are neither special nor extreme */
00145   if (ay<ax) {
00146     u=ay/ax;
00147     EMULV(ax,u,v,vv,t1,t2,t3,t4,t5)
00148     du=((ay-v)-vv)/ax; }
00149   else {
00150     u=ax/ay;
00151     EMULV(ay,u,v,vv,t1,t2,t3,t4,t5)
00152     du=((ax-v)-vv)/ay; }
00153 
00154   if (x>ZERO) {
00155 
00156     /* (i)   x>0, abs(y)< abs(x):  atan(ay/ax) */
00157     if (ay<ax) {
00158       if (u<inv16.d) {
00159         v=u*u;  zz=du+u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
00160         if ((z=u+(zz-u1.d*u)) == u+(zz+u1.d*u))  return signArctan2(y,z);
00161 
00162         MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
00163         s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
00164         ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
00165         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00166         ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
00167         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00168         ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
00169         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00170         ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
00171         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00172         MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
00173         ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
00174         if ((z=s1+(ss1-u5.d*s1)) == s1+(ss1+u5.d*s1))  return signArctan2(y,z);
00175         return atan2Mp(x,y,pr);
00176       }
00177       else {
00178         i=(TWO52+TWO8*u)-TWO52;  i-=16;
00179         t3=u-cij[i][0].d;
00180         EADD(t3,du,v,dv)
00181         t1=cij[i][1].d;  t2=cij[i][2].d;
00182         zz=v*t2+(dv*t2+v*v*(cij[i][3].d+v*(cij[i][4].d+
00183                          v*(cij[i][5].d+v* cij[i][6].d))));
00184         if (i<112) {
00185           if (i<48)  u9=u91.d;    /* u < 1/4        */
00186           else       u9=u92.d; }  /* 1/4 <= u < 1/2 */
00187         else {
00188           if (i<176) u9=u93.d;    /* 1/2 <= u < 3/4 */
00189           else       u9=u94.d; }  /* 3/4 <= u <= 1  */
00190         if ((z=t1+(zz-u9*t1)) == t1+(zz+u9*t1))  return signArctan2(y,z);
00191 
00192         t1=u-hij[i][0].d;
00193         EADD(t1,du,v,vv)
00194         s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
00195            v*(hij[i][14].d+v* hij[i][15].d))));
00196         ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
00197         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00198         ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
00199         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00200         ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
00201         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00202         ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
00203         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00204         ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
00205         if ((z=s2+(ss2-ub.d*s2)) == s2+(ss2+ub.d*s2))  return signArctan2(y,z);
00206         return atan2Mp(x,y,pr);
00207       }
00208     }
00209 
00210     /* (ii)  x>0, abs(x)<=abs(y):  pi/2-atan(ax/ay) */
00211     else {
00212       if (u<inv16.d) {
00213         v=u*u;
00214         zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
00215         ESUB(hpi.d,u,t2,cor)
00216         t3=((hpi1.d+cor)-du)-zz;
00217         if ((z=t2+(t3-u2.d)) == t2+(t3+u2.d))  return signArctan2(y,z);
00218 
00219         MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
00220         s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
00221         ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
00222         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00223         ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
00224         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00225         ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
00226         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00227         ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
00228         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00229         MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
00230         ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
00231         SUB2(hpi.d,hpi1.d,s1,ss1,s2,ss2,t1,t2)
00232         if ((z=s2+(ss2-u6.d)) == s2+(ss2+u6.d))  return signArctan2(y,z);
00233         return atan2Mp(x,y,pr);
00234       }
00235       else {
00236         i=(TWO52+TWO8*u)-TWO52;  i-=16;
00237         v=(u-cij[i][0].d)+du;
00238         zz=hpi1.d-v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
00239                                  v*(cij[i][5].d+v* cij[i][6].d))));
00240         t1=hpi.d-cij[i][1].d;
00241         if (i<112)  ua=ua1.d;  /* w <  1/2 */
00242         else        ua=ua2.d;  /* w >= 1/2 */
00243         if ((z=t1+(zz-ua)) == t1+(zz+ua))  return signArctan2(y,z);
00244 
00245         t1=u-hij[i][0].d;
00246         EADD(t1,du,v,vv)
00247         s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
00248            v*(hij[i][14].d+v* hij[i][15].d))));
00249         ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
00250         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00251         ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
00252         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00253         ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
00254         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00255         ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
00256         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00257         ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
00258         SUB2(hpi.d,hpi1.d,s2,ss2,s1,ss1,t1,t2)
00259         if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d))  return signArctan2(y,z);
00260         return atan2Mp(x,y,pr);
00261       }
00262     }
00263   }
00264   else {
00265 
00266     /* (iii) x<0, abs(x)< abs(y):  pi/2+atan(ax/ay) */
00267     if (ax<ay) {
00268       if (u<inv16.d) {
00269         v=u*u;
00270         zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
00271         EADD(hpi.d,u,t2,cor)
00272         t3=((hpi1.d+cor)+du)+zz;
00273         if ((z=t2+(t3-u3.d)) == t2+(t3+u3.d))  return signArctan2(y,z);
00274 
00275         MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
00276         s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
00277         ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
00278         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00279         ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
00280         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00281         ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
00282         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00283         ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
00284         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00285         MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
00286         ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
00287         ADD2(hpi.d,hpi1.d,s1,ss1,s2,ss2,t1,t2)
00288         if ((z=s2+(ss2-u7.d)) == s2+(ss2+u7.d))  return signArctan2(y,z);
00289         return atan2Mp(x,y,pr);
00290       }
00291       else {
00292         i=(TWO52+TWO8*u)-TWO52;  i-=16;
00293         v=(u-cij[i][0].d)+du;
00294         zz=hpi1.d+v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
00295                                  v*(cij[i][5].d+v* cij[i][6].d))));
00296         t1=hpi.d+cij[i][1].d;
00297         if (i<112)  ua=ua1.d;  /* w <  1/2 */
00298         else        ua=ua2.d;  /* w >= 1/2 */
00299         if ((z=t1+(zz-ua)) == t1+(zz+ua))  return signArctan2(y,z);
00300 
00301         t1=u-hij[i][0].d;
00302         EADD(t1,du,v,vv)
00303         s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
00304            v*(hij[i][14].d+v* hij[i][15].d))));
00305         ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
00306         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00307         ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
00308         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00309         ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
00310         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00311         ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
00312         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00313         ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
00314         ADD2(hpi.d,hpi1.d,s2,ss2,s1,ss1,t1,t2)
00315         if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d))  return signArctan2(y,z);
00316         return atan2Mp(x,y,pr);
00317       }
00318     }
00319 
00320     /* (iv)  x<0, abs(y)<=abs(x):  pi-atan(ax/ay) */
00321     else {
00322       if (u<inv16.d) {
00323         v=u*u;
00324         zz=u*v*(d3.d+v*(d5.d+v*(d7.d+v*(d9.d+v*(d11.d+v*d13.d)))));
00325         ESUB(opi.d,u,t2,cor)
00326         t3=((opi1.d+cor)-du)-zz;
00327         if ((z=t2+(t3-u4.d)) == t2+(t3+u4.d))  return signArctan2(y,z);
00328 
00329         MUL2(u,du,u,du,v,vv,t1,t2,t3,t4,t5,t6,t7,t8)
00330         s1=v*(f11.d+v*(f13.d+v*(f15.d+v*(f17.d+v*f19.d))));
00331         ADD2(f9.d,ff9.d,s1,ZERO,s2,ss2,t1,t2)
00332         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00333         ADD2(f7.d,ff7.d,s1,ss1,s2,ss2,t1,t2)
00334         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00335         ADD2(f5.d,ff5.d,s1,ss1,s2,ss2,t1,t2)
00336         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00337         ADD2(f3.d,ff3.d,s1,ss1,s2,ss2,t1,t2)
00338         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00339         MUL2(u,du,s1,ss1,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8)
00340         ADD2(u,du,s2,ss2,s1,ss1,t1,t2)
00341         SUB2(opi.d,opi1.d,s1,ss1,s2,ss2,t1,t2)
00342         if ((z=s2+(ss2-u8.d)) == s2+(ss2+u8.d))  return signArctan2(y,z);
00343         return atan2Mp(x,y,pr);
00344       }
00345       else {
00346         i=(TWO52+TWO8*u)-TWO52;  i-=16;
00347         v=(u-cij[i][0].d)+du;
00348         zz=opi1.d-v*(cij[i][2].d+v*(cij[i][3].d+v*(cij[i][4].d+
00349                                  v*(cij[i][5].d+v* cij[i][6].d))));
00350         t1=opi.d-cij[i][1].d;
00351         if (i<112)  ua=ua1.d;  /* w <  1/2 */
00352         else        ua=ua2.d;  /* w >= 1/2 */
00353         if ((z=t1+(zz-ua)) == t1+(zz+ua))  return signArctan2(y,z);
00354 
00355         t1=u-hij[i][0].d;
00356         EADD(t1,du,v,vv)
00357         s1=v*(hij[i][11].d+v*(hij[i][12].d+v*(hij[i][13].d+
00358            v*(hij[i][14].d+v* hij[i][15].d))));
00359         ADD2(hij[i][9].d,hij[i][10].d,s1,ZERO,s2,ss2,t1,t2)
00360         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00361         ADD2(hij[i][7].d,hij[i][8].d,s1,ss1,s2,ss2,t1,t2)
00362         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00363         ADD2(hij[i][5].d,hij[i][6].d,s1,ss1,s2,ss2,t1,t2)
00364         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00365         ADD2(hij[i][3].d,hij[i][4].d,s1,ss1,s2,ss2,t1,t2)
00366         MUL2(v,vv,s2,ss2,s1,ss1,t1,t2,t3,t4,t5,t6,t7,t8)
00367         ADD2(hij[i][1].d,hij[i][2].d,s1,ss1,s2,ss2,t1,t2)
00368         SUB2(opi.d,opi1.d,s2,ss2,s1,ss1,t1,t2)
00369         if ((z=s1+(ss1-uc.d)) == s1+(ss1+uc.d))  return signArctan2(y,z);
00370         return atan2Mp(x,y,pr);
00371       }
00372     }
00373   }
00374 }
00375   /* Treat the Denormalized case */
00376 static double  normalized(double ax,double ay,double y, double z)
00377     { int p;
00378       mp_no mpx,mpy,mpz,mperr,mpz2,mpt1;
00379   p=6;
00380   __dbl_mp(ax,&mpx,p);  __dbl_mp(ay,&mpy,p);  __dvd(&mpy,&mpx,&mpz,p);
00381   __dbl_mp(ue.d,&mpt1,p);   __mul(&mpz,&mpt1,&mperr,p);
00382   __sub(&mpz,&mperr,&mpz2,p);  __mp_dbl(&mpz2,&z,p);
00383   return signArctan2(y,z);
00384 }
00385   /* Fix the sign and return after stage 1 or stage 2 */
00386 static double signArctan2(double y,double z)
00387 {
00388   return ((y<ZERO) ? -z : z);
00389 }
00390   /* Stage 3: Perform a multi-Precision computation */
00391 static double  atan2Mp(double x,double y,const int pr[])
00392 {
00393   double z1,z2;
00394   int i,p;
00395   mp_no mpx,mpy,mpz,mpz1,mpz2,mperr,mpt1;
00396   for (i=0; i<MM; i++) {
00397     p = pr[i];
00398     __dbl_mp(x,&mpx,p);  __dbl_mp(y,&mpy,p);
00399     __mpatan2(&mpy,&mpx,&mpz,p);
00400     __dbl_mp(ud[i].d,&mpt1,p);   __mul(&mpz,&mpt1,&mperr,p);
00401     __add(&mpz,&mperr,&mpz1,p);  __sub(&mpz,&mperr,&mpz2,p);
00402     __mp_dbl(&mpz1,&z1,p);       __mp_dbl(&mpz2,&z2,p);
00403     if (z1==z2)   return z1;
00404   }
00405   return z1; /*if unpossible to do exact computing */
00406 }