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