1 /* 2 * ==================================================== 3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 4 * 5 * Developed at SunPro, a Sun Microsystems, Inc. business. 6 * Permission to use, copy, modify, and distribute this 7 * software is freely granted, provided that this notice 8 * is preserved. 9 * ==================================================== 10 */ 11 12 /* 13 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> 14 * 15 * Permission to use, copy, modify, and distribute this software for any 16 * purpose with or without fee is hereby granted, provided that the above 17 * copyright notice and this permission notice appear in all copies. 18 * 19 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES 20 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF 21 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR 22 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 23 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN 24 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 25 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 26 */ 27 28 /* coshl(x) 29 * Method : 30 * mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2 31 * 1. Replace x by |x| (coshl(x) = coshl(-x)). 32 * 2. 33 * [ exp(x) - 1 ]^2 34 * 0 <= x <= ln2/2 : coshl(x) := 1 + ------------------- 35 * 2*exp(x) 36 * 37 * exp(x) + 1/exp(x) 38 * ln2/2 <= x <= 22 : coshl(x) := ------------------- 39 * 2 40 * 22 <= x <= lnovft : coshl(x) := expl(x)/2 41 * lnovft <= x <= ln2ovft: coshl(x) := expl(x/2)/2 * expl(x/2) 42 * ln2ovft < x : coshl(x) := huge*huge (overflow) 43 * 44 * Special cases: 45 * coshl(x) is |x| if x is +INF, -INF, or NaN. 46 * only coshl(0)=1 is exact for finite x. 47 */ 48 49 #include <math.h> 50 51 #include "math_private.h" 52 53 static const long double one = 1.0, half = 0.5, huge = 1.0e4900L, 54 ovf_thresh = 1.1357216553474703894801348310092223067821E4L; 55 56 long double 57 coshl(long double x) 58 { 59 long double t, w; 60 int32_t ex; 61 ieee_quad_shape_type u; 62 63 u.value = x; 64 ex = u.parts32.mswhi & 0x7fffffff; 65 66 /* Absolute value of x. */ 67 u.parts32.mswhi = ex; 68 69 /* x is INF or NaN */ 70 if (ex >= 0x7fff0000) 71 return x * x; 72 73 /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */ 74 if (ex < 0x3ffd62e4) /* 0.3465728759765625 */ 75 { 76 t = expm1l (u.value); 77 w = one + t; 78 if (ex < 0x3fb80000) /* |x| < 2^-116 */ 79 return w; /* cosh(tiny) = 1 */ 80 81 return one + (t * t) / (w + w); 82 } 83 84 /* |x| in [0.5*ln2,40], return (exp(|x|)+1/exp(|x|)/2; */ 85 if (ex < 0x40044000) 86 { 87 t = expl (u.value); 88 return half * t + half / t; 89 } 90 91 /* |x| in [22, ln(maxdouble)] return half*exp(|x|) */ 92 if (ex <= 0x400c62e3) /* 11356.375 */ 93 return half * expl (u.value); 94 95 /* |x| in [log(maxdouble), overflowthresold] */ 96 if (u.value <= ovf_thresh) 97 { 98 w = expl (half * u.value); 99 t = half * w; 100 return t * w; 101 } 102 103 /* |x| > overflowthresold, cosh(x) overflow */ 104 return huge * huge; 105 } 106