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