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