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 40 xatan(double arg) 41 { 42 double argsq, value; 43 44 argsq = arg*arg; 45 value = ((((p4*argsq + p3)*argsq + p2)*argsq + p1)*argsq + p0); 46 value = value/(((((argsq + q4)*argsq + q3)*argsq + q2)*argsq + q1)*argsq + q0); 47 return value*arg; 48 } 49 50 /* 51 satan reduces its argument (known to be positive) 52 to the range [0,0.414...] and calls xatan. 53 */ 54 55 static 56 double 57 satan(double arg) 58 { 59 60 if(arg < sq2m1) 61 return xatan(arg); 62 if(arg > sq2p1) 63 return pio2 - xatan(1.0/arg); 64 return pio4 + xatan((arg-1.0)/(arg+1.0)); 65 } 66 67 /* 68 atan makes its argument positive and 69 calls the inner routine satan. 70 */ 71 72 double 73 atan(double arg) 74 { 75 76 if(arg > 0) 77 return satan(arg); 78 return -satan(-arg); 79 } 80