/*
 *				s i n h . c
 */

/*)LIBRARY
*/

#ifdef	DOCUMENTATION
title	sinh	hyperbolic sine function
index	hyperbolic sine function
usage
	.s
	double x, f, sinh();
	.br
	f = sinh(x);
	.s
description
	.s
	Returns hyperbolic sine of argument.
	.s
diagnostics
	.s
	If the argument is large enough to cause overflow in the result
	the message 'sinh arg too large', 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 pp. 217-238 of Cody and Waite.  The constant v which
	appears in 2 define statements is chosen so that log(v) is an exact
	octal constant slightly larger than log(2) (See C and W pp. 218 and
	227.
	.s
author
	.s
	Hamish Ross.
	.s
date
	.s
	12-Jan-85
#endif

#include <math.h>

#define WMAX 88.0265309203708668	/* internal overflow limit	*/
#define CON1 0.6931610107421875		/* log(v)			*/
#define CON2 0.138302778796019026e-4	/* v / 2 - 1			*/

static double p0 = -0.351812834301771179e6;
static double p1 = -0.115635211968517683e5;
static double p2 = -0.163757982026307514e3;
static double p3 = -0.789661274173570995e0;
static double q0 = -0.211087700581062712e7;
static double q1 =  0.361627231094218365e5;
static double q2 = -0.277735231196507017e3;

double sinh(x)
double x;
{
    double f, p, q, y, exp();
    int sgn;

    if (x < 0.0) {
	y = -x;
	sgn = -1;
    }
    else {
	y = x;
	sgn = 1;
    }
    if (y > 1.0) {
	if (y > LOG_HUGE) {
	    y -= CON1;
	    if (y > WMAX) {
		cmemsg(FP_SINH, &x);
		return(sgn < 0 ? -HUGE : HUGE);
	    }
	    y = exp(y);
	    y += CON2;
	}
	else {
	    y = exp(y);
	    y = 0.5 * (y - 1.0 / y);
	}
    if (sgn < 0)
	y = -y;
    }
    else {
	if (y > ROOT_EPS) {
	    f = x * x;
	    p = (((p3 * f + p2) * f + p1) * f + p0) * f;
	    q = ((f + q2) * f + q1) * f + q0;
	    f = p / q;
	    y = x + x * f;
	}
    }
    return(y);
}