/********************************************************************** * Copyright (c) Digital Equipment Corporation 1984, 1985, 1986. * * All Rights Reserved. * * Reference "/usr/src/COPYRIGHT" for applicable restrictions. * **********************************************************************/ /* SCCSID: @(#)sinh.c 3.0 4/22/86 */ /* (System 5) sinh.c 1.9 */ /*LINTLIBRARY*/ /* * sinh returns the hyperbolic sine of its double-precision argument. * A series is used for arguments smaller in magnitude than 1. * The exponential function is used for arguments * greater in magnitude than 1. * * cosh returns the hyperbolic cosine of its double-precision argument. * cosh is computed from the exponential function for * all arguments. * * Returns ERANGE error and value HUGE (or -HUGE for sinh of * negative argument) if the correct result would overflow. * * Algorithm and coefficients from Cody and Waite (1980). */ #include <math.h> #include <values.h> #include <errno.h> #define X_MAX (LN_MAXDOUBLE + M_LN2) #define LNV 0.6931610107421875 #define V2M1 0.13830277879601902638e-4 extern double sinh_exc(); double sinh(x) register double x; { register double y = _ABS(x); if (y <= 1) { static double p[] = { -0.78966127417357099479e0, -0.16375798202630751372e3, -0.11563521196851768270e5, -0.35181283430177117881e6, }, q[] = { 1.0, -0.27773523119650701667e3, 0.36162723109421836460e5, -0.21108770058106271242e7, }; if (y < X_EPS) /* for efficiency and to prevent underflow */ return (x); y = x * x; return (x + x * y * _POLY3(y, p)/_POLY3(y, q)); } if (y > LN_MAXDOUBLE) /* exp(x) would overflow */ return (sinh_exc(x, y, 1)); x = exp(x); return (0.5 * (x - 1.0/x)); } double cosh(x) register double x; { register double y = _ABS(x); if (y > LN_MAXDOUBLE) /* exp(x) would overflow */ return (sinh_exc(x, y, 0)); x = exp(y); return (0.5 * (x + 1.0/x)); } static double sinh_exc(x, y, sinhflag) register double x, y; register int sinhflag; { int neg = (x < 0 && sinhflag); /* sinh of negative argument */ struct exception exc; if (y < X_MAX) { /* result is still representable */ x = exp(y - LNV); x += V2M1 * x; return (neg ? -x : x); } exc.type = OVERFLOW; exc.name = sinhflag ? "sinh" : "cosh"; exc.arg1 = x; exc.retval = neg ? -HUGE : HUGE; if (!matherr(&exc)) errno = ERANGE; return (exc.retval); }