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

00001 /* @(#)e_cosh.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_cosh.c,v 1.7 1995/05/10 20:44:58 jtc Exp $";
00015 #endif
00016 
00017 /* __ieee754_cosh(x)
00018  * Method : 
00019  * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
00020  *  1. Replace x by |x| (cosh(x) = cosh(-x)). 
00021  *  2. 
00022  *                                              [ exp(x) - 1 ]^2 
00023  *      0        <= x <= ln2/2  :  cosh(x) := 1 + -------------------
00024  *                                     2*exp(x)
00025  *
00026  *                                        exp(x) +  1/exp(x)
00027  *      ln2/2    <= x <= 22     :  cosh(x) := -------------------
00028  *                                    2
00029  *      22       <= x <= lnovft :  cosh(x) := exp(x)/2 
00030  *      lnovft   <= x <= ln2ovft:  cosh(x) := exp(x/2)/2 * exp(x/2)
00031  *      ln2ovft  <  x       :  cosh(x) := huge*huge (overflow)
00032  *
00033  * Special cases:
00034  *  cosh(x) is |x| if x is +INF, -INF, or NaN.
00035  *  only cosh(0)=1 is exact for finite x.
00036  */
00037 
00038 #include "math.h"
00039 #include "mathP.h"
00040 
00041 #ifdef __STDC__
00042 static const double one = 1.0, half=0.5, huge = 1.0e300;
00043 #else
00044 static double one = 1.0, half=0.5, huge = 1.0e300;
00045 #endif
00046 
00047 #ifdef __STDC__
00048     double __ieee754_cosh(double x)
00049 #else
00050     double __ieee754_cosh(x)
00051     double x;
00052 #endif
00053 {   
00054     double t,w;
00055     int32_t ix;
00056     u_int32_t lx;
00057 
00058     /* High word of |x|. */
00059     GET_HIGH_WORD(ix,x);
00060     ix &= 0x7fffffff;
00061 
00062     /* x is INF or NaN */
00063     if(ix>=0x7ff00000) return x*x;  
00064 
00065     /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
00066     if(ix<0x3fd62e43) {
00067         t = expm1(fabs(x));
00068         w = one+t;
00069         if (ix<0x3c800000) return w;    /* cosh(tiny) = 1 */
00070         return one+(t*t)/(w+w);
00071     }
00072 
00073     /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
00074     if (ix < 0x40360000) {
00075         t = __ieee754_exp(fabs(x));
00076         return half*t+half/t;
00077     }
00078 
00079     /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
00080     if (ix < 0x40862E42)  return half*__ieee754_exp(fabs(x));
00081 
00082     /* |x| in [log(maxdouble), overflowthresold] */
00083     GET_LOW_WORD(lx,x);
00084     if (ix<0x408633CE || 
00085           ((ix==0x408633ce)&&(lx<=(u_int32_t)0x8fb9f87d))) {
00086         w = __ieee754_exp(half*fabs(x));
00087         t = half*w;
00088         return t*w;
00089     }
00090 
00091     /* |x| > overflowthresold, cosh(x) overflow */
00092     return huge*huge;
00093 }