e_cosh.c revision 141296
1141296Sdas 2141296Sdas/* @(#)e_cosh.c 1.3 95/01/18 */ 32116Sjkh/* 42116Sjkh * ==================================================== 52116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 62116Sjkh * 7141296Sdas * Developed at SunSoft, a Sun Microsystems, Inc. business. 82116Sjkh * Permission to use, copy, modify, and distribute this 9141296Sdas * software is freely granted, provided that this notice 102116Sjkh * is preserved. 112116Sjkh * ==================================================== 122116Sjkh */ 132116Sjkh 142116Sjkh#ifndef lint 1550476Speterstatic char rcsid[] = "$FreeBSD: head/lib/msun/src/e_cosh.c 141296 2005-02-04 18:26:06Z das $"; 162116Sjkh#endif 172116Sjkh 182116Sjkh/* __ieee754_cosh(x) 19141296Sdas * Method : 202116Sjkh * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2 21141296Sdas * 1. Replace x by |x| (cosh(x) = cosh(-x)). 22141296Sdas * 2. 23141296Sdas * [ exp(x) - 1 ]^2 242116Sjkh * 0 <= x <= ln2/2 : cosh(x) := 1 + ------------------- 252116Sjkh * 2*exp(x) 262116Sjkh * 272116Sjkh * exp(x) + 1/exp(x) 282116Sjkh * ln2/2 <= x <= 22 : cosh(x) := ------------------- 292116Sjkh * 2 30141296Sdas * 22 <= x <= lnovft : cosh(x) := exp(x)/2 312116Sjkh * lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2) 322116Sjkh * ln2ovft < x : cosh(x) := huge*huge (overflow) 332116Sjkh * 342116Sjkh * Special cases: 352116Sjkh * cosh(x) is |x| if x is +INF, -INF, or NaN. 362116Sjkh * only cosh(0)=1 is exact for finite x. 372116Sjkh */ 382116Sjkh 392116Sjkh#include "math.h" 402116Sjkh#include "math_private.h" 412116Sjkh 422116Sjkhstatic const double one = 1.0, half=0.5, huge = 1.0e300; 432116Sjkh 4497407Salfreddouble 4597407Salfred__ieee754_cosh(double x) 468870Srgrimes{ 472116Sjkh double t,w; 482116Sjkh int32_t ix; 492116Sjkh u_int32_t lx; 502116Sjkh 512116Sjkh /* High word of |x|. */ 522116Sjkh GET_HIGH_WORD(ix,x); 532116Sjkh ix &= 0x7fffffff; 542116Sjkh 552116Sjkh /* x is INF or NaN */ 56141296Sdas if(ix>=0x7ff00000) return x*x; 572116Sjkh 582116Sjkh /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ 592116Sjkh if(ix<0x3fd62e43) { 602116Sjkh t = expm1(fabs(x)); 612116Sjkh w = one+t; 622116Sjkh if (ix<0x3c800000) return w; /* cosh(tiny) = 1 */ 632116Sjkh return one+(t*t)/(w+w); 642116Sjkh } 652116Sjkh 662116Sjkh /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */ 672116Sjkh if (ix < 0x40360000) { 682116Sjkh t = __ieee754_exp(fabs(x)); 692116Sjkh return half*t+half/t; 702116Sjkh } 712116Sjkh 722116Sjkh /* |x| in [22, log(maxdouble)] return half*exp(|x|) */ 732116Sjkh if (ix < 0x40862E42) return half*__ieee754_exp(fabs(x)); 742116Sjkh 752116Sjkh /* |x| in [log(maxdouble), overflowthresold] */ 762116Sjkh GET_LOW_WORD(lx,x); 778870Srgrimes if (ix<0x408633CE || 7817141Sjkh ((ix==0x408633ce)&&(lx<=(u_int32_t)0x8fb9f87d))) { 792116Sjkh w = __ieee754_exp(half*fabs(x)); 802116Sjkh t = half*w; 812116Sjkh return t*w; 822116Sjkh } 832116Sjkh 842116Sjkh /* |x| > overflowthresold, cosh(x) overflow */ 852116Sjkh return huge*huge; 862116Sjkh} 87