Commit | Line | Data |
---|---|---|
7fde0e2e KB |
1 | /*- |
2 | * Copyright (c) 1990 The Regents of the University of California. | |
3 | * All rights reserved. | |
4 | * | |
5 | * This code is derived from software contributed to Berkeley by | |
6 | * the Systems Programming Group of the University of Utah Computer | |
7 | * Science Department. | |
8 | * | |
9 | * %sccs.include.redist.c% | |
10 | */ | |
11 | ||
12 | #ifndef lint | |
13 | static char sccsid[] = "@(#)atan2.c 5.1 (Berkeley) %G%"; | |
14 | #endif /* not lint */ | |
15 | ||
16 | /* | |
17 | * ATAN2(Y,X) | |
18 | * RETURN ARG (X+iY) | |
19 | * DOUBLE PRECISION (IEEE DOUBLE 53 BITS) | |
20 | * | |
21 | * Scaled down version to weed out special cases. "Normal" cases are | |
22 | * handled by calling atan2__A(), an assembly coded support routine in | |
23 | * support.s. | |
24 | * | |
25 | * Required system supported functions : | |
26 | * copysign(x,y) | |
27 | * atan2__A(y,x) | |
28 | * | |
29 | * Method : | |
30 | * 1. Deal with special cases | |
31 | * 2. Call atan2__A() to do the others | |
32 | * | |
33 | * Special cases: | |
34 | * Notations: atan2(y,x) == ARG (x+iy) == ARG(x,y). | |
35 | * | |
36 | * ARG( NAN , (anything) ) is NaN; | |
37 | * ARG( (anything), NaN ) is NaN; | |
38 | * ARG(+(anything but NaN), +-0) is +-0 ; | |
39 | * ARG(-(anything but NaN), +-0) is +-PI ; | |
40 | * ARG( 0, +-(anything but 0 and NaN) ) is +-PI/2; | |
41 | * ARG( +INF,+-(anything but INF and NaN) ) is +-0 ; | |
42 | * ARG( -INF,+-(anything but INF and NaN) ) is +-PI; | |
43 | * ARG( +INF,+-INF ) is +-PI/4 ; | |
44 | * ARG( -INF,+-INF ) is +-3PI/4; | |
45 | * ARG( (anything but,0,NaN, and INF),+-INF ) is +-PI/2; | |
46 | * | |
47 | * Accuracy: | |
48 | * atan2(y,x) returns (PI/pi) * the exact ARG (x+iy) nearly rounded, | |
49 | * where | |
50 | * | |
51 | * in decimal: | |
52 | * pi = 3.141592653589793 23846264338327 ..... | |
53 | * 53 bits PI = 3.141592653589793 115997963 ..... , | |
54 | * 56 bits PI = 3.141592653589793 227020265 ..... , | |
55 | * | |
56 | * in hexadecimal: | |
57 | * pi = 3.243F6A8885A308D313198A2E.... | |
58 | * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps | |
59 | * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps | |
60 | * | |
61 | * In a test run with 356,000 random argument on [-1,1] * [-1,1] on a | |
62 | * VAX, the maximum observed error was 1.41 ulps (units of the last place) | |
63 | * compared with (PI/pi)*(the exact ARG(x+iy)). | |
64 | * | |
65 | * Note: | |
66 | * We use machine PI (the true pi rounded) in place of the actual | |
67 | * value of pi for all the trig and inverse trig functions. In general, | |
68 | * if trig is one of sin, cos, tan, then computed trig(y) returns the | |
69 | * exact trig(y*pi/PI) nearly rounded; correspondingly, computed arctrig | |
70 | * returns the exact arctrig(y)*PI/pi nearly rounded. These guarantee the | |
71 | * trig functions have period PI, and trig(arctrig(x)) returns x for | |
72 | * all critical values x. | |
73 | * | |
74 | * Constants: | |
75 | * The hexadecimal values are the intended ones for the following constants. | |
76 | * The decimal values may be used, provided that the compiler will convert | |
77 | * from decimal to binary accurately enough to produce the hexadecimal values | |
78 | * shown. | |
79 | */ | |
80 | ||
81 | static double | |
82 | PIo4 = 7.8539816339744827900E-1 , /*Hex 2^ -1 * 1.921FB54442D18 */ | |
83 | PIo2 = 1.5707963267948965580E0 , /*Hex 2^ 0 * 1.921FB54442D18 */ | |
84 | PI = 3.1415926535897931160E0 ; /*Hex 2^ 1 * 1.921FB54442D18 */ | |
85 | ||
86 | double atan2(y,x) | |
87 | double y,x; | |
88 | { | |
89 | static double zero=0, one=1; | |
90 | double copysign(),atan2__A(),signy,signx; | |
91 | int finite(); | |
92 | ||
93 | /* if x or y is NAN */ | |
94 | if(x!=x) return(x); if(y!=y) return(y); | |
95 | ||
96 | /* copy down the sign of y and x */ | |
97 | signy = copysign(one,y); | |
98 | signx = copysign(one,x); | |
99 | ||
100 | /* when y = 0 */ | |
101 | if(y==zero) return((signx==one)?y:copysign(PI,signy)); | |
102 | ||
103 | /* when x = 0 */ | |
104 | if(x==zero) return(copysign(PIo2,signy)); | |
105 | ||
106 | /* when x is INF */ | |
107 | if(!finite(x)) | |
108 | if(!finite(y)) | |
109 | return(copysign((signx==one)?PIo4:3*PIo4,signy)); | |
110 | else | |
111 | return(copysign((signx==one)?zero:PI,signy)); | |
112 | ||
113 | /* when y is INF */ | |
114 | if(!finite(y)) return(copysign(PIo2,signy)); | |
115 | ||
116 | /* else let atan2__A do the work */ | |
117 | return(atan2__A(y,x)); | |
118 | } |