Back to index

glibc  2.9
e_coshl.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 long double one = 1.0L, half=0.5L, huge = 1.0e300L;
00043 #else
00044 static long double one = 1.0L, half=0.5L, huge = 1.0e300L;
00045 #endif
00046 
00047 #ifdef __STDC__
00048        long double __ieee754_coshl(long double x)
00049 #else
00050        long double __ieee754_coshl(x)
00051        long double x;
00052 #endif
00053 {
00054        long double t,w;
00055        int64_t ix;
00056 
00057     /* High word of |x|. */
00058        GET_LDOUBLE_MSW64(ix,x);
00059        ix &= 0x7fffffffffffffffLL;
00060 
00061     /* x is INF or NaN */
00062        if(ix>=0x7ff0000000000000LL) return x*x;
00063 
00064     /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
00065        if(ix<0x3fd62e42fefa39efLL) {
00066            t = __expm1l(fabsl(x));
00067            w = one+t;
00068            if (ix<0x3c80000000000000LL) return w;       /* cosh(tiny) = 1 */
00069            return one+(t*t)/(w+w);
00070        }
00071 
00072     /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
00073        if (ix < 0x4036000000000000LL) {
00074               t = __ieee754_expl(fabsl(x));
00075               return half*t+half/t;
00076        }
00077 
00078     /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
00079        if (ix < 0x40862e42fefa39efLL)  return half*__ieee754_expl(fabsl(x));
00080 
00081     /* |x| in [log(maxdouble), overflowthresold] */
00082         if (ix < 0x408633ce8fb9f87dLL) {
00083            w = __ieee754_expl(half*fabsl(x));
00084            t = half*w;
00085            return t*w;
00086        }
00087 
00088     /* |x| > overflowthresold, cosh(x) overflow */
00089        return huge*huge;
00090 }