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 }