/* sin - sine and cosine functions */

/*
Based on the Yale CS "C" library.  Coefficients are #3370 from
Hart and Cheney (18.80D)

Calls MODF
*/
extern	double	modf();		/* floating point modulus		*/
extern	double	itof();		/* integer to floating conversion	*/

static	double	twoopi	= 0.63661977236758134308;
static	double	p0	= 0.1357884097877375669092680e8;
static	double	p1	= -.4942908100902844161158627e7;
static	double	p2	= 0.4401030535375266501944918e6;
static	double	p3	= -.1384727249982452873054457e5;
static	double	p4	= 0.1459688406665768722226959e3;
static	double	q0	= 0.8644558652922534429915149e7;
static	double	q1	= 0.4081792252343299749395779e6;
static	double	q2	= 0.9463096101538208180571257e4;
static	double	q3	= 0.1326534908786136358911494e3;
/* sinus - does the work for sine and cosine */

static	double	sinus(arg,quad)

	double	arg;
	int	quad;
{
	double	e,f;
	double	ysq,x,y;
	int	k;
	double	temp1,temp2;

	x = arg;
	if (x < 0)
	{
	    x = -x;
	    quad += 2;
	}
	x = x*twoopi;
	if (x > 32764.)
	{
	    y = modf(x,&e);
	    e = e+ftoi(quad);
	    modf(0.25*e,&f);
	    quad = ftoi(e)-ftoi(4.0*f);
	}
	else
	{
	    k = ftoi(x);
	    y = x-itof(k);
	    quad = (quad+k)&3;
	}
	if (quad&01) y = 1.0-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);
}
/* the sine and cosine calls */

double	cos(arg)

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

double	sin(arg)

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