1*9938Ssam /* @(#)sin.c 4.1 12/25/82 */ 2*9938Ssam 3*9938Ssam /* 4*9938Ssam C program for floating point sin/cos. 5*9938Ssam Calls modf. 6*9938Ssam There are no error exits. 7*9938Ssam Coefficients are #3370 from Hart & Cheney (18.80D). 8*9938Ssam */ 9*9938Ssam 10*9938Ssam static double twoopi = 0.63661977236758134308; 11*9938Ssam static double p0 = .1357884097877375669092680e8; 12*9938Ssam static double p1 = -.4942908100902844161158627e7; 13*9938Ssam static double p2 = .4401030535375266501944918e6; 14*9938Ssam static double p3 = -.1384727249982452873054457e5; 15*9938Ssam static double p4 = .1459688406665768722226959e3; 16*9938Ssam static double q0 = .8644558652922534429915149e7; 17*9938Ssam static double q1 = .4081792252343299749395779e6; 18*9938Ssam static double q2 = .9463096101538208180571257e4; 19*9938Ssam static double q3 = .1326534908786136358911494e3; 20*9938Ssam 21*9938Ssam double cos(arg)22*9938Ssamcos(arg) 23*9938Ssam double arg; 24*9938Ssam { 25*9938Ssam double sinus(); 26*9938Ssam if(arg<0) 27*9938Ssam arg = -arg; 28*9938Ssam return(sinus(arg, 1)); 29*9938Ssam } 30*9938Ssam 31*9938Ssam double sin(arg)32*9938Ssamsin(arg) 33*9938Ssam double arg; 34*9938Ssam { 35*9938Ssam double sinus(); 36*9938Ssam return(sinus(arg, 0)); 37*9938Ssam } 38*9938Ssam 39*9938Ssam static double sinus(arg,quad)40*9938Ssamsinus(arg, quad) 41*9938Ssam double arg; 42*9938Ssam int quad; 43*9938Ssam { 44*9938Ssam double modf(); 45*9938Ssam double e, f; 46*9938Ssam double ysq; 47*9938Ssam double x,y; 48*9938Ssam int k; 49*9938Ssam double temp1, temp2; 50*9938Ssam 51*9938Ssam x = arg; 52*9938Ssam if(x<0) { 53*9938Ssam x = -x; 54*9938Ssam quad = quad + 2; 55*9938Ssam } 56*9938Ssam x = x*twoopi; /*underflow?*/ 57*9938Ssam if(x>32764){ 58*9938Ssam y = modf(x,&e); 59*9938Ssam e = e + quad; 60*9938Ssam modf(0.25*e,&f); 61*9938Ssam quad = e - 4*f; 62*9938Ssam }else{ 63*9938Ssam k = x; 64*9938Ssam y = x - k; 65*9938Ssam quad = (quad + k) & 03; 66*9938Ssam } 67*9938Ssam if (quad & 01) 68*9938Ssam y = 1-y; 69*9938Ssam if(quad > 1) 70*9938Ssam y = -y; 71*9938Ssam 72*9938Ssam ysq = y*y; 73*9938Ssam temp1 = ((((p4*ysq+p3)*ysq+p2)*ysq+p1)*ysq+p0)*y; 74*9938Ssam temp2 = ((((ysq+q3)*ysq+q2)*ysq+q1)*ysq+q0); 75*9938Ssam return(temp1/temp2); 76*9938Ssam } 77