149393c00Smartynas /*
249393c00Smartynas * ====================================================
349393c00Smartynas * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
449393c00Smartynas *
549393c00Smartynas * Developed at SunPro, a Sun Microsystems, Inc. business.
649393c00Smartynas * Permission to use, copy, modify, and distribute this
749393c00Smartynas * software is freely granted, provided that this notice
849393c00Smartynas * is preserved.
949393c00Smartynas * ====================================================
1049393c00Smartynas */
1149393c00Smartynas
1249393c00Smartynas /*
1349393c00Smartynas * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
1449393c00Smartynas *
1549393c00Smartynas * Permission to use, copy, modify, and distribute this software for any
1649393c00Smartynas * purpose with or without fee is hereby granted, provided that the above
1749393c00Smartynas * copyright notice and this permission notice appear in all copies.
1849393c00Smartynas *
1949393c00Smartynas * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
2049393c00Smartynas * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
2149393c00Smartynas * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
2249393c00Smartynas * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
2349393c00Smartynas * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
2449393c00Smartynas * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
2549393c00Smartynas * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
2649393c00Smartynas */
2749393c00Smartynas
2849393c00Smartynas /* coshl(x)
2949393c00Smartynas * Method :
3049393c00Smartynas * mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2
3149393c00Smartynas * 1. Replace x by |x| (coshl(x) = coshl(-x)).
3249393c00Smartynas * 2.
3349393c00Smartynas * [ exp(x) - 1 ]^2
3449393c00Smartynas * 0 <= x <= ln2/2 : coshl(x) := 1 + -------------------
3549393c00Smartynas * 2*exp(x)
3649393c00Smartynas *
3749393c00Smartynas * exp(x) + 1/exp(x)
3849393c00Smartynas * ln2/2 <= x <= 22 : coshl(x) := -------------------
3949393c00Smartynas * 2
4049393c00Smartynas * 22 <= x <= lnovft : coshl(x) := expl(x)/2
4149393c00Smartynas * lnovft <= x <= ln2ovft: coshl(x) := expl(x/2)/2 * expl(x/2)
4249393c00Smartynas * ln2ovft < x : coshl(x) := huge*huge (overflow)
4349393c00Smartynas *
4449393c00Smartynas * Special cases:
4549393c00Smartynas * coshl(x) is |x| if x is +INF, -INF, or NaN.
4649393c00Smartynas * only coshl(0)=1 is exact for finite x.
4749393c00Smartynas */
4849393c00Smartynas
4949393c00Smartynas #include <math.h>
5049393c00Smartynas
5149393c00Smartynas #include "math_private.h"
5249393c00Smartynas
5349393c00Smartynas static const long double one = 1.0, half = 0.5, huge = 1.0e4900L,
5449393c00Smartynas ovf_thresh = 1.1357216553474703894801348310092223067821E4L;
5549393c00Smartynas
5649393c00Smartynas long double
coshl(long double x)5749393c00Smartynas coshl(long double x)
5849393c00Smartynas {
5949393c00Smartynas long double t, w;
6049393c00Smartynas int32_t ex;
6149393c00Smartynas ieee_quad_shape_type u;
6249393c00Smartynas
6349393c00Smartynas u.value = x;
6449393c00Smartynas ex = u.parts32.mswhi & 0x7fffffff;
6549393c00Smartynas
6649393c00Smartynas /* Absolute value of x. */
6749393c00Smartynas u.parts32.mswhi = ex;
6849393c00Smartynas
6949393c00Smartynas /* x is INF or NaN */
7049393c00Smartynas if (ex >= 0x7fff0000)
7149393c00Smartynas return x * x;
7249393c00Smartynas
7349393c00Smartynas /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
7449393c00Smartynas if (ex < 0x3ffd62e4) /* 0.3465728759765625 */
7549393c00Smartynas {
7649393c00Smartynas t = expm1l (u.value);
7749393c00Smartynas w = one + t;
7849393c00Smartynas if (ex < 0x3fb80000) /* |x| < 2^-116 */
7949393c00Smartynas return w; /* cosh(tiny) = 1 */
8049393c00Smartynas
8149393c00Smartynas return one + (t * t) / (w + w);
8249393c00Smartynas }
8349393c00Smartynas
8449393c00Smartynas /* |x| in [0.5*ln2,40], return (exp(|x|)+1/exp(|x|)/2; */
8549393c00Smartynas if (ex < 0x40044000)
8649393c00Smartynas {
8749393c00Smartynas t = expl (u.value);
8849393c00Smartynas return half * t + half / t;
8949393c00Smartynas }
9049393c00Smartynas
9149393c00Smartynas /* |x| in [22, ln(maxdouble)] return half*exp(|x|) */
9249393c00Smartynas if (ex <= 0x400c62e3) /* 11356.375 */
9349393c00Smartynas return half * expl (u.value);
9449393c00Smartynas
9549393c00Smartynas /* |x| in [log(maxdouble), overflowthresold] */
9649393c00Smartynas if (u.value <= ovf_thresh)
9749393c00Smartynas {
9849393c00Smartynas w = expl (half * u.value);
9949393c00Smartynas t = half * w;
10049393c00Smartynas return t * w;
10149393c00Smartynas }
10249393c00Smartynas
10349393c00Smartynas /* |x| > overflowthresold, cosh(x) overflow */
10449393c00Smartynas return huge * huge;
10549393c00Smartynas }
106*2f2c0062Sguenther DEF_STD(coshl);
107