newlib/libm/mathfp/sf_sineh.c

   1
   2 /* @(#)z_sinehf.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  * Hyperbolic Sine
  11  *
  12  * Input:
  13  *   x - floating point value
  14  *
  15  * Output:
  16  *   hyperbolic sine of x
  17  *
  18  * Description:
  19  *   This routine calculates hyperbolic sines.
  20  *
  21  *****************************************************************/
  22
  23 #include <float.h>
  24 #include "fdlibm.h"
  25 #include "zmath.h"
  26
  27 static const float q[] = { -0.428277109e+2 };
  28 static const float p[] = { -0.713793159e+1,
  29                            -0.190333399 };
  30 static const float LNV = 0.6931610107;
  31 static const float INV_V2 = 0.2499930850;
  32 static const float V_OVER2_MINUS1 = 0.1383027787e-4;
  33
  34 float
  35 sinehf (float x,
  36         int cosineh)
  37 {
  38   float y, f, P, Q, R, res, z, w;
  39   int sgn = 1;
  40   float WBAR = 18.55;
  41
  42   /* Check for special values. */
  43   switch (numtestf (x))
  44     {
  45       case NAN:
  46         errno = EDOM;
  47         return (x);
  48       case INF:
  49         errno = ERANGE;
  50         return (ispos (x) ? z_infinity_f.f : -z_infinity_f.f);
  51     }
  52
  53   y = fabs (x);
  54
  55   if (!cosineh && x < 0.0)
  56     sgn = -1;
  57
  58   if ((y > 1.0 && !cosineh) || cosineh)
  59     {
  60       if (y > BIGX)
  61         {
  62           w = y - LNV;
  63
  64           /* Check for w > maximum here. */
  65           if (w > BIGX)
  66             {
  67               errno = ERANGE;
  68               return (x);
  69             }
  70
  71           z = exp (w);
  72
  73           if (w > WBAR)
  74             res = z * (V_OVER2_MINUS1 + 1.0);
  75         }
  76
  77       else
  78         {
  79           z = exp (y);
  80           if (cosineh)
  81             res = (z + 1 / z) / 2.0;
  82           else
  83             res = (z - 1 / z) / 2.0;
  84         }
  85
  86       if (sgn < 0)
  87         res = -res;
  88     }
  89   else
  90     {
  91       /* Check for y being too small. */
  92       if (y < z_rooteps_f)
  93         {
  94           res = x;
  95         }
  96       /* Calculate the Taylor series. */
  97       else
  98         {
  99           f = x * x;
 100           Q = f + q[0];
 101           P = p[1] * f + p[0];
 102           R = f * (P / Q);
 103
 104           res = x + x * R;
 105         }
 106     }
 107
 108   return (res);
 109 }