| @@ -1,123 +1,123 @@ | | | @@ -1,123 +1,123 @@ |
1 | /* @(#)e_atan2.c 5.1 93/09/24 */ | | 1 | /* @(#)e_atan2.c 5.1 93/09/24 */ |
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 SunPro, a Sun Microsystems, Inc. business. | | 6 | * Developed at SunPro, 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 | #include <sys/cdefs.h> | | 13 | #include <sys/cdefs.h> |
14 | #if defined(LIBM_SCCS) && !defined(lint) | | 14 | #if defined(LIBM_SCCS) && !defined(lint) |
15 | __RCSID("$NetBSD: e_atan2.c,v 1.12 2002/05/26 22:01:48 wiz Exp $"); | | 15 | __RCSID("$NetBSD: e_atan2.c,v 1.13 2018/03/10 09:44:47 eadler Exp $"); |
16 | #endif | | 16 | #endif |
17 | | | 17 | |
18 | /* __ieee754_atan2(y,x) | | 18 | /* __ieee754_atan2(y,x) |
19 | * Method : | | 19 | * Method : |
20 | * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). | | 20 | * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). |
21 | * 2. Reduce x to positive by (if x and y are unexceptional): | | 21 | * 2. Reduce x to positive by (if x and y are unexceptional): |
22 | * ARG (x+iy) = arctan(y/x) ... if x > 0, | | 22 | * ARG (x+iy) = arctan(y/x) ... if x > 0, |
23 | * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, | | 23 | * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, |
24 | * | | 24 | * |
25 | * Special cases: | | 25 | * Special cases: |
26 | * | | 26 | * |
27 | * ATAN2((anything), NaN ) is NaN; | | 27 | * ATAN2((anything), NaN ) is NaN; |
28 | * ATAN2(NAN , (anything) ) is NaN; | | 28 | * ATAN2(NAN , (anything) ) is NaN; |
29 | * ATAN2(+-0, +(anything but NaN)) is +-0 ; | | 29 | * ATAN2(+-0, +(anything but NaN)) is +-0 ; |
30 | * ATAN2(+-0, -(anything but NaN)) is +-pi ; | | 30 | * ATAN2(+-0, -(anything but NaN)) is +-pi ; |
31 | * ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2; | | 31 | * ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2; |
32 | * ATAN2(+-(anything but INF and NaN), +INF) is +-0 ; | | 32 | * ATAN2(+-(anything but INF and NaN), +INF) is +-0 ; |
33 | * ATAN2(+-(anything but INF and NaN), -INF) is +-pi; | | 33 | * ATAN2(+-(anything but INF and NaN), -INF) is +-pi; |
34 | * ATAN2(+-INF,+INF ) is +-pi/4 ; | | 34 | * ATAN2(+-INF,+INF ) is +-pi/4 ; |
35 | * ATAN2(+-INF,-INF ) is +-3pi/4; | | 35 | * ATAN2(+-INF,-INF ) is +-3pi/4; |
36 | * ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2; | | 36 | * ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2; |
37 | * | | 37 | * |
38 | * Constants: | | 38 | * Constants: |
39 | * The hexadecimal values are the intended ones for the following | | 39 | * The hexadecimal values are the intended ones for the following |
40 | * constants. The decimal values may be used, provided that the | | 40 | * constants. The decimal values may be used, provided that the |
41 | * compiler will convert from decimal to binary accurately enough | | 41 | * compiler will convert from decimal to binary accurately enough |
42 | * to produce the hexadecimal values shown. | | 42 | * to produce the hexadecimal values shown. |
43 | */ | | 43 | */ |
44 | | | 44 | |
45 | #include "math.h" | | 45 | #include "math.h" |
46 | #include "math_private.h" | | 46 | #include "math_private.h" |
47 | | | 47 | |
48 | static const double | | 48 | static const double |
49 | tiny = 1.0e-300, | | 49 | tiny = 1.0e-300, |
50 | zero = 0.0, | | 50 | zero = 0.0, |
51 | pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */ | | 51 | pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */ |
52 | pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */ | | 52 | pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */ |
53 | pi = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */ | | 53 | pi = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */ |
54 | pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */ | | 54 | pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */ |
55 | | | 55 | |
56 | double | | 56 | double |
57 | __ieee754_atan2(double y, double x) | | 57 | __ieee754_atan2(double y, double x) |
58 | { | | 58 | { |
59 | double z; | | 59 | double z; |
60 | int32_t k,m,hx,hy,ix,iy; | | 60 | int32_t k,m,hx,hy,ix,iy; |
61 | u_int32_t lx,ly; | | 61 | u_int32_t lx,ly; |
62 | | | 62 | |
63 | EXTRACT_WORDS(hx,lx,x); | | 63 | EXTRACT_WORDS(hx,lx,x); |
64 | ix = hx&0x7fffffff; | | 64 | ix = hx&0x7fffffff; |
65 | EXTRACT_WORDS(hy,ly,y); | | 65 | EXTRACT_WORDS(hy,ly,y); |
66 | iy = hy&0x7fffffff; | | 66 | iy = hy&0x7fffffff; |
67 | if(((ix|((lx|-lx)>>31))>0x7ff00000)|| | | 67 | if(((ix|((lx|-lx)>>31))>0x7ff00000)|| |
68 | ((iy|((ly|-ly)>>31))>0x7ff00000)) /* x or y is NaN */ | | 68 | ((iy|((ly|-ly)>>31))>0x7ff00000)) /* x or y is NaN */ |
69 | return x+y; | | 69 | return x+y; |
70 | if(((hx-0x3ff00000)|lx)==0) return atan(y); /* x=1.0 */ | | 70 | if(hx==0x3ff00000&&lx==0) return atan(y); /* x=1.0 */ |
71 | m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */ | | 71 | m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */ |
72 | | | 72 | |
73 | /* when y = 0 */ | | 73 | /* when y = 0 */ |
74 | if((iy|ly)==0) { | | 74 | if((iy|ly)==0) { |
75 | switch(m) { | | 75 | switch(m) { |
76 | case 0: | | 76 | case 0: |
77 | case 1: return y; /* atan(+-0,+anything)=+-0 */ | | 77 | case 1: return y; /* atan(+-0,+anything)=+-0 */ |
78 | case 2: return pi+tiny;/* atan(+0,-anything) = pi */ | | 78 | case 2: return pi+tiny;/* atan(+0,-anything) = pi */ |
79 | case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */ | | 79 | case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */ |
80 | } | | 80 | } |
81 | } | | 81 | } |
82 | /* when x = 0 */ | | 82 | /* when x = 0 */ |
83 | if((ix|lx)==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny; | | 83 | if((ix|lx)==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny; |
84 | | | 84 | |
85 | /* when x is INF */ | | 85 | /* when x is INF */ |
86 | if(ix==0x7ff00000) { | | 86 | if(ix==0x7ff00000) { |
87 | if(iy==0x7ff00000) { | | 87 | if(iy==0x7ff00000) { |
88 | switch(m) { | | 88 | switch(m) { |
89 | case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */ | | 89 | case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */ |
90 | case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */ | | 90 | case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */ |
91 | case 2: return 3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/ | | 91 | case 2: return 3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/ |
92 | case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/ | | 92 | case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/ |
93 | } | | 93 | } |
94 | } else { | | 94 | } else { |
95 | switch(m) { | | 95 | switch(m) { |
96 | case 0: return zero ; /* atan(+...,+INF) */ | | 96 | case 0: return zero ; /* atan(+...,+INF) */ |
97 | case 1: return -zero ; /* atan(-...,+INF) */ | | 97 | case 1: return -zero ; /* atan(-...,+INF) */ |
98 | case 2: return pi+tiny ; /* atan(+...,-INF) */ | | 98 | case 2: return pi+tiny ; /* atan(+...,-INF) */ |
99 | case 3: return -pi-tiny ; /* atan(-...,-INF) */ | | 99 | case 3: return -pi-tiny ; /* atan(-...,-INF) */ |
100 | } | | 100 | } |
101 | } | | 101 | } |
102 | } | | 102 | } |
103 | /* when y is INF */ | | 103 | /* when y is INF */ |
104 | if(iy==0x7ff00000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny; | | 104 | if(iy==0x7ff00000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny; |
105 | | | 105 | |
106 | /* compute y/x */ | | 106 | /* compute y/x */ |
107 | k = (iy-ix)>>20; | | 107 | k = (iy-ix)>>20; |
108 | if(k > 60) z=pi_o_2+0.5*pi_lo; /* |y/x| > 2**60 */ | | 108 | if(k > 60) z=pi_o_2+0.5*pi_lo; /* |y/x| > 2**60 */ |
109 | else if(hx<0&&k<-60) z=0.0; /* |y|/x < -2**60 */ | | 109 | else if(hx<0&&k<-60) z=0.0; /* |y|/x < -2**60 */ |
110 | else z=atan(fabs(y/x)); /* safe to do y/x */ | | 110 | else z=atan(fabs(y/x)); /* safe to do y/x */ |
111 | switch (m) { | | 111 | switch (m) { |
112 | case 0: return z ; /* atan(+,+) */ | | 112 | case 0: return z ; /* atan(+,+) */ |
113 | case 1: { | | 113 | case 1: { |
114 | u_int32_t zh; | | 114 | u_int32_t zh; |
115 | GET_HIGH_WORD(zh,z); | | 115 | GET_HIGH_WORD(zh,z); |
116 | SET_HIGH_WORD(z,zh ^ 0x80000000); | | 116 | SET_HIGH_WORD(z,zh ^ 0x80000000); |
117 | } | | 117 | } |
118 | return z ; /* atan(-,+) */ | | 118 | return z ; /* atan(-,+) */ |
119 | case 2: return pi-(z-pi_lo);/* atan(+,-) */ | | 119 | case 2: return pi-(z-pi_lo);/* atan(+,-) */ |
120 | default: /* case 3 */ | | 120 | default: /* case 3 */ |
121 | return (z-pi_lo)-pi;/* atan(-,-) */ | | 121 | return (z-pi_lo)-pi;/* atan(-,-) */ |
122 | } | | 122 | } |
123 | } | | 123 | } |