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 <u.h> 9*3e12c5d1SDavid du Colombier #include <libc.h> 10*3e12c5d1SDavid du Colombier 11*3e12c5d1SDavid du Colombier static double p0 = -0.1306820264754825668269611177e+5; 12*3e12c5d1SDavid du Colombier static double p1 = 0.1055970901714953193602353981e+4; 13*3e12c5d1SDavid du Colombier static double p2 = -0.1550685653483266376941705728e+2; 14*3e12c5d1SDavid du Colombier static double p3 = 0.3422554387241003435328470489e-1; 15*3e12c5d1SDavid du Colombier static double p4 = 0.3386638642677172096076369e-4; 16*3e12c5d1SDavid du Colombier static double q0 = -0.1663895238947119001851464661e+5; 17*3e12c5d1SDavid du Colombier static double q1 = 0.4765751362916483698926655581e+4; 18*3e12c5d1SDavid du Colombier static double q2 = -0.1555033164031709966900124574e+3; 19*3e12c5d1SDavid du Colombier 20*3e12c5d1SDavid du Colombier double tan(double arg)21*3e12c5d1SDavid du Colombiertan(double arg) 22*3e12c5d1SDavid du Colombier { 23*3e12c5d1SDavid du Colombier double temp, e, x, xsq; 24*3e12c5d1SDavid du Colombier int flag, sign, i; 25*3e12c5d1SDavid du Colombier 26*3e12c5d1SDavid du Colombier flag = 0; 27*3e12c5d1SDavid du Colombier sign = 0; 28*3e12c5d1SDavid du Colombier if(arg < 0){ 29*3e12c5d1SDavid du Colombier arg = -arg; 30*3e12c5d1SDavid du Colombier sign++; 31*3e12c5d1SDavid du Colombier } 32*3e12c5d1SDavid du Colombier arg = 2*arg/PIO2; /* overflow? */ 33*3e12c5d1SDavid du Colombier x = modf(arg, &e); 34*3e12c5d1SDavid du Colombier i = e; 35*3e12c5d1SDavid du Colombier switch(i%4) { 36*3e12c5d1SDavid du Colombier case 1: 37*3e12c5d1SDavid du Colombier x = 1 - x; 38*3e12c5d1SDavid du Colombier flag = 1; 39*3e12c5d1SDavid du Colombier break; 40*3e12c5d1SDavid du Colombier 41*3e12c5d1SDavid du Colombier case 2: 42*3e12c5d1SDavid du Colombier sign = !sign; 43*3e12c5d1SDavid du Colombier flag = 1; 44*3e12c5d1SDavid du Colombier break; 45*3e12c5d1SDavid du Colombier 46*3e12c5d1SDavid du Colombier case 3: 47*3e12c5d1SDavid du Colombier x = 1 - x; 48*3e12c5d1SDavid du Colombier sign = !sign; 49*3e12c5d1SDavid du Colombier break; 50*3e12c5d1SDavid du Colombier 51*3e12c5d1SDavid du Colombier case 0: 52*3e12c5d1SDavid du Colombier break; 53*3e12c5d1SDavid du Colombier } 54*3e12c5d1SDavid du Colombier 55*3e12c5d1SDavid du Colombier xsq = x*x; 56*3e12c5d1SDavid du Colombier temp = ((((p4*xsq+p3)*xsq+p2)*xsq+p1)*xsq+p0)*x; 57*3e12c5d1SDavid du Colombier temp = temp/(((xsq+q2)*xsq+q1)*xsq+q0); 58*3e12c5d1SDavid du Colombier 59*3e12c5d1SDavid du Colombier if(flag) { 60*3e12c5d1SDavid du Colombier if(temp == 0) 61*3e12c5d1SDavid du Colombier return NaN(); 62*3e12c5d1SDavid du Colombier temp = 1/temp; 63*3e12c5d1SDavid du Colombier } 64*3e12c5d1SDavid du Colombier if(sign) 65*3e12c5d1SDavid du Colombier temp = -temp; 66*3e12c5d1SDavid du Colombier return temp; 67*3e12c5d1SDavid du Colombier } 68