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