• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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__
35 static const double one = 1.0, shuge = 1.0e307;
36 #else
37 static double one = 1.0, shuge = 1.0e307;
38 #endif
39 
40 #ifdef __STDC__
__ieee754_sinh(double x)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