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