/*
 *				s i n c o s . c
 */

/*)LIBRARY
*/

#ifdef	DOCUMENTATION
title	sincos	sin and cosine functions
index	sin function
index	cosine function
usage
	.s
	double x, f, sin();
	.br
	double g, cos();
	.br
	f = sin(x);
	.br
	g = cos(x);
	.s
description
	.s
	sin returns sin of x.
	.br
	cos returns cosine of x.
	.s
diagnostics
	.s
	If the absolute value of the argument is greater than TRIG_MAX the
	error message 'sin, cos arg too large', followed by the value of
	the argument, is written to stderr.  The
	algorithm returns a result of lower accuracy.  The number of reliable
	decimal digits will be roughly 20 - log10(|x / pi|).
	.s
internal
	.s
	Algorithm from Cody and Waite pp. 125-149.  A small modification to
	the published algorithm has meant that the value of the argument
	is not restricted to values whose integer part can be represented
	by a long integer.  The accuracy of the result in those cases will,
	however, be very low (c.f. above).
	.s
author
	.s
	Hamish Ross.
	.s
date
	.s
	8-Jan-85
#endif

#include <math.h>

double sin(x)			
double x;
{
    double sin_cos();

    if (x < 0.0)
	return(sin_cos(x, -x, -1.0));
    else
	return(sin_cos(x, x, 1.0));
}

double cos(x)
double x;
{
    double fabs(), sin_cos();

    return(sin_cos(x, fabs(x) + HALF_PI, 1.0));
}

#define C1 3.1416015625			/* C1 + C2 = pi to 8 bits more	*/
#define C2 -8.90891020676153736e-6	/* than machine precision	*/

static double r1 = -0.166666666666666651e0;
static double r2 =  0.833333333333316503e-2;
static double r3 = -0.198412698412018405e-3;
static double r4 =  0.275573192101527561e-5;
static double r5 = -0.250521067982745845e-7;
static double r6 =  0.160589364903715891e-9;
static double r7 = -0.764291780689104677e-12;
static double r8 =  0.272047909578888462e-14;
static double sin_cos(x, y, sgn)
double x, y, sgn;
{
    double f, g, modf(), r, xn;

    if (y > TRIG_MAX) {
	cmemsg(FP_TRIG, &x);
    }
    modf(REC_PI * y + 0.5, &xn);
    if ((int)(2.0 * modf(0.5 * xn, &r)) & 1 == 1)
	sgn = -sgn;
    r = x < 0.0 ? -x : x;
    if (r != y)
	xn -= 0.5;
    f = (r - xn * C1) - xn * C2;
    f *= sgn;
    if (f < ROOT_EPS && f > -ROOT_EPS)
	return(f);
    g = f * f;
    r = (((((((r8*g + r7)*g + r6)*g + r5)*g + r4)*g + r3)*g + r2)*g + r1)*g;
    return(f + f * r);
}