e_cosh.c revision 260067
1152286Semax 2152286Semax/* @(#)e_cosh.c 1.3 95/01/18 */ 3152286Semax/* 4152286Semax * ==================================================== 5152286Semax * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 6152286Semax * 7152286Semax * Developed at SunSoft, a Sun Microsystems, Inc. business. 8152286Semax * Permission to use, copy, modify, and distribute this 9152286Semax * software is freely granted, provided that this notice 10152286Semax * is preserved. 11152286Semax * ==================================================== 12152286Semax */ 13152286Semax 14152286Semax#include <sys/cdefs.h> 15152286Semax__FBSDID("$FreeBSD: head/lib/msun/src/e_cosh.c 260067 2013-12-30 01:06:21Z kargl $"); 16152286Semax 17152286Semax/* __ieee754_cosh(x) 18152286Semax * Method : 19152286Semax * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2 20152286Semax * 1. Replace x by |x| (cosh(x) = cosh(-x)). 21152286Semax * 2. 22152286Semax * [ exp(x) - 1 ]^2 23152286Semax * 0 <= x <= ln2/2 : cosh(x) := 1 + ------------------- 24152286Semax * 2*exp(x) 25152286Semax * 26152286Semax * exp(x) + 1/exp(x) 27152286Semax * ln2/2 <= x <= 22 : cosh(x) := ------------------- 28152286Semax * 2 29152286Semax * 22 <= x <= lnovft : cosh(x) := exp(x)/2 30152286Semax * lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2) 31152286Semax * ln2ovft < x : cosh(x) := huge*huge (overflow) 32152286Semax * 33152286Semax * Special cases: 34152286Semax * cosh(x) is |x| if x is +INF, -INF, or NaN. 35152286Semax * only cosh(0)=1 is exact for finite x. 36152286Semax */ 37152286Semax 38152286Semax#include <float.h> 39165683Syar 40152286Semax#include "math.h" 41152286Semax#include "math_private.h" 42152286Semax 43152286Semaxstatic const double one = 1.0, half=0.5, huge = 1.0e300; 44152286Semax 45152286Semaxdouble 46152286Semax__ieee754_cosh(double x) 47165664Syar{ 48165664Syar double t,w; 49152286Semax int32_t ix; 50152286Semax 51152286Semax /* High word of |x|. */ 52152286Semax GET_HIGH_WORD(ix,x); 53152286Semax ix &= 0x7fffffff; 54152286Semax 55152286Semax /* x is INF or NaN */ 56152286Semax if(ix>=0x7ff00000) return x*x; 57152286Semax 58152286Semax /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ 59152286Semax if(ix<0x3fd62e43) { 60152286Semax t = expm1(fabs(x)); 61152286Semax w = one+t; 62152286Semax if (ix<0x3c800000) return w; /* cosh(tiny) = 1 */ 63152286Semax return one+(t*t)/(w+w); 64152286Semax } 65152286Semax 66152286Semax /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */ 67152286Semax if (ix < 0x40360000) { 68152286Semax t = __ieee754_exp(fabs(x)); 69152286Semax return half*t+half/t; 70152286Semax } 71152286Semax 72152286Semax /* |x| in [22, log(maxdouble)] return half*exp(|x|) */ 73152286Semax if (ix < 0x40862E42) return half*__ieee754_exp(fabs(x)); 74152286Semax 75152286Semax /* |x| in [log(maxdouble), overflowthresold] */ 76152286Semax if (ix<=0x408633CE) 77152286Semax return __ldexp_exp(fabs(x), -1); 78152286Semax 79152286Semax /* |x| > overflowthresold, cosh(x) overflow */ 80152286Semax return huge*huge; 81152286Semax} 82152286Semax 83152286Semax#if (LDBL_MANT_DIG == 53) 84152286Semax__weak_reference(cosh, coshl); 85152286Semax#endif 86152286Semax