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