128 bit sinh implementation.
[kopensolaris-gnu/glibc.git] / sysdeps / ieee754 / ldbl-128 / e_sinhl.c
1 /* e_sinhl.c -- long double version of e_sinh.c.
2  * Conversion to long double by Ulrich Drepper,
3  * Cygnus Support, drepper@cygnus.com.
4  */
5
6 /*
7  * ====================================================
8  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
9  *
10  * Developed at SunPro, a Sun Microsystems, Inc. business.
11  * Permission to use, copy, modify, and distribute this
12  * software is freely granted, provided that this notice
13  * is preserved.
14  * ====================================================
15  */
16
17 /* Changes for 128-bit long double contributed by
18    Stephen L. Moshier <moshier@na-net.ornl.gov> */
19
20 /* __ieee754_sinhl(x)
21  * Method :
22  * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
23  *      1. Replace x by |x| (sinhl(-x) = -sinhl(x)).
24  *      2.
25  *                                                   E + E/(E+1)
26  *          0        <= x <= 25     :  sinhl(x) := --------------, E=expm1l(x)
27  *                                                       2
28  *
29  *          25       <= x <= lnovft :  sinhl(x) := expl(x)/2
30  *          lnovft   <= x <= ln2ovft:  sinhl(x) := expl(x/2)/2 * expl(x/2)
31  *          ln2ovft  <  x           :  sinhl(x) := x*shuge (overflow)
32  *
33  * Special cases:
34  *      sinhl(x) is |x| if x is +INF, -INF, or NaN.
35  *      only sinhl(0)=0 is exact for finite x.
36  */
37
38 #include "math.h"
39 #include "math_private.h"
40
41 #ifdef __STDC__
42 static const long double one = 1.0, shuge = 1.0e4931L,
43 ovf_thresh = 1.1357216553474703894801348310092223067821E4L;
44 #else
45 static long double one = 1.0, shuge = 1.0e4931L,
46 ovf_thresh = 1.1357216553474703894801348310092223067821E4L;
47 #endif
48
49 #ifdef __STDC__
50 long double
51 __ieee754_sinhl (long double x)
52 #else
53 long double
54 __ieee754_sinhl (x)
55      long double x;
56 #endif
57 {
58   long double t, w, h;
59   u_int32_t jx, ix;
60   ieee854_long_double_shape_type u;
61
62   /* Words of |x|. */
63   u.value = x;
64   jx = u.parts32.w0;
65   ix = jx & 0x7fffffff;
66
67   /* x is INF or NaN */
68   if (ix >= 0x7fff0000)
69     return x + x;
70
71   h = 0.5;
72   if (jx & 0x80000000)
73     h = -h;
74
75   /* Absolute value of x.  */
76   u.parts32.w0 = ix;
77
78   /* |x| in [0,40], return sign(x)*0.5*(E+E/(E+1))) */
79   if (ix <= 0x40044000)
80     {
81       if (ix < 0x3fc60000) /* |x| < 2^-57 */
82         if (shuge + x > one)
83           return x;             /* sinh(tiny) = tiny with inexact */
84       t = __expm1l (u.value);
85       if (ix < 0x3fff0000)
86         return h * (2.0 * t - t * t / (t + one));
87       return h * (t + t / (t + one));
88     }
89
90   /* |x| in [40, log(maxdouble)] return 0.5*exp(|x|) */
91   if (ix <= 0x400c62e3) /* 11356.375 */
92     return h * __ieee754_expl (u.value);
93
94   /* |x| in [log(maxdouble), overflowthreshold]
95      Overflow threshold is log(2 * maxdouble).  */
96   if (u.value <= ovf_thresh)
97     {
98       w = __ieee754_expl (0.5 * u.value);
99       t = h * w;
100       return t * w;
101     }
102
103   /* |x| > overflowthreshold, sinhl(x) overflow */
104   return x * shuge;
105 }