13e12c5d1SDavid du Colombier /*
23e12c5d1SDavid du Colombier C program for floating point sin/cos.
33e12c5d1SDavid du Colombier Calls modf.
43e12c5d1SDavid du Colombier There are no error exits.
53e12c5d1SDavid du Colombier Coefficients are #3370 from Hart & Cheney (18.80D).
63e12c5d1SDavid du Colombier */
73e12c5d1SDavid du Colombier
83e12c5d1SDavid du Colombier #include <math.h>
93e12c5d1SDavid du Colombier
10*7dd7cddfSDavid du Colombier #define PIO2 1.570796326794896619231e0
113e12c5d1SDavid du Colombier #define p0 .1357884097877375669092680e8
123e12c5d1SDavid du Colombier #define p1 -.4942908100902844161158627e7
133e12c5d1SDavid du Colombier #define p2 .4401030535375266501944918e6
143e12c5d1SDavid du Colombier #define p3 -.1384727249982452873054457e5
153e12c5d1SDavid du Colombier #define p4 .1459688406665768722226959e3
163e12c5d1SDavid du Colombier #define q0 .8644558652922534429915149e7
173e12c5d1SDavid du Colombier #define q1 .4081792252343299749395779e6
183e12c5d1SDavid du Colombier #define q2 .9463096101538208180571257e4
193e12c5d1SDavid du Colombier #define q3 .1326534908786136358911494e3
203e12c5d1SDavid du Colombier
213e12c5d1SDavid du Colombier static
223e12c5d1SDavid du Colombier double
sinus(double arg,int quad)233e12c5d1SDavid du Colombier sinus(double arg, int quad)
243e12c5d1SDavid du Colombier {
253e12c5d1SDavid du Colombier double e, f, ysq, x, y, temp1, temp2;
263e12c5d1SDavid du Colombier int k;
273e12c5d1SDavid du Colombier
283e12c5d1SDavid du Colombier x = arg;
293e12c5d1SDavid du Colombier if(x < 0) {
303e12c5d1SDavid du Colombier x = -x;
313e12c5d1SDavid du Colombier quad += 2;
323e12c5d1SDavid du Colombier }
33*7dd7cddfSDavid du Colombier x *= 1/PIO2; /* underflow? */
343e12c5d1SDavid du Colombier if(x > 32764) {
353e12c5d1SDavid du Colombier y = modf(x, &e);
363e12c5d1SDavid du Colombier e += quad;
373e12c5d1SDavid du Colombier modf(0.25*e, &f);
383e12c5d1SDavid du Colombier quad = e - 4*f;
393e12c5d1SDavid du Colombier } else {
403e12c5d1SDavid du Colombier k = x;
413e12c5d1SDavid du Colombier y = x - k;
423e12c5d1SDavid du Colombier quad += k;
433e12c5d1SDavid du Colombier quad &= 3;
443e12c5d1SDavid du Colombier }
453e12c5d1SDavid du Colombier if(quad & 1)
463e12c5d1SDavid du Colombier y = 1-y;
473e12c5d1SDavid du Colombier if(quad > 1)
483e12c5d1SDavid du Colombier y = -y;
493e12c5d1SDavid du Colombier
503e12c5d1SDavid du Colombier ysq = y*y;
513e12c5d1SDavid du Colombier temp1 = ((((p4*ysq+p3)*ysq+p2)*ysq+p1)*ysq+p0)*y;
523e12c5d1SDavid du Colombier temp2 = ((((ysq+q3)*ysq+q2)*ysq+q1)*ysq+q0);
533e12c5d1SDavid du Colombier return temp1/temp2;
543e12c5d1SDavid du Colombier }
553e12c5d1SDavid du Colombier
563e12c5d1SDavid du Colombier double
cos(double arg)573e12c5d1SDavid du Colombier cos(double arg)
583e12c5d1SDavid du Colombier {
593e12c5d1SDavid du Colombier if(arg < 0)
603e12c5d1SDavid du Colombier arg = -arg;
613e12c5d1SDavid du Colombier return sinus(arg, 1);
623e12c5d1SDavid du Colombier }
633e12c5d1SDavid du Colombier
643e12c5d1SDavid du Colombier double
sin(double arg)653e12c5d1SDavid du Colombier sin(double arg)
663e12c5d1SDavid du Colombier {
673e12c5d1SDavid du Colombier return sinus(arg, 0);
683e12c5d1SDavid du Colombier }
69