00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #if defined(LIBM_SCCS) && !defined(lint)
00014 static char rcsid[] = "$NetBSD: e_atan2.c,v 1.8 1995/05/10 20:44:51 jtc Exp $";
00015 #endif
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 #include "math.h"
00045 #include "mathP.h"
00046
00047 #ifdef __STDC__
00048 static const double
00049 #else
00050 static double
00051 #endif
00052 tiny = 1.0e-300,
00053 zero = 0.0,
00054 pi_o_4 = 7.8539816339744827900E-01,
00055 pi_o_2 = 1.5707963267948965580E+00,
00056 pi = 3.1415926535897931160E+00,
00057 pi_lo = 1.2246467991473531772E-16;
00058
00059 #ifdef __STDC__
00060 double __ieee754_atan2(double y, double x)
00061 #else
00062 double __ieee754_atan2(y,x)
00063 double y,x;
00064 #endif
00065 {
00066 double z;
00067 int32_t k,m,hx,hy,ix,iy;
00068 u_int32_t lx,ly;
00069
00070 EXTRACT_WORDS(hx,lx,x);
00071 ix = hx&0x7fffffff;
00072 EXTRACT_WORDS(hy,ly,y);
00073 iy = hy&0x7fffffff;
00074 if(((ix|((lx|-lx)>>31))>0x7ff00000)||
00075 ((iy|((ly|-ly)>>31))>0x7ff00000))
00076 return x+y;
00077 if((hx-(0x3ff00000|lx))==0) return atan(y);
00078 m = ((hy>>31)&1)|((hx>>30)&2);
00079
00080
00081 if((iy|ly)==0) {
00082 switch(m) {
00083 case 0:
00084 case 1: return y;
00085 case 2: return pi+tiny;
00086 case 3: return -pi-tiny;
00087 }
00088 }
00089
00090 if((ix|lx)==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
00091
00092
00093 if(ix==0x7ff00000) {
00094 if(iy==0x7ff00000) {
00095 switch(m) {
00096 case 0: return pi_o_4+tiny;
00097 case 1: return -pi_o_4-tiny;
00098 case 2: return 3.0*pi_o_4+tiny;
00099 case 3: return -3.0*pi_o_4-tiny;
00100 }
00101 } else {
00102 switch(m) {
00103 case 0: return zero ;
00104 case 1: return -zero ;
00105 case 2: return pi+tiny ;
00106 case 3: return -pi-tiny ;
00107 }
00108 }
00109 }
00110
00111 if(iy==0x7ff00000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
00112
00113
00114 k = (iy-ix)>>20;
00115 if(k > 60) z=pi_o_2+0.5*pi_lo;
00116 else if(hx<0&&k<-60) z=0.0;
00117 else z=atan(fabs(y/x));
00118 switch (m) {
00119 case 0: return z ;
00120 case 1: {
00121 u_int32_t zh;
00122 GET_HIGH_WORD(zh,z);
00123 SET_HIGH_WORD(z,zh ^ 0x80000000);
00124 }
00125 return z ;
00126 case 2: return pi-(z-pi_lo);
00127 default:
00128 return (z-pi_lo)-pi;
00129 }
00130 }