Subversion Repositories planix.SVN

Rev

Rev 2 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 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 <u.h>
9 
#include <libc.h>
10 
 
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
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
57 
cos(double arg)
58 
{
59 
	if(arg < 0)
60 
		arg = -arg;
61 
	return sinus(arg, 1);
62 
}
63 
 
64 
double
65 
sin(double arg)
66 
{
67 
	return sinus(arg, 0);
68 
}