Back to index

glibc  2.9
e_cosh.c
Go to the documentation of this 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 "math_private.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 || ((ix==0x408633ce)&&(lx<=(u_int32_t)0x8fb9f87d))) {
00085            w = __ieee754_exp(half*fabs(x));
00086            t = half*w;
00087            return t*w;
00088        }
00089 
00090     /* |x| > overflowthresold, cosh(x) overflow */
00091        return huge*huge;
00092 }