xref: /plan9/sys/src/ape/lib/ap/math/sinh.c (revision 3e12c5d1bb89fc02707907988834ef147769ddaf)
1*3e12c5d1SDavid du Colombier /*
2*3e12c5d1SDavid du Colombier 	sinh(arg) returns the hyperbolic sine of its floating-
3*3e12c5d1SDavid du Colombier 	point argument.
4*3e12c5d1SDavid du Colombier 
5*3e12c5d1SDavid du Colombier 	The exponential function is called for arguments
6*3e12c5d1SDavid du Colombier 	greater in magnitude than 0.5.
7*3e12c5d1SDavid du Colombier 
8*3e12c5d1SDavid du Colombier 	A series is used for arguments smaller in magnitude than 0.5.
9*3e12c5d1SDavid du Colombier 	The coefficients are #2029 from Hart & Cheney. (20.36D)
10*3e12c5d1SDavid du Colombier 
11*3e12c5d1SDavid du Colombier 	cosh(arg) is computed from the exponential function for
12*3e12c5d1SDavid du Colombier 	all arguments.
13*3e12c5d1SDavid du Colombier */
14*3e12c5d1SDavid du Colombier 
15*3e12c5d1SDavid du Colombier #include <math.h>
16*3e12c5d1SDavid du Colombier #include <errno.h>
17*3e12c5d1SDavid du Colombier 
18*3e12c5d1SDavid du Colombier static double p0  = -0.6307673640497716991184787251e+6;
19*3e12c5d1SDavid du Colombier static double p1  = -0.8991272022039509355398013511e+5;
20*3e12c5d1SDavid du Colombier static double p2  = -0.2894211355989563807284660366e+4;
21*3e12c5d1SDavid du Colombier static double p3  = -0.2630563213397497062819489e+2;
22*3e12c5d1SDavid du Colombier static double q0  = -0.6307673640497716991212077277e+6;
23*3e12c5d1SDavid du Colombier static double q1   = 0.1521517378790019070696485176e+5;
24*3e12c5d1SDavid du Colombier static double q2  = -0.173678953558233699533450911e+3;
25*3e12c5d1SDavid du Colombier 
26*3e12c5d1SDavid du Colombier double
sinh(double arg)27*3e12c5d1SDavid du Colombier sinh(double arg)
28*3e12c5d1SDavid du Colombier {
29*3e12c5d1SDavid du Colombier 	double temp, argsq;
30*3e12c5d1SDavid du Colombier 	int sign;
31*3e12c5d1SDavid du Colombier 
32*3e12c5d1SDavid du Colombier 	sign = 1;
33*3e12c5d1SDavid du Colombier 	if(arg < 0) {
34*3e12c5d1SDavid du Colombier 		arg = - arg;
35*3e12c5d1SDavid du Colombier 		sign = -1;
36*3e12c5d1SDavid du Colombier 	}
37*3e12c5d1SDavid du Colombier 	if(arg > 21) {
38*3e12c5d1SDavid du Colombier 		if(arg >= HUGE_VAL){
39*3e12c5d1SDavid du Colombier 			errno = ERANGE;
40*3e12c5d1SDavid du Colombier 			temp = HUGE_VAL;
41*3e12c5d1SDavid du Colombier 		} else
42*3e12c5d1SDavid du Colombier 			temp = exp(arg)/2;
43*3e12c5d1SDavid du Colombier 		if(sign > 0)
44*3e12c5d1SDavid du Colombier 			return temp;
45*3e12c5d1SDavid du Colombier 		else
46*3e12c5d1SDavid du Colombier 			return -temp;
47*3e12c5d1SDavid du Colombier 	}
48*3e12c5d1SDavid du Colombier 
49*3e12c5d1SDavid du Colombier 	if(arg > 0.5)
50*3e12c5d1SDavid du Colombier 		return sign*(exp(arg) - exp(-arg))/2;
51*3e12c5d1SDavid du Colombier 
52*3e12c5d1SDavid du Colombier 	argsq = arg*arg;
53*3e12c5d1SDavid du Colombier 	temp = (((p3*argsq+p2)*argsq+p1)*argsq+p0)*arg;
54*3e12c5d1SDavid du Colombier 	temp /= (((argsq+q2)*argsq+q1)*argsq+q0);
55*3e12c5d1SDavid du Colombier 	return sign*temp;
56*3e12c5d1SDavid du Colombier }
57*3e12c5d1SDavid du Colombier 
58*3e12c5d1SDavid du Colombier double
cosh(double arg)59*3e12c5d1SDavid du Colombier cosh(double arg)
60*3e12c5d1SDavid du Colombier {
61*3e12c5d1SDavid du Colombier 	if(arg < 0)
62*3e12c5d1SDavid du Colombier 		arg = - arg;
63*3e12c5d1SDavid du Colombier 	if(arg > 21) {
64*3e12c5d1SDavid du Colombier 		if(arg >= HUGE_VAL){
65*3e12c5d1SDavid du Colombier 			errno = ERANGE;
66*3e12c5d1SDavid du Colombier 			return HUGE_VAL;
67*3e12c5d1SDavid du Colombier 		} else
68*3e12c5d1SDavid du Colombier 			return(exp(arg)/2);
69*3e12c5d1SDavid du Colombier 	}
70*3e12c5d1SDavid du Colombier 	return (exp(arg) + exp(-arg))/2;
71*3e12c5d1SDavid du Colombier }
72