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