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