OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/src/e_atan2.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_atan2.c */ |
2 /* | 2 /* |
3 * ==================================================== | 3 * ==================================================== |
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
5 * | 5 * |
6 * Developed at SunSoft, a Sun Microsystems, Inc. business. | 6 * Developed at SunSoft, a Sun Microsystems, Inc. business. |
7 * Permission to use, copy, modify, and distribute this | 7 * Permission to use, copy, modify, and distribute this |
8 * software is freely granted, provided that this notice | 8 * software is freely granted, provided that this notice |
9 * is preserved. | 9 * is preserved. |
10 * ==================================================== | 10 * ==================================================== |
(...skipping 21 matching lines...) Expand all Loading... |
32 * | 32 * |
33 * Constants: | 33 * Constants: |
34 * The hexadecimal values are the intended ones for the following | 34 * The hexadecimal values are the intended ones for the following |
35 * constants. The decimal values may be used, provided that the | 35 * constants. The decimal values may be used, provided that the |
36 * compiler will convert from decimal to binary accurately enough | 36 * compiler will convert from decimal to binary accurately enough |
37 * to produce the hexadecimal values shown. | 37 * to produce the hexadecimal values shown. |
38 */ | 38 */ |
39 | 39 |
40 #include "libm.h" | 40 #include "libm.h" |
41 | 41 |
42 static const double | 42 static const double pi = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */ |
43 pi = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */ | 43 pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */ |
44 pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */ | |
45 | 44 |
46 double atan2(double y, double x) | 45 double atan2(double y, double x) { |
47 { | 46 double z; |
48 » double z; | 47 uint32_t m, lx, ly, ix, iy; |
49 » uint32_t m,lx,ly,ix,iy; | |
50 | 48 |
51 » if (isnan(x) || isnan(y)) | 49 if (isnan(x) || isnan(y)) |
52 » » return x+y; | 50 return x + y; |
53 » EXTRACT_WORDS(ix, lx, x); | 51 EXTRACT_WORDS(ix, lx, x); |
54 » EXTRACT_WORDS(iy, ly, y); | 52 EXTRACT_WORDS(iy, ly, y); |
55 » if ((ix-0x3ff00000 | lx) == 0) /* x = 1.0 */ | 53 if ((ix - 0x3ff00000 | lx) == 0) /* x = 1.0 */ |
56 » » return atan(y); | 54 return atan(y); |
57 » m = ((iy>>31)&1) | ((ix>>30)&2); /* 2*sign(x)+sign(y) */ | 55 m = ((iy >> 31) & 1) | ((ix >> 30) & 2); /* 2*sign(x)+sign(y) */ |
58 » ix = ix & 0x7fffffff; | 56 ix = ix & 0x7fffffff; |
59 » iy = iy & 0x7fffffff; | 57 iy = iy & 0x7fffffff; |
60 | 58 |
61 » /* when y = 0 */ | 59 /* when y = 0 */ |
62 » if ((iy|ly) == 0) { | 60 if ((iy | ly) == 0) { |
63 » » switch(m) { | 61 switch (m) { |
64 » » case 0: | 62 case 0: |
65 » » case 1: return y; /* atan(+-0,+anything)=+-0 */ | 63 case 1: |
66 » » case 2: return pi; /* atan(+0,-anything) = pi */ | 64 return y; /* atan(+-0,+anything)=+-0 */ |
67 » » case 3: return -pi; /* atan(-0,-anything) =-pi */ | 65 case 2: |
68 » » } | 66 return pi; /* atan(+0,-anything) = pi */ |
69 » } | 67 case 3: |
70 » /* when x = 0 */ | 68 return -pi; /* atan(-0,-anything) =-pi */ |
71 » if ((ix|lx) == 0) | 69 } |
72 » » return m&1 ? -pi/2 : pi/2; | 70 } |
73 » /* when x is INF */ | 71 /* when x = 0 */ |
74 » if (ix == 0x7ff00000) { | 72 if ((ix | lx) == 0) |
75 » » if (iy == 0x7ff00000) { | 73 return m & 1 ? -pi / 2 : pi / 2; |
76 » » » switch(m) { | 74 /* when x is INF */ |
77 » » » case 0: return pi/4; /* atan(+INF,+INF) */ | 75 if (ix == 0x7ff00000) { |
78 » » » case 1: return -pi/4; /* atan(-INF,+INF) */ | 76 if (iy == 0x7ff00000) { |
79 » » » case 2: return 3*pi/4; /* atan(+INF,-INF) */ | 77 switch (m) { |
80 » » » case 3: return -3*pi/4; /* atan(-INF,-INF) */ | 78 case 0: |
81 » » » } | 79 return pi / 4; /* atan(+INF,+INF) */ |
82 » » } else { | 80 case 1: |
83 » » » switch(m) { | 81 return -pi / 4; /* atan(-INF,+INF) */ |
84 » » » case 0: return 0.0; /* atan(+...,+INF) */ | 82 case 2: |
85 » » » case 1: return -0.0; /* atan(-...,+INF) */ | 83 return 3 * pi / 4; /* atan(+INF,-INF) */ |
86 » » » case 2: return pi; /* atan(+...,-INF) */ | 84 case 3: |
87 » » » case 3: return -pi; /* atan(-...,-INF) */ | 85 return -3 * pi / 4; /* atan(-INF,-INF) */ |
88 » » » } | 86 } |
89 » » } | 87 } else { |
90 » } | 88 switch (m) { |
91 » /* |y/x| > 0x1p64 */ | 89 case 0: |
92 » if (ix+(64<<20) < iy || iy == 0x7ff00000) | 90 return 0.0; /* atan(+...,+INF) */ |
93 » » return m&1 ? -pi/2 : pi/2; | 91 case 1: |
| 92 return -0.0; /* atan(-...,+INF) */ |
| 93 case 2: |
| 94 return pi; /* atan(+...,-INF) */ |
| 95 case 3: |
| 96 return -pi; /* atan(-...,-INF) */ |
| 97 } |
| 98 } |
| 99 } |
| 100 /* |y/x| > 0x1p64 */ |
| 101 if (ix + (64 << 20) < iy || iy == 0x7ff00000) |
| 102 return m & 1 ? -pi / 2 : pi / 2; |
94 | 103 |
95 » /* z = atan(|y/x|) without spurious underflow */ | 104 /* z = atan(|y/x|) without spurious underflow */ |
96 » if ((m&2) && iy+(64<<20) < ix) /* |y/x| < 0x1p-64, x<0 */ | 105 if ((m & 2) && iy + (64 << 20) < ix) /* |y/x| < 0x1p-64, x<0 */ |
97 » » z = 0; | 106 z = 0; |
98 » else | 107 else |
99 » » z = atan(fabs(y/x)); | 108 z = atan(fabs(y / x)); |
100 » switch (m) { | 109 switch (m) { |
101 » case 0: return z; /* atan(+,+) */ | 110 case 0: |
102 » case 1: return -z; /* atan(-,+) */ | 111 return z; /* atan(+,+) */ |
103 » case 2: return pi - (z-pi_lo); /* atan(+,-) */ | 112 case 1: |
104 » default: /* case 3 */ | 113 return -z; /* atan(-,+) */ |
105 » » return (z-pi_lo) - pi; /* atan(-,-) */ | 114 case 2: |
106 » } | 115 return pi - (z - pi_lo); /* atan(+,-) */ |
| 116 default: /* case 3 */ |
| 117 return (z - pi_lo) - pi; /* atan(-,-) */ |
| 118 } |
107 } | 119 } |
OLD | NEW |