/*
 *				t a n c o t . c
 */

/*)LIBRARY
*/

#ifdef	DOCUMENTATION
title	tancot	tangent and cotangent functions
index	tangent function
index	cotangent function
usage
	.s
	double x, f, tan();
	double g, cotan();
	.br
	f = tan(x);
	g = cotan(x);
	.s
description
	.s
	tan returns tan of x.
	cotan returns cotangent of x.
	.s
diagnostics
	.s
	If the absolute value of the argument is greater than TRIG_MAX the
	error message 'tan, cotan 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|).
	.br
	If the argument of cotan is less in magnitude than 1 / HUGE
	the error message 'cotan arg too small', followed by the value of the
	argument, is written to stderr.  A value of HUGE with the same sign as
	the argument is returned.
	.s
internal
	.s
	Algorithm from Cody and Waite pp. 150-173.  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>

#define C1 1.57080078125	/* C1 + C2 represent pi / 2 to 12 more	*/
#define C2 -4.45445510338076868e-6	/* bits than machine precision	*/

static int flag;
double tan(x)			
double x;
{
    double tan_cot();

    flag = 0;
    return(tan_cot(x));
}

double cotan(x)
double x;
{
    double tan_cot();

    flag = 1;
    if ((x < REC_HUGE) && (x > -REC_HUGE)) {
	cmemsg(FP_COTE, &x);
	return(x < 0.0 ? -HUGE : HUGE);
    }
    return(tan_cot(x));
}

static double p1 = -0.133383500064219607e0;
static double p2 =  0.342488782358905900e-2;
static double p3 = -0.178617073422544267e-4;
static double q1 = -0.466716833397552942e0; 
static double q2 =  0.256638322894401129e-1;
static double q3 = -0.311815319070100273e-3;
static double q4 =  0.498194339937865123e-6;
static double tan_cot(x)
double x;
{
    double f, g, p, q;
    double modf();
    int n;

    f = x < 0.0 ? -x : x;
    if (f > TRIG_MAX / 2.0)
	cmemsg(FP_TANE, &x);
    if (x < 0.0)
	modf(x * RPIBY2 - 0.5, &g);
    else
	modf(x * RPIBY2 + 0.5, &g);
    n = 2.0 * modf(0.5 * g, &p);
    f = (x - g * C1) - g * C2;
    if (f < ROOT_EPS && f > -ROOT_EPS) {
	p = f;
	q = 1.0;
    }
    else {
	g = f * f;
	p = ((p3 * g + p2) * g + p1) * g * f + f;
	q = (((q4 * g + q3) * g + q2) * g + q1) * g + 1.0;
    }
    if (flag == 0) {
	if (n == 0)
	    return(p / q);
	else
	    return(-q / p);
    }
    else {
	if (n == 0)
	    return(q / p);
	else
	    return(-p / q);
    }
}