sys/src/ape/lib/ap/math/sinh.c

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