Bell 32V development

[unix-history] / usr / src / libm / sin.c

/*
	C program for floating point sin/cos.
	Calls modf.
	There are no error exits.
	Coefficients are #3370 from Hart & Cheney (18.80D).
*/

static double twoopi	= 0.63661977236758134308;
static double p0	=  .1357884097877375669092680e8;
static double p1	= -.4942908100902844161158627e7;
static double p2	=  .4401030535375266501944918e6;
static double p3	= -.1384727249982452873054457e5;
static double p4	=  .1459688406665768722226959e3;
static double q0	=  .8644558652922534429915149e7;
static double q1	=  .4081792252343299749395779e6;
static double q2	=  .9463096101538208180571257e4;
static double q3	=  .1326534908786136358911494e3;

double
cos(arg)
double arg;
{
	double sinus();
	if(arg<0)
		arg = -arg;
	return(sinus(arg, 1));
}

double
sin(arg)
double arg;
{
	double sinus();
	return(sinus(arg, 0));
}

static double
sinus(arg, quad)
double arg;
int quad;
{
	double modf();
	double e, f;
	double ysq;
	double x,y;
	int k;
	double temp1, temp2;

	x = arg;
	if(x<0) {
		x = -x;
		quad = quad + 2;
	}
	x = x*twoopi;	/*underflow?*/
	if(x>32764){
		y = modf(x,&e);
		e = e + quad;
		modf(0.25*e,&f);
		quad = e - 4*f;
	}else{
		k = x;
		y = x - k;
		quad = (quad + k) & 03;
	}
	if (quad & 01)
		y = 1-y;
	if(quad > 1)
		y = -y;

	ysq = y*y;
	temp1 = ((((p4*ysq+p3)*ysq+p2)*ysq+p1)*ysq+p0)*y;
	temp2 = ((((ysq+q3)*ysq+q2)*ysq+q1)*ysq+q0);
	return(temp1/temp2);
}