lib/math/sinh.c

   1 /*
   2  * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
   3  * See the copyright notice in the ACK home directory, in the file "Copyright".
   4  *
   5  * Author: Ceriel J.H. Jacobs
   6  */
   7 /* $Header$ */
   8
   9 #include        <math.h>
  10 #include        <float.h>
  11 #include        <errno.h>
  12 #include        "localmath.h"
  13
  14 static double
  15 sinh_cosh(double x, int cosh_flag)
  16 {
  17         /*      Algorithm and coefficients from:
  18                         "Software manual for the elementary functions"
  19                         by W.J. Cody and W. Waite, Prentice-Hall, 1980
  20         */
  21
  22         static double p[] = {
  23                 -0.35181283430177117881e+6,
  24                 -0.11563521196851768270e+5,
  25                 -0.16375798202630751372e+3,
  26                 -0.78966127417357099479e+0
  27         };
  28         static double q[] = {
  29                 -0.21108770058106271242e+7,
  30                  0.36162723109421836460e+5,
  31                 -0.27773523119650701167e+3,
  32                  1.0
  33         };
  34         int     negative = x < 0;
  35         double  y = negative ? -x : x;
  36
  37         if (__IsNan(x)) {
  38                 errno = EDOM;
  39                 return x;
  40         }
  41         if (! cosh_flag && y <= 1.0) {
  42                 /* ??? check for underflow ??? */
  43                 y = y * y;
  44                 return x + x * y * POLYNOM3(y, p)/POLYNOM3(y,q);
  45         }
  46
  47         if (y >= M_LN_MAX_D) {
  48                 /* exp(y) would cause overflow */
  49 #define LNV     0.69316101074218750000e+0
  50 #define VD2M1   0.52820835025874852469e-4
  51                 double  w = y - LNV;
  52
  53                 if (w < M_LN_MAX_D+M_LN2-LNV) {
  54                         x = exp(w);
  55                         x += VD2M1 * x;
  56                 }
  57                 else {
  58                         errno = ERANGE;
  59                         x = HUGE_VAL;
  60                 }
  61         }
  62         else {
  63                 double  z = exp(y);
  64
  65                 x = 0.5 * (z + (cosh_flag ? 1.0 : -1.0)/z);
  66         }
  67         return negative ? -x : x;
  68 }
  69
  70 double
  71 sinh(double x)
  72 {
  73         return sinh_cosh(x, 0);
  74 }
  75
  76 double
  77 cosh(double x)
  78 {
  79         if (x < 0) x = -x;
  80         return sinh_cosh(x, 1);
  81 }