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 _DEFUN (sinehf, (float, int),
36         float x _AND
37         int cosineh)
38 {
39   float y, f, P, Q, R, res, z, w;
40   int sgn = 1;
41   float WBAR = 18.55;
42 
43   /* Check for special values. */
44   switch (numtestf (x))
45     {
46       case NAN:
47         errno = EDOM;
48         return (x);
49       case INF:
50         errno = ERANGE;
51         return (ispos (x) ? z_infinity_f.f : -z_infinity_f.f);
52     }
53 
54   y = fabs (x);
55 
56   if (!cosineh && x < 0.0)
57     sgn = -1;
58 
59   if ((y > 1.0 && !cosineh) || cosineh)
60     {
61       if (y > BIGX)
62         {
63           w = y - LNV;
64 
65           /* Check for w > maximum here. */
66           if (w > BIGX)
67             {
68               errno = ERANGE;
69               return (x);
70             }
71 
72           z = exp (w);
73 
74           if (w > WBAR)
75             res = z * (V_OVER2_MINUS1 + 1.0);
76         }
77 
78       else
79         {
80           z = exp (y);
81           if (cosineh)
82             res = (z + 1 / z) / 2.0;
83           else
84             res = (z - 1 / z) / 2.0;
85         }
86 
87       if (sgn < 0)
88         res = -res;
89     }
90   else
91     {
92       /* Check for y being too small. */
93       if (y < z_rooteps_f)
94         {
95           res = x;
96         }
97       /* Calculate the Taylor series. */
98       else
99         {
100           f = x * x;
101           Q = f + q[0];
102           P = p[1] * f + p[0];
103           R = f * (P / Q);
104 
105           res = x + x * R;
106         }
107     }
108 
109   return (res);
110 }
111