1
2/* @(#)e_sinh.c 1.3 95/01/18 */
3/*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
7 * Developed at SunSoft, a Sun Microsystems, Inc. business.
8 * Permission to use, copy, modify, and distribute this
9 * software is freely granted, provided that this notice
10 * is preserved.
11 * ====================================================
12 */
13
14/* __ieee754_sinh(x)
15 * Method :
16 * mathematically ieee_sinh(x) if defined to be (ieee_exp(x)-ieee_exp(-x))/2
17 *	1. Replace x by |x| (ieee_sinh(-x) = -ieee_sinh(x)).
18 *	2.
19 *		                                    E + E/(E+1)
20 *	    0        <= x <= 22     :  ieee_sinh(x) := --------------, E=ieee_expm1(x)
21 *			       			        2
22 *
23 *	    22       <= x <= lnovft :  ieee_sinh(x) := ieee_exp(x)/2
24 *	    lnovft   <= x <= ln2ovft:  ieee_sinh(x) := ieee_exp(x/2)/2 * ieee_exp(x/2)
25 *	    ln2ovft  <  x	    :  ieee_sinh(x) := x*shuge (overflow)
26 *
27 * Special cases:
28 *	sinh(x) is |x| if x is +INF, -INF, or NaN.
29 *	only ieee_sinh(0)=0 is exact for finite x.
30 */
31
32#include "fdlibm.h"
33
34#ifdef __STDC__
35static const double one = 1.0, shuge = 1.0e307;
36#else
37static double one = 1.0, shuge = 1.0e307;
38#endif
39
40#ifdef __STDC__
41	double __ieee754_sinh(double x)
42#else
43	double __ieee754_sinh(x)
44	double x;
45#endif
46{
47	double t,w,h;
48	int ix,jx;
49	unsigned lx;
50
51    /* High word of |x|. */
52	jx = __HI(x);
53	ix = jx&0x7fffffff;
54
55    /* x is INF or NaN */
56	if(ix>=0x7ff00000) return x+x;
57
58	h = 0.5;
59	if (jx<0) h = -h;
60    /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
61	if (ix < 0x40360000) {		/* |x|<22 */
62	    if (ix<0x3e300000) 		/* |x|<2**-28 */
63		if(shuge+x>one) return x;/* ieee_sinh(tiny) = tiny with inexact */
64	    t = ieee_expm1(ieee_fabs(x));
65	    if(ix<0x3ff00000) return h*(2.0*t-t*t/(t+one));
66	    return h*(t+t/(t+one));
67	}
68
69    /* |x| in [22, ieee_log(maxdouble)] return 0.5*ieee_exp(|x|) */
70	if (ix < 0x40862E42)  return h*__ieee754_exp(ieee_fabs(x));
71
72    /* |x| in [log(maxdouble), overflowthresold] */
73	lx = *( (((*(unsigned*)&one)>>29)) + (unsigned*)&x);
74	if (ix<0x408633CE || (ix==0x408633ce)&&(lx<=(unsigned)0x8fb9f87d)) {
75	    w = __ieee754_exp(0.5*ieee_fabs(x));
76	    t = h*w;
77	    return t*w;
78	}
79
80    /* |x| > overflowthresold, ieee_sinh(x) overflow */
81	return x*shuge;
82}
83