1*1e651e1eSRoland Levillain 2*1e651e1eSRoland Levillain /* @(#)e_cosh.c 1.3 95/01/18 */ 3*1e651e1eSRoland Levillain /* 4*1e651e1eSRoland Levillain * ==================================================== 5*1e651e1eSRoland Levillain * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 6*1e651e1eSRoland Levillain * 7*1e651e1eSRoland Levillain * Developed at SunSoft, a Sun Microsystems, Inc. business. 8*1e651e1eSRoland Levillain * Permission to use, copy, modify, and distribute this 9*1e651e1eSRoland Levillain * software is freely granted, provided that this notice 10*1e651e1eSRoland Levillain * is preserved. 11*1e651e1eSRoland Levillain * ==================================================== 12*1e651e1eSRoland Levillain */ 13*1e651e1eSRoland Levillain 14*1e651e1eSRoland Levillain /* __ieee754_cosh(x) 15*1e651e1eSRoland Levillain * Method : 16*1e651e1eSRoland Levillain * mathematically ieee_cosh(x) if defined to be (ieee_exp(x)+ieee_exp(-x))/2 17*1e651e1eSRoland Levillain * 1. Replace x by |x| (ieee_cosh(x) = ieee_cosh(-x)). 18*1e651e1eSRoland Levillain * 2. 19*1e651e1eSRoland Levillain * [ ieee_exp(x) - 1 ]^2 20*1e651e1eSRoland Levillain * 0 <= x <= ln2/2 : ieee_cosh(x) := 1 + ------------------- 21*1e651e1eSRoland Levillain * 2*ieee_exp(x) 22*1e651e1eSRoland Levillain * 23*1e651e1eSRoland Levillain * ieee_exp(x) + 1/ieee_exp(x) 24*1e651e1eSRoland Levillain * ln2/2 <= x <= 22 : ieee_cosh(x) := ------------------- 25*1e651e1eSRoland Levillain * 2 26*1e651e1eSRoland Levillain * 22 <= x <= lnovft : ieee_cosh(x) := ieee_exp(x)/2 27*1e651e1eSRoland Levillain * lnovft <= x <= ln2ovft: ieee_cosh(x) := ieee_exp(x/2)/2 * ieee_exp(x/2) 28*1e651e1eSRoland Levillain * ln2ovft < x : ieee_cosh(x) := huge*huge (overflow) 29*1e651e1eSRoland Levillain * 30*1e651e1eSRoland Levillain * Special cases: 31*1e651e1eSRoland Levillain * cosh(x) is |x| if x is +INF, -INF, or NaN. 32*1e651e1eSRoland Levillain * only ieee_cosh(0)=1 is exact for finite x. 33*1e651e1eSRoland Levillain */ 34*1e651e1eSRoland Levillain 35*1e651e1eSRoland Levillain #include "fdlibm.h" 36*1e651e1eSRoland Levillain 37*1e651e1eSRoland Levillain #ifdef __STDC__ 38*1e651e1eSRoland Levillain static const double one = 1.0, half=0.5, huge = 1.0e300; 39*1e651e1eSRoland Levillain #else 40*1e651e1eSRoland Levillain static double one = 1.0, half=0.5, huge = 1.0e300; 41*1e651e1eSRoland Levillain #endif 42*1e651e1eSRoland Levillain 43*1e651e1eSRoland Levillain #ifdef __STDC__ __ieee754_cosh(double x)44*1e651e1eSRoland Levillain double __ieee754_cosh(double x) 45*1e651e1eSRoland Levillain #else 46*1e651e1eSRoland Levillain double __ieee754_cosh(x) 47*1e651e1eSRoland Levillain double x; 48*1e651e1eSRoland Levillain #endif 49*1e651e1eSRoland Levillain { 50*1e651e1eSRoland Levillain double t,w; 51*1e651e1eSRoland Levillain int ix; 52*1e651e1eSRoland Levillain unsigned lx; 53*1e651e1eSRoland Levillain 54*1e651e1eSRoland Levillain /* High word of |x|. */ 55*1e651e1eSRoland Levillain ix = __HI(x); 56*1e651e1eSRoland Levillain ix &= 0x7fffffff; 57*1e651e1eSRoland Levillain 58*1e651e1eSRoland Levillain /* x is INF or NaN */ 59*1e651e1eSRoland Levillain if(ix>=0x7ff00000) return x*x; 60*1e651e1eSRoland Levillain 61*1e651e1eSRoland Levillain /* |x| in [0,0.5*ln2], return 1+ieee_expm1(|x|)^2/(2*ieee_exp(|x|)) */ 62*1e651e1eSRoland Levillain if(ix<0x3fd62e43) { 63*1e651e1eSRoland Levillain t = ieee_expm1(ieee_fabs(x)); 64*1e651e1eSRoland Levillain w = one+t; 65*1e651e1eSRoland Levillain if (ix<0x3c800000) return w; /* ieee_cosh(tiny) = 1 */ 66*1e651e1eSRoland Levillain return one+(t*t)/(w+w); 67*1e651e1eSRoland Levillain } 68*1e651e1eSRoland Levillain 69*1e651e1eSRoland Levillain /* |x| in [0.5*ln2,22], return (ieee_exp(|x|)+1/ieee_exp(|x|)/2; */ 70*1e651e1eSRoland Levillain if (ix < 0x40360000) { 71*1e651e1eSRoland Levillain t = __ieee754_exp(ieee_fabs(x)); 72*1e651e1eSRoland Levillain return half*t+half/t; 73*1e651e1eSRoland Levillain } 74*1e651e1eSRoland Levillain 75*1e651e1eSRoland Levillain /* |x| in [22, ieee_log(maxdouble)] return half*ieee_exp(|x|) */ 76*1e651e1eSRoland Levillain if (ix < 0x40862E42) return half*__ieee754_exp(ieee_fabs(x)); 77*1e651e1eSRoland Levillain 78*1e651e1eSRoland Levillain /* |x| in [log(maxdouble), overflowthresold] */ 79*1e651e1eSRoland Levillain lx = *( (((*(unsigned*)&one)>>29)) + (unsigned*)&x); 80*1e651e1eSRoland Levillain if (ix<0x408633CE || 81*1e651e1eSRoland Levillain (ix==0x408633ce)&&(lx<=(unsigned)0x8fb9f87d)) { 82*1e651e1eSRoland Levillain w = __ieee754_exp(half*ieee_fabs(x)); 83*1e651e1eSRoland Levillain t = half*w; 84*1e651e1eSRoland Levillain return t*w; 85*1e651e1eSRoland Levillain } 86*1e651e1eSRoland Levillain 87*1e651e1eSRoland Levillain /* |x| > overflowthresold, ieee_cosh(x) overflow */ 88*1e651e1eSRoland Levillain return huge*huge; 89*1e651e1eSRoland Levillain } 90