e_sinh.c revision 2116
1196212Sscottl/* @(#)e_sinh.c 5.1 93/09/24 */ 2196212Sscottl/* 3196212Sscottl * ==================================================== 4196212Sscottl * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5196212Sscottl * 6196212Sscottl * Developed at SunPro, a Sun Microsystems, Inc. business. 7196212Sscottl * Permission to use, copy, modify, and distribute this 8196212Sscottl * software is freely granted, provided that this notice 9196212Sscottl * is preserved. 10196212Sscottl * ==================================================== 11196212Sscottl */ 12196212Sscottl 13196212Sscottl#ifndef lint 14196212Sscottlstatic char rcsid[] = "$Id: e_sinh.c,v 1.5 1994/08/18 23:06:03 jtc Exp $"; 15196212Sscottl#endif 16196212Sscottl 17196212Sscottl/* __ieee754_sinh(x) 18196212Sscottl * Method : 19196212Sscottl * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2 20196212Sscottl * 1. Replace x by |x| (sinh(-x) = -sinh(x)). 21196212Sscottl * 2. 22196212Sscottl * E + E/(E+1) 23196212Sscottl * 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x) 24196212Sscottl * 2 25196212Sscottl * 26196212Sscottl * 22 <= x <= lnovft : sinh(x) := exp(x)/2 27196212Sscottl * lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2) 28196212Sscottl * ln2ovft < x : sinh(x) := x*shuge (overflow) 29196212Sscottl * 30196212Sscottl * Special cases: 31196212Sscottl * sinh(x) is |x| if x is +INF, -INF, or NaN. 32196212Sscottl * only sinh(0)=0 is exact for finite x. 33196212Sscottl */ 34196212Sscottl 35196212Sscottl#include "math.h" 36196212Sscottl#include "math_private.h" 37196212Sscottl 38196212Sscottl#ifdef __STDC__ 39196212Sscottlstatic const double one = 1.0, shuge = 1.0e307; 40196212Sscottl#else 41196212Sscottlstatic double one = 1.0, shuge = 1.0e307; 42196212Sscottl#endif 43196212Sscottl 44196212Sscottl#ifdef __STDC__ 45196212Sscottl double __ieee754_sinh(double x) 46196212Sscottl#else 47196212Sscottl double __ieee754_sinh(x) 48196212Sscottl double x; 49196212Sscottl#endif 50196212Sscottl{ 51196212Sscottl double t,w,h; 52196212Sscottl int32_t ix,jx; 53196212Sscottl u_int32_t lx; 54196212Sscottl 55196212Sscottl /* High word of |x|. */ 56196212Sscottl GET_HIGH_WORD(jx,x); 57196212Sscottl ix = jx&0x7fffffff; 58196212Sscottl 59196212Sscottl /* x is INF or NaN */ 60196212Sscottl if(ix>=0x7ff00000) return x+x; 61196212Sscottl 62196212Sscottl h = 0.5; 63196212Sscottl if (jx<0) h = -h; 64196212Sscottl /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */ 65196212Sscottl if (ix < 0x40360000) { /* |x|<22 */ 66196212Sscottl if (ix<0x3e300000) /* |x|<2**-28 */ 67196212Sscottl if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */ 68196212Sscottl t = expm1(fabs(x)); 69196212Sscottl if(ix<0x3ff00000) return h*(2.0*t-t*t/(t+one)); 70196212Sscottl return h*(t+t/(t+one)); 71196212Sscottl } 72196212Sscottl 73196212Sscottl /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */ 74196212Sscottl if (ix < 0x40862E42) return h*__ieee754_exp(fabs(x)); 75196212Sscottl 76196212Sscottl /* |x| in [log(maxdouble), overflowthresold] */ 77196212Sscottl GET_LOW_WORD(lx,x); 78196212Sscottl if (ix<0x408633CE || (ix==0x408633ce)&&(lx<=(u_int32_t)0x8fb9f87d)) { 79196212Sscottl w = __ieee754_exp(0.5*fabs(x)); 80196212Sscottl t = h*w; 81196212Sscottl return t*w; 82196212Sscottl } 83196212Sscottl 84196212Sscottl /* |x| > overflowthresold, sinh(x) overflow */ 85196212Sscottl return x*shuge; 86196212Sscottl} 87196212Sscottl