/********************************************************************** * Copyright (c) Digital Equipment Corporation 1984, 1985, 1986. * * All Rights Reserved. * * Reference "/usr/src/COPYRIGHT" for applicable restrictions. * **********************************************************************/ /* SCCSID: @(#)atan.c 3.0 4/22/86 */ /* (System 5) atan.c 1.10 */ /*LINTLIBRARY*/ /* * atan returns the arctangent of its double-precision argument, * in the range [-pi/2, pi/2]. * There are no error returns. * * atan2(y, x) returns the arctangent of y/x, * in the range [-pi, pi]. * atan2 discovers what quadrant the angle is in and calls atan. * atan2 returns EDOM error and value 0 if both arguments are zero. * * Coefficients are #5077 from Hart & Cheney (19.56D). */ #include <math.h> #include <values.h> #include <errno.h> #define SQ2_M_1 0.41421356237309504880 #define SQ2_P_1 2.41421356237309504880 #define NEG 1 #define INV 2 #define MID 4 double atan(x) register double x; { register int type = 0; if (x < 0) { x = -x; type = NEG; } if (x > SQ2_P_1) { x = 1/x; type |= INV; } else if (x > SQ2_M_1) { x = (x - 1)/(x + 1); type |= MID; } if (x < -X_EPS || x > X_EPS) { /* skip for efficiency and to prevent underflow */ static double p[] = { 0.161536412982230228262e2, 0.26842548195503973794141e3, 0.11530293515404850115428136e4, 0.178040631643319697105464587e4, 0.89678597403663861959987488e3, }, q[] = { 1.0, 0.5895697050844462222791e2, 0.536265374031215315104235e3, 0.16667838148816337184521798e4, 0.207933497444540981287275926e4, 0.89678597403663861962481162e3, }; register double xsq = x * x; x *= _POLY4(xsq, p)/_POLY5(xsq, q); } if (type & INV) x = M_PI_2 - x; else if (type & MID) x += M_PI_4; return (type & NEG ? -x : x); } double atan2(y, x) register double y, x; { register int neg_y = 0; struct exception exc; exc.arg1 = y; if (!x && !y) { exc.type = DOMAIN; exc.name = "atan2"; exc.arg2 = x; exc.retval = 0.0; if (!matherr(&exc)) { (void) write(2, "atan2: DOMAIN error\n", 20); errno = EDOM; } return (exc.retval); } /* * The next lines determine if |x| is negligible compared to |y|, * without dividing, and without adding values of the same sign. */ if (exc.arg1 < 0) { exc.arg1 = -exc.arg1; neg_y++; } if (exc.arg1 - _ABS(x) == exc.arg1) return (neg_y ? -M_PI_2 : M_PI_2); /* * The next line assumes that if y/x underflows the result * is zero with no error indication, so it's safe to divide. */ return ((y = atan(y/x)) == 0 || x > 0 ? y : neg_y ? y - M_PI : y + M_PI); }