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
sinus(double arg,int quad)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
cos(double arg)57 cos(double arg)
58 {
59 if(arg < 0)
60 arg = -arg;
61 return sinus(arg, 1);
62 }
63
64 double
sin(double arg)65 sin(double arg)
66 {
67 return sinus(arg, 0);
68 }
69