1*9941Ssam /* @(#)tan.c 4.1 12/25/82 */ 2*9941Ssam 3*9941Ssam /* 4*9941Ssam floating point tangent 5*9941Ssam 6*9941Ssam A series is used after range reduction. 7*9941Ssam Coefficients are #4285 from Hart & Cheney. (19.74D) 8*9941Ssam */ 9*9941Ssam 10*9941Ssam #include <errno.h> 11*9941Ssam #include <math.h> 12*9941Ssam 13*9941Ssam int errno; 14*9941Ssam static double invpi = 1.27323954473516268; 15*9941Ssam static double p0 = -0.1306820264754825668269611177e+5; 16*9941Ssam static double p1 = 0.1055970901714953193602353981e+4; 17*9941Ssam static double p2 = -0.1550685653483266376941705728e+2; 18*9941Ssam static double p3 = 0.3422554387241003435328470489e-1; 19*9941Ssam static double p4 = 0.3386638642677172096076369e-4; 20*9941Ssam static double q0 = -0.1663895238947119001851464661e+5; 21*9941Ssam static double q1 = 0.4765751362916483698926655581e+4; 22*9941Ssam static double q2 = -0.1555033164031709966900124574e+3; 23*9941Ssam 24*9941Ssam double tan(arg)25*9941Ssamtan(arg) 26*9941Ssam double arg; 27*9941Ssam { 28*9941Ssam double modf(); 29*9941Ssam double sign, temp, e, x, xsq; 30*9941Ssam int flag, i; 31*9941Ssam 32*9941Ssam flag = 0; 33*9941Ssam sign = 1.; 34*9941Ssam if(arg < 0.){ 35*9941Ssam arg = -arg; 36*9941Ssam sign = -1.; 37*9941Ssam } 38*9941Ssam arg = arg*invpi; /*overflow?*/ 39*9941Ssam x = modf(arg,&e); 40*9941Ssam i = e; 41*9941Ssam switch(i%4) { 42*9941Ssam case 1: 43*9941Ssam x = 1. - x; 44*9941Ssam flag = 1; 45*9941Ssam break; 46*9941Ssam 47*9941Ssam case 2: 48*9941Ssam sign = - sign; 49*9941Ssam flag = 1; 50*9941Ssam break; 51*9941Ssam 52*9941Ssam case 3: 53*9941Ssam x = 1. - x; 54*9941Ssam sign = - sign; 55*9941Ssam break; 56*9941Ssam 57*9941Ssam case 0: 58*9941Ssam break; 59*9941Ssam } 60*9941Ssam 61*9941Ssam xsq = x*x; 62*9941Ssam temp = ((((p4*xsq+p3)*xsq+p2)*xsq+p1)*xsq+p0)*x; 63*9941Ssam temp = temp/(((1.0*xsq+q2)*xsq+q1)*xsq+q0); 64*9941Ssam 65*9941Ssam if(flag == 1) { 66*9941Ssam if(temp == 0.) { 67*9941Ssam errno = ERANGE; 68*9941Ssam if (sign>0) 69*9941Ssam return(HUGE); 70*9941Ssam return(-HUGE); 71*9941Ssam } 72*9941Ssam temp = 1./temp; 73*9941Ssam } 74*9941Ssam return(sign*temp); 75*9941Ssam } 76