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