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