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 <u.h> 16 #include <libc.h> 17 18 #define sq2p1 2.414213562373095048802e0 19 #define sq2m1 .414213562373095048802e0 20 #define p4 .161536412982230228262e2 21 #define p3 .26842548195503973794141e3 22 #define p2 .11530293515404850115428136e4 23 #define p1 .178040631643319697105464587e4 24 #define p0 .89678597403663861959987488e3 25 #define q4 .5895697050844462222791e2 26 #define q3 .536265374031215315104235e3 27 #define q2 .16667838148816337184521798e4 28 #define q1 .207933497444540981287275926e4 29 #define q0 .89678597403663861962481162e3 30 31 32 /* 33 xatan evaluates a series valid in the 34 range [-0.414...,+0.414...]. (tan(pi/8)) 35 */ 36 37 static 38 double xatan(double arg)39xatan(double arg) 40 { 41 double argsq, value; 42 43 argsq = arg*arg; 44 value = ((((p4*argsq + p3)*argsq + p2)*argsq + p1)*argsq + p0); 45 value = value/(((((argsq + q4)*argsq + q3)*argsq + q2)*argsq + q1)*argsq + q0); 46 return value*arg; 47 } 48 49 /* 50 satan reduces its argument (known to be positive) 51 to the range [0,0.414...] and calls xatan. 52 */ 53 54 static 55 double satan(double arg)56satan(double arg) 57 { 58 59 if(arg < sq2m1) 60 return xatan(arg); 61 if(arg > sq2p1) 62 return PIO2 - xatan(1/arg); 63 return PIO2/2 + xatan((arg-1)/(arg+1)); 64 } 65 66 /* 67 atan makes its argument positive and 68 calls the inner routine satan. 69 */ 70 71 double atan(double arg)72atan(double arg) 73 { 74 75 if(arg > 0) 76 return satan(arg); 77 return -satan(-arg); 78 } 79