RTAI API: base/math/e_sinh.c Source File

00001 /* @(#)e_sinh.c 5.1 93/09/24 */ 00002 /* 00003 * ==================================================== 00004 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 00005 * 00006 * Developed at SunPro, a Sun Microsystems, Inc. business. 00007 * Permission to use, copy, modify, and distribute this 00008 * software is freely granted, provided that this notice 00009 * is preserved. 00010 * ==================================================== 00011 */ 00012 00013 #if defined(LIBM_SCCS) && !defined(lint) 00014 static char rcsid[] = "$NetBSD: e_sinh.c,v 1.7 1995/05/10 20:46:13 jtc Exp $"; 00015 #endif 00016 00017 /* __ieee754_sinh(x) 00018 * Method : 00019 * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2 00020 * 1. Replace x by |x| (sinh(-x) = -sinh(x)). 00021 * 2. 00022 * E + E/(E+1) 00023 * 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x) 00024 * 2 00025 * 00026 * 22 <= x <= lnovft : sinh(x) := exp(x)/2 00027 * lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2) 00028 * ln2ovft < x : sinh(x) := x*shuge (overflow) 00029 * 00030 * Special cases: 00031 * sinh(x) is |x| if x is +INF, -INF, or NaN. 00032 * only sinh(0)=0 is exact for finite x. 00033 */ 00034 00035 #include "math.h" 00036 #include "mathP.h" 00037 00038 #ifdef __STDC__ 00039 static const double one = 1.0, shuge = 1.0e307; 00040 #else 00041 static double one = 1.0, shuge = 1.0e307; 00042 #endif 00043 00044 #ifdef __STDC__ 00045 double __ieee754_sinh(double x) 00046 #else 00047 double __ieee754_sinh(x) 00048 double x; 00049 #endif 00050 { 00051 double t,w,h; 00052 int32_t ix,jx; 00053 u_int32_t lx; 00054 00055 /* High word of |x|. */ 00056 GET_HIGH_WORD(jx,x); 00057 ix = jx&0x7fffffff; 00058 00059 /* x is INF or NaN */ 00060 if(ix>=0x7ff00000) return x+x; 00061 00062 h = 0.5; 00063 if (jx<0) h = -h; 00064 /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */ 00065 if (ix < 0x40360000) { /* |x|<22 */ 00066 if (ix<0x3e300000) /* |x|<2**-28 */ 00067 if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */ 00068 t = expm1(fabs(x)); 00069 if(ix<0x3ff00000) return h*(2.0*t-t*t/(t+one)); 00070 return h*(t+t/(t+one)); 00071 } 00072 00073 /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */ 00074 if (ix < 0x40862E42) return h*__ieee754_exp(fabs(x)); 00075 00076 /* |x| in [log(maxdouble), overflowthresold] */ 00077 GET_LOW_WORD(lx,x); 00078 if (ix<0x408633CE || ((ix==0x408633ce)&&(lx<=(u_int32_t)0x8fb9f87d))) { 00079 w = __ieee754_exp(0.5*fabs(x)); 00080 t = h*w; 00081 return t*w; 00082 } 00083 00084 /* |x| > overflowthresold, sinh(x) overflow */ 00085 return x*shuge; 00086 }