1 /* 2 floating-point arctangent 3 4 atan returns the value of the arctangent of its 5 argument in the range [-pi/2,pi/2]. 6 7 atan2 returns the arctangent of arg1/arg2 8 in the range [-pi,pi]. 9 10 there are no error returns. 11 12 coefficients are #5077 from Hart & Cheney. (19.56D) 13 */ 14 15 #include <math.h> 16 17 #define sq2p1 2.414213562373095048802e0 18 #define sq2m1 .414213562373095048802e0 19 #define pio2 1.570796326794896619231e0 20 #define pio4 .785398163397448309615e0 21 #define p4 .161536412982230228262e2 22 #define p3 .26842548195503973794141e3 23 #define p2 .11530293515404850115428136e4 24 #define p1 .178040631643319697105464587e4 25 #define p0 .89678597403663861959987488e3 26 #define q4 .5895697050844462222791e2 27 #define q3 .536265374031215315104235e3 28 #define q2 .16667838148816337184521798e4 29 #define q1 .207933497444540981287275926e4 30 #define q0 .89678597403663861962481162e3 31 32 33 /* 34 xatan evaluates a series valid in the 35 range [-0.414...,+0.414...]. 36 */ 37 38 static 39 double xatan(double arg)40xatan(double arg) 41 { 42 double argsq, value; 43 44 /* get denormalized add in following if range arg**10 is much smaller 45 than q1, so check for that case 46 */ 47 if(-.01 < arg && arg < .01) 48 value = p0/q0; 49 else { 50 argsq = arg*arg; 51 value = ((((p4*argsq + p3)*argsq + p2)*argsq + p1)*argsq + p0); 52 value = value/(((((argsq + q4)*argsq + q3)*argsq + q2)*argsq + q1)*argsq + q0); 53 } 54 return value*arg; 55 } 56 57 /* 58 satan reduces its argument (known to be positive) 59 to the range [0,0.414...] and calls xatan. 60 */ 61 62 static 63 double satan(double arg)64satan(double arg) 65 { 66 67 if(arg < sq2m1) 68 return xatan(arg); 69 if(arg > sq2p1) 70 return pio2 - xatan(1.0/arg); 71 return pio4 + xatan((arg-1.0)/(arg+1.0)); 72 } 73 74 /* 75 atan makes its argument positive and 76 calls the inner routine satan. 77 */ 78 79 double atan(double arg)80atan(double arg) 81 { 82 83 if(arg > 0) 84 return satan(arg); 85 return -satan(-arg); 86 } 87