1 /* origin: FreeBSD /usr/src/lib/msun/src/e_atan2f.c */ 2 /* 3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 4 */ 5 /* 6 * ==================================================== 7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 8 * 9 * Developed at SunPro, a Sun Microsystems, Inc. business. 10 * Permission to use, copy, modify, and distribute this 11 * software is freely granted, provided that this notice 12 * is preserved. 13 * ==================================================== 14 */ 15 16 #include "libm.h" 17 18 static const float 19 pi = 3.1415927410e+00, /* 0x40490fdb */ 20 pi_lo = -8.7422776573e-08; /* 0xb3bbbd2e */ 21 22 float atan2f(float y, float x) 23 { 24 float z; 25 uint32_t m,ix,iy; 26 27 if (isnan(x) || isnan(y)) 28 return x+y; 29 GET_FLOAT_WORD(ix, x); 30 GET_FLOAT_WORD(iy, y); 31 if (ix == 0x3f800000) /* x=1.0 */ 32 return atanf(y); 33 m = ((iy>>31)&1) | ((ix>>30)&2); /* 2*sign(x)+sign(y) */ 34 ix &= 0x7fffffff; 35 iy &= 0x7fffffff; 36 37 /* when y = 0 */ 38 if (iy == 0) { 39 switch (m) { 40 case 0: 41 case 1: return y; /* atan(+-0,+anything)=+-0 */ 42 case 2: return pi; /* atan(+0,-anything) = pi */ 43 case 3: return -pi; /* atan(-0,-anything) =-pi */ 44 } 45 } 46 /* when x = 0 */ 47 if (ix == 0) 48 return m&1 ? -pi/2 : pi/2; 49 /* when x is INF */ 50 if (ix == 0x7f800000) { 51 if (iy == 0x7f800000) { 52 switch (m) { 53 case 0: return pi/4; /* atan(+INF,+INF) */ 54 case 1: return -pi/4; /* atan(-INF,+INF) */ 55 case 2: return 3*pi/4; /*atan(+INF,-INF)*/ 56 case 3: return -3*pi/4; /*atan(-INF,-INF)*/ 57 } 58 } else { 59 switch (m) { 60 case 0: return 0.0f; /* atan(+...,+INF) */ 61 case 1: return -0.0f; /* atan(-...,+INF) */ 62 case 2: return pi; /* atan(+...,-INF) */ 63 case 3: return -pi; /* atan(-...,-INF) */ 64 } 65 } 66 } 67 /* |y/x| > 0x1p26 */ 68 if (ix+(26<<23) < iy || iy == 0x7f800000) 69 return m&1 ? -pi/2 : pi/2; 70 71 /* z = atan(|y/x|) with correct underflow */ 72 if ((m&2) && iy+(26<<23) < ix) /*|y/x| < 0x1p-26, x < 0 */ 73 z = 0.0; 74 else 75 z = atanf(fabsf(y/x)); 76 switch (m) { 77 case 0: return z; /* atan(+,+) */ 78 case 1: return -z; /* atan(-,+) */ 79 case 2: return pi - (z-pi_lo); /* atan(+,-) */ 80 default: /* case 3 */ 81 return (z-pi_lo) - pi; /* atan(-,-) */ 82 } 83 } 84