1 /*- 2 * Copyright (c) 1990, 1993 3 * The Regents of the University of California. 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 8.1 (Berkeley) 06/04/93"; 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 } 119