1
2 /* @(#)e_atan2.c 5.1 93/09/24 */
3 /*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
7 * Developed at SunPro, a Sun Microsystems, Inc. business.
8 * Permission to use, copy, modify, and distribute this
9 * software is freely granted, provided that this notice
10 * is preserved.
11 * ====================================================
12 *
13 */
14
15 /* atan2(y,x)
16 * Method :
17 * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
18 * 2. Reduce x to positive by (if x and y are unexceptional):
19 * ARG (x+iy) = arctan(y/x) ... if x > 0,
20 * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0,
21 *
22 * Special cases:
23 *
24 * ATAN2((anything), NaN ) is NaN;
25 * ATAN2(NAN , (anything) ) is NaN;
26 * ATAN2(+-0, +(anything but NaN)) is +-0 ;
27 * ATAN2(+-0, -(anything but NaN)) is +-pi ;
28 * ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
29 * ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
30 * ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
31 * ATAN2(+-INF,+INF ) is +-pi/4 ;
32 * ATAN2(+-INF,-INF ) is +-3pi/4;
33 * ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
34 *
35 * Constants:
36 * The hexadecimal values are the intended ones for the following
37 * constants. The decimal values may be used, provided that the
38 * compiler will convert from decimal to binary accurately enough
39 * to produce the hexadecimal values shown.
40 */
41
42 #include "fdlibm.h"
43
44 #ifdef _NEED_FLOAT64
45
46 static const __float64
47 tiny = _F_64(1.0e-300), zero = _F_64(0.0),
48 pi_o_4 = _F_64(7.8539816339744827900E-01), /* 0x3FE921FB, 0x54442D18 */
49 pi_o_2 = _F_64(1.5707963267948965580E+00), /* 0x3FF921FB, 0x54442D18 */
50 pi = _F_64(3.1415926535897931160E+00), /* 0x400921FB, 0x54442D18 */
51 pi_lo = _F_64(1.2246467991473531772E-16); /* 0x3CA1A626, 0x33145C07 */
52
53 __float64
atan264(__float64 y,__float64 x)54 atan264(__float64 y, __float64 x)
55 {
56 __float64 z;
57 __int32_t k, m, hx, hy, ix, iy;
58 __uint32_t lx, ly;
59
60 EXTRACT_WORDS(hx, lx, x);
61 ix = hx & 0x7fffffff;
62 EXTRACT_WORDS(hy, ly, y);
63 iy = hy & 0x7fffffff;
64 if (((ix | ((lx | -lx) >> 31)) > 0x7ff00000) ||
65 ((iy | ((ly | -ly) >> 31)) > 0x7ff00000)) /* x or y is NaN */
66 return x + y;
67 if (((hx - 0x3ff00000) | lx) == 0)
68 return atan64(y); /* x=1.0 */
69 m = ((hy >> 31) & 1) | ((hx >> 30) & 2); /* 2*sign(x)+sign(y) */
70
71 /* when y = 0 */
72 if ((iy | ly) == 0) {
73 switch (m) {
74 case 0:
75 case 1:
76 return y; /* atan(+-0,+anything)=+-0 */
77 case 2:
78 return pi + tiny; /* atan(+0,-anything) = pi */
79 case 3:
80 return -pi - tiny; /* atan(-0,-anything) =-pi */
81 }
82 }
83 /* when x = 0 */
84 if ((ix | lx) == 0)
85 return (hy < 0) ? -pi_o_2 - tiny : pi_o_2 + tiny;
86
87 /* when x is INF */
88 if (ix == 0x7ff00000) {
89 if (iy == 0x7ff00000) {
90 switch (m) {
91 case 0:
92 return pi_o_4 + tiny; /* atan(+INF,+INF) */
93 case 1:
94 return -pi_o_4 - tiny; /* atan(-INF,+INF) */
95 case 2:
96 return _F_64(3.0) * pi_o_4 + tiny; /*atan(+INF,-INF)*/
97 case 3:
98 return _F_64(-3.0) * pi_o_4 - tiny; /*atan(-INF,-INF)*/
99 }
100 } else {
101 switch (m) {
102 case 0:
103 return zero; /* atan(+...,+INF) */
104 case 1:
105 return -zero; /* atan(-...,+INF) */
106 case 2:
107 return pi + tiny; /* atan(+...,-INF) */
108 case 3:
109 return -pi - tiny; /* atan(-...,-INF) */
110 }
111 }
112 }
113 /* when y is INF */
114 if (iy == 0x7ff00000)
115 return (hy < 0) ? -pi_o_2 - tiny : pi_o_2 + tiny;
116
117 /* compute y/x */
118 k = (iy - ix) >> 20;
119 if (k > 60)
120 z = pi_o_2 + _F_64(0.5) * pi_lo; /* |y/x| > 2**60 */
121 else if (hx < 0 && k < -60)
122 z = _F_64(0.0); /* |y|/x < -2**60 */
123 else
124 z = atan64(check_uflow(fabs64(y / x))); /* safe to do y/x */
125 switch (m) {
126 case 0:
127 return z; /* atan(+,+) */
128 case 1: {
129 __uint32_t zh;
130 GET_HIGH_WORD(zh, z);
131 SET_HIGH_WORD(z, zh ^ 0x80000000);
132 }
133 return z; /* atan(-,+) */
134 case 2:
135 return pi - (z - pi_lo); /* atan(+,-) */
136 default: /* case 3 */
137 return (z - pi_lo) - pi; /* atan(-,-) */
138 }
139 }
140
141 _MATH_ALIAS_d_dd(atan2)
142
143 #endif /* _NEED_FLOAT64 */
144
