OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/src/e_atan2l.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_atan2l.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 * ==================================================== |
11 * | 11 * |
12 */ | 12 */ |
13 /* | 13 /* |
14 * See comments in atan2.c. | 14 * See comments in atan2.c. |
15 * Converted to long double by David Schultz <das@FreeBSD.ORG>. | 15 * Converted to long double by David Schultz <das@FreeBSD.ORG>. |
16 */ | 16 */ |
17 | 17 |
18 #include "libm.h" | 18 #include "libm.h" |
19 | 19 |
20 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 | 20 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 |
21 long double atan2l(long double y, long double x) | 21 long double atan2l(long double y, long double x) { |
22 { | 22 return atan2(y, x); |
23 » return atan2(y, x); | |
24 } | 23 } |
25 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 | 24 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 |
26 #include "__invtrigl.h" | 25 #include "__invtrigl.h" |
27 | 26 |
28 long double atan2l(long double y, long double x) | 27 long double atan2l(long double y, long double x) { |
29 { | 28 union ldshape ux, uy; |
30 » union ldshape ux, uy; | 29 long double z; |
31 » long double z; | 30 int m, ex, ey; |
32 » int m, ex, ey; | |
33 | 31 |
34 » if (isnan(x) || isnan(y)) | 32 if (isnan(x) || isnan(y)) |
35 » » return x+y; | 33 return x + y; |
36 » if (x == 1) | 34 if (x == 1) |
37 » » return atanl(y); | 35 return atanl(y); |
38 » ux.f = x; | 36 ux.f = x; |
39 » uy.f = y; | 37 uy.f = y; |
40 » ex = ux.i.se & 0x7fff; | 38 ex = ux.i.se & 0x7fff; |
41 » ey = uy.i.se & 0x7fff; | 39 ey = uy.i.se & 0x7fff; |
42 » m = 2*(ux.i.se>>15) | uy.i.se>>15; | 40 m = 2 * (ux.i.se >> 15) | uy.i.se >> 15; |
43 » if (y == 0) { | 41 if (y == 0) { |
44 » » switch(m) { | 42 switch (m) { |
45 » » case 0: | 43 case 0: |
46 » » case 1: return y; /* atan(+-0,+anything)=+-0 */ | 44 case 1: |
47 » » case 2: return 2*pio2_hi; /* atan(+0,-anything) = pi */ | 45 return y; /* atan(+-0,+anything)=+-0 */ |
48 » » case 3: return -2*pio2_hi; /* atan(-0,-anything) =-pi */ | 46 case 2: |
49 » » } | 47 return 2 * pio2_hi; /* atan(+0,-anything) = pi */ |
50 » } | 48 case 3: |
51 » if (x == 0) | 49 return -2 * pio2_hi; /* atan(-0,-anything) =-pi */ |
52 » » return m&1 ? -pio2_hi : pio2_hi; | 50 } |
53 » if (ex == 0x7fff) { | 51 } |
54 » » if (ey == 0x7fff) { | 52 if (x == 0) |
55 » » » switch(m) { | 53 return m & 1 ? -pio2_hi : pio2_hi; |
56 » » » case 0: return pio2_hi/2; /* atan(+INF,+INF) */ | 54 if (ex == 0x7fff) { |
57 » » » case 1: return -pio2_hi/2; /* atan(-INF,+INF) */ | 55 if (ey == 0x7fff) { |
58 » » » case 2: return 1.5*pio2_hi; /* atan(+INF,-INF) */ | 56 switch (m) { |
59 » » » case 3: return -1.5*pio2_hi; /* atan(-INF,-INF) */ | 57 case 0: |
60 » » » } | 58 return pio2_hi / 2; /* atan(+INF,+INF) */ |
61 » » } else { | 59 case 1: |
62 » » » switch(m) { | 60 return -pio2_hi / 2; /* atan(-INF,+INF) */ |
63 » » » case 0: return 0.0; /* atan(+...,+INF) */ | 61 case 2: |
64 » » » case 1: return -0.0; /* atan(-...,+INF) */ | 62 return 1.5 * pio2_hi; /* atan(+INF,-INF) */ |
65 » » » case 2: return 2*pio2_hi; /* atan(+...,-INF) */ | 63 case 3: |
66 » » » case 3: return -2*pio2_hi; /* atan(-...,-INF) */ | 64 return -1.5 * pio2_hi; /* atan(-INF,-INF) */ |
67 » » » } | 65 } |
68 » » } | 66 } else { |
69 » } | 67 switch (m) { |
70 » if (ex+120 < ey || ey == 0x7fff) | 68 case 0: |
71 » » return m&1 ? -pio2_hi : pio2_hi; | 69 return 0.0; /* atan(+...,+INF) */ |
72 » /* z = atan(|y/x|) without spurious underflow */ | 70 case 1: |
73 » if ((m&2) && ey+120 < ex) /* |y/x| < 0x1p-120, x<0 */ | 71 return -0.0; /* atan(-...,+INF) */ |
74 » » z = 0.0; | 72 case 2: |
75 » else | 73 return 2 * pio2_hi; /* atan(+...,-INF) */ |
76 » » z = atanl(fabsl(y/x)); | 74 case 3: |
77 » switch (m) { | 75 return -2 * pio2_hi; /* atan(-...,-INF) */ |
78 » case 0: return z; /* atan(+,+) */ | 76 } |
79 » case 1: return -z; /* atan(-,+) */ | 77 } |
80 » case 2: return 2*pio2_hi-(z-2*pio2_lo); /* atan(+,-) */ | 78 } |
81 » default: /* case 3 */ | 79 if (ex + 120 < ey || ey == 0x7fff) |
82 » » return (z-2*pio2_lo)-2*pio2_hi; /* atan(-,-) */ | 80 return m & 1 ? -pio2_hi : pio2_hi; |
83 » } | 81 /* z = atan(|y/x|) without spurious underflow */ |
| 82 if ((m & 2) && ey + 120 < ex) /* |y/x| < 0x1p-120, x<0 */ |
| 83 z = 0.0; |
| 84 else |
| 85 z = atanl(fabsl(y / x)); |
| 86 switch (m) { |
| 87 case 0: |
| 88 return z; /* atan(+,+) */ |
| 89 case 1: |
| 90 return -z; /* atan(-,+) */ |
| 91 case 2: |
| 92 return 2 * pio2_hi - (z - 2 * pio2_lo); /* atan(+,-) */ |
| 93 default: /* case 3 */ |
| 94 return (z - 2 * pio2_lo) - 2 * pio2_hi; /* atan(-,-) */ |
| 95 } |
84 } | 96 } |
85 #endif | 97 #endif |
OLD | NEW |