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 <u.h>
9*3e12c5d1SDavid du Colombier #include <libc.h>
10*3e12c5d1SDavid du Colombier
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
sinus(double arg,int quad)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 *= 1/PIO2; /* 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
cos(double arg)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
sin(double arg)65*3e12c5d1SDavid du Colombier sin(double arg)
66*3e12c5d1SDavid du Colombier {
67*3e12c5d1SDavid du Colombier return sinus(arg, 0);
68*3e12c5d1SDavid du Colombier }
69