newlib/libm/mathfp/s_sineh.c

   1
   2 /* @(#)z_sineh.c 1.0 98/08/13 */
   3 /******************************************************************
   4  * The following routines are coded directly from the algorithms
   5  * and coefficients given in "Software Manual for the Elementary
   6  * Functions" by William J. Cody, Jr. and William Waite, Prentice
   7  * Hall, 1980.
   8  ******************************************************************/
   9
  10 /*
  11 FUNCTION
  12         <<sinh>>, <<sinhf>>, <<cosh>>, <<coshf>>, <<sineh>>---hyperbolic sine or cosine
  13
  14 INDEX
  15         sinh
  16 INDEX
  17         sinhf
  18 INDEX
  19         cosh
  20 INDEX
  21         coshf
  22
  23 SYNOPSIS
  24         #include <math.h>
  25         double sinh(double <[x]>);
  26         float  sinhf(float <[x]>);
  27         double cosh(double <[x]>);
  28         float  coshf(float <[x]>);
  29
  30 DESCRIPTION
  31         <<sinh>> and <<cosh>> compute the hyperbolic sine or cosine
  32         of the argument <[x]>.
  33         Angles are specified in radians.   <<sinh>>(<[x]>) is defined as
  34         @ifnottex
  35         . (exp(<[x]>) - exp(-<[x]>))/2
  36         @end ifnottex
  37         @tex
  38         $${e^x - e^{-x}}\over 2$$
  39         @end tex
  40         <<cosh>> is defined as
  41         @ifnottex
  42         . (exp(<[x]>) - exp(-<[x]>))/2
  43         @end ifnottex
  44         @tex
  45         $${e^x + e^{-x}}\over 2$$
  46         @end tex
  47
  48         <<sinhf>> and <<coshf>> are identical, save that they take
  49         and returns <<float>> values.
  50
  51 RETURNS
  52         The hyperbolic sine or cosine of <[x]> is returned.
  53
  54         When the correct result is too large to be representable (an
  55         overflow),  the functions return <<HUGE_VAL>> with the
  56         appropriate sign, and sets the global value <<errno>> to
  57         <<ERANGE>>.
  58
  59 PORTABILITY
  60         <<sinh>> is ANSI C.
  61         <<sinhf>> is an extension.
  62         <<cosh>> is ANSI C.
  63         <<coshf>> is an extension.
  64
  65 */
  66
  67 /******************************************************************
  68  * Hyperbolic Sine
  69  *
  70  * Input:
  71  *   x - floating point value
  72  *
  73  * Output:
  74  *   hyperbolic sine of x
  75  *
  76  * Description:
  77  *   This routine calculates hyperbolic sines.
  78  *
  79  *****************************************************************/
  80
  81 #include <float.h>
  82 #include "fdlibm.h"
  83 #include "zmath.h"
  84
  85 static const double q[] = { -0.21108770058106271242e+7,
  86                              0.36162723109421836460e+5,
  87                             -0.27773523119650701667e+3 };
  88 static const double p[] = { -0.35181283430177117881e+6,
  89                             -0.11563521196851768270e+5,
  90                             -0.16375798202630751372e+3,
  91                             -0.78966127417357099479 };
  92 static const double LNV = 0.6931610107421875000;
  93 static const double INV_V2 = 0.24999308500451499336;
  94 static const double V_OVER2_MINUS1 = 0.13830277879601902638e-4;
  95
  96 double
  97 sineh (double x,
  98         int cosineh)
  99 {
 100   double y, f, P, Q, R, res, z, w;
 101   int sgn = 1;
 102   double WBAR = 18.55;
 103
 104   /* Check for special values. */
 105   switch (numtest (x))
 106     {
 107       case NAN:
 108         errno = EDOM;
 109         return (x);
 110       case INF:
 111         errno = ERANGE;
 112         return (ispos (x) ? z_infinity.d : -z_infinity.d);
 113     }
 114
 115   y = fabs (x);
 116
 117   if (!cosineh && x < 0.0)
 118     sgn = -1;
 119
 120   if ((y > 1.0 && !cosineh) || cosineh)
 121     {
 122       if (y > BIGX)
 123         {
 124           w = y - LNV;
 125
 126           /* Check for w > maximum here. */
 127           if (w > BIGX)
 128             {
 129               errno = ERANGE;
 130               return (x);
 131             }
 132
 133           z = exp (w);
 134
 135           if (w > WBAR)
 136             res = z * (V_OVER2_MINUS1 + 1.0);
 137         }
 138
 139       else
 140         {
 141           z = exp (y);
 142           if (cosineh)
 143             res = (z + 1 / z) / 2.0;
 144           else
 145             res = (z - 1 / z) / 2.0;
 146         }
 147
 148       if (sgn < 0)
 149         res = -res;
 150     }
 151   else
 152     {
 153       /* Check for y being too small. */
 154       if (y < z_rooteps)
 155         {
 156           res = x;
 157         }
 158       /* Calculate the Taylor series. */
 159       else
 160         {
 161           f = x * x;
 162           Q = ((f + q[2]) * f + q[1]) * f + q[0];
 163           P = ((p[3] * f + p[2]) * f + p[1]) * f + p[0];
 164           R = f * (P / Q);
 165
 166           res = x + x * R;
 167         }
 168     }
 169
 170   return (res);
 171 }