/*
   FUNCTIONS:	double sin(x)
		double cos(x)

   PURPOSE:
	To provide sine and cosine of angles expressed in radians.

   INPUT PARAMETERS:
	double x - The given angle in radians (maximum |x| is 3.37E+9)

   OUTPUT PARAMETERS:
	Returned - Sine(x) or Cosine(x), except where |x| > 3.37E+9.  In this
		   case, 0 is returned.
   DESCRIPTION:
	The algorithm uses a polynomial evaluation technique.  See algorithm
	3354 of Hart's "Computer Approximations".  Note: the cosine routine
	simply calls the sine routine via sin(x+pi/2).

   WRITTEN BY:  Scott D. Roth		DATE:  March 27, 1981
*/


#define ODD(x)		(((x) & 1) != 0)

static double coeff[9] = {
	 0.15707963267948966E1d,
	-0.64596409750624619E0d,
	 0.79692626246165436E-1d,
	-0.46817541353026426E-2d,
	 0.16044118470682207E-3d,
	-0.35988430072086927E-5d,
	 0.5692134872719023E-7d,
	-0.66843217206396E-9d,
	 0.587061098171E-11d};

double sin (x)			/* And "double cos (x)" follows */
    double x;
	{
	double abs_x, s, c;
	register double *c_ptr;
	long int_part;
	register int neg_x;			/* Boolean */
						/* For small x, x = sin(x) */
	if (1.0e-13d < (abs_x = ((neg_x = (x < 0)) ? -x : x)))
	    if (3.37E9d < abs_x)		/* Too big to figure */
		x = 0;
	    else
		{
		abs_x /= coeff[0];		/* x / (pi/2) */
		if (ODD((int) (int_part = abs_x)))
		    int_part++;
		c = abs_x - int_part;
		if (neg_x ^ ((((int) int_part) & 3) == 2))
		    c = -c;
		s = c * c;
		x = *(c_ptr = &coeff[8]);
		do    x = s * x + *--c_ptr;
		while (c_ptr != coeff);
		x *= c;
		}
	return(x);
	}	


double cos (x)
    double x;
	{
	return(sin(x + 1.57079632679489661));	/* sin(x + pi/2) */
	}
