1*3e12c5d1SDavid du Colombier #include <u.h> 2*3e12c5d1SDavid du Colombier #include <libc.h> 3*3e12c5d1SDavid du Colombier 4*3e12c5d1SDavid du Colombier /* 5*3e12c5d1SDavid du Colombier * sinh(arg) returns the hyperbolic sine of its floating- 6*3e12c5d1SDavid du Colombier * point argument. 7*3e12c5d1SDavid du Colombier * 8*3e12c5d1SDavid du Colombier * The exponential function is called for arguments 9*3e12c5d1SDavid du Colombier * greater in magnitude than 0.5. 10*3e12c5d1SDavid du Colombier * 11*3e12c5d1SDavid du Colombier * A series is used for arguments smaller in magnitude than 0.5. 12*3e12c5d1SDavid du Colombier * The coefficients are #2029 from Hart & Cheney. (20.36D) 13*3e12c5d1SDavid du Colombier * 14*3e12c5d1SDavid du Colombier * cosh(arg) is computed from the exponential function for 15*3e12c5d1SDavid du Colombier * all arguments. 16*3e12c5d1SDavid du Colombier */ 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 Colombiersinh(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 = 0; 33*3e12c5d1SDavid du Colombier if(arg < 0) { 34*3e12c5d1SDavid du Colombier arg = -arg; 35*3e12c5d1SDavid du Colombier sign++; 36*3e12c5d1SDavid du Colombier } 37*3e12c5d1SDavid du Colombier if(arg > 21) { 38*3e12c5d1SDavid du Colombier temp = exp(arg)/2; 39*3e12c5d1SDavid du Colombier goto out; 40*3e12c5d1SDavid du Colombier } 41*3e12c5d1SDavid du Colombier if(arg > 0.5) { 42*3e12c5d1SDavid du Colombier temp = (exp(arg) - exp(-arg))/2; 43*3e12c5d1SDavid du Colombier goto out; 44*3e12c5d1SDavid du Colombier } 45*3e12c5d1SDavid du Colombier argsq = arg*arg; 46*3e12c5d1SDavid du Colombier temp = (((p3*argsq+p2)*argsq+p1)*argsq+p0)*arg; 47*3e12c5d1SDavid du Colombier temp /= (((argsq+q2)*argsq+q1)*argsq+q0); 48*3e12c5d1SDavid du Colombier out: 49*3e12c5d1SDavid du Colombier if(sign) 50*3e12c5d1SDavid du Colombier temp = -temp; 51*3e12c5d1SDavid du Colombier return temp; 52*3e12c5d1SDavid du Colombier } 53*3e12c5d1SDavid du Colombier 54*3e12c5d1SDavid du Colombier double cosh(double arg)55*3e12c5d1SDavid du Colombiercosh(double arg) 56*3e12c5d1SDavid du Colombier { 57*3e12c5d1SDavid du Colombier if(arg < 0) 58*3e12c5d1SDavid du Colombier arg = - arg; 59*3e12c5d1SDavid du Colombier if(arg > 21) 60*3e12c5d1SDavid du Colombier return exp(arg)/2; 61*3e12c5d1SDavid du Colombier return (exp(arg) + exp(-arg))/2; 62*3e12c5d1SDavid du Colombier } 63