1207753Smm/* 2207753Smm * ==================================================== 3207753Smm * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 4207753Smm * 5207753Smm * Developed at SunPro, a Sun Microsystems, Inc. business. 6207753Smm * Permission to use, copy, modify, and distribute this 7207753Smm * software is freely granted, provided that this notice 8207753Smm * is preserved. 9207753Smm * ==================================================== 10207753Smm */ 11207753Smm 12207753Smm/* Changes for 128-bit long double are 13207753Smm Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov> 14207753Smm and are incorporated herein by permission of the author. The author 15207753Smm reserves the right to distribute this material elsewhere under different 16207753Smm copying permissions. These modifications are distributed here under 17207753Smm the following terms: 18207753Smm 19207753Smm This library is free software; you can redistribute it and/or 20207753Smm modify it under the terms of the GNU Lesser General Public 21207753Smm License as published by the Free Software Foundation; either 22207753Smm version 2.1 of the License, or (at your option) any later version. 23207753Smm 24207753Smm This library is distributed in the hope that it will be useful, 25207753Smm but WITHOUT ANY WARRANTY; without even the implied warranty of 26207753Smm MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 27207753Smm Lesser General Public License for more details. 28207753Smm 29207753Smm You should have received a copy of the GNU Lesser General Public 30207753Smm License along with this library; if not, see 31207753Smm <http://www.gnu.org/licenses/>. */ 32207753Smm 33207753Smm/* coshq(x) 34207753Smm * Method : 35207753Smm * mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2 36207753Smm * 1. Replace x by |x| (coshl(x) = coshl(-x)). 37207753Smm * 2. 38207753Smm * [ exp(x) - 1 ]^2 39207753Smm * 0 <= x <= ln2/2 : coshl(x) := 1 + ------------------- 40207753Smm * 2*exp(x) 41207753Smm * 42207753Smm * exp(x) + 1/exp(x) 43207753Smm * ln2/2 <= x <= 22 : coshl(x) := ------------------- 44207753Smm * 2 45207753Smm * 22 <= x <= lnovft : coshl(x) := expq(x)/2 46207753Smm * lnovft <= x <= ln2ovft: coshl(x) := expq(x/2)/2 * expq(x/2) 47207753Smm * ln2ovft < x : coshl(x) := huge*huge (overflow) 48207753Smm * 49207753Smm * Special cases: 50207753Smm * coshl(x) is |x| if x is +INF, -INF, or NaN. 51207753Smm * only coshl(0)=1 is exact for finite x. 52207753Smm */ 53207753Smm 54207753Smm#include "quadmath-imp.h" 55207753Smm 56207753Smmstatic const __float128 one = 1.0, half = 0.5, huge = 1.0e4900Q, 57207753Smmovf_thresh = 1.1357216553474703894801348310092223067821E4Q; 58207753Smm 59207753Smm__float128 60207753Smmcoshq (__float128 x) 61207753Smm{ 62207753Smm __float128 t, w; 63207753Smm int32_t ex; 64207753Smm ieee854_float128 u; 65207753Smm 66207753Smm u.value = x; 67207753Smm ex = u.words32.w0 & 0x7fffffff; 68207753Smm 69207753Smm /* Absolute value of x. */ 70223935Smm u.words32.w0 = ex; 71223935Smm 72223935Smm /* x is INF or NaN */ 73223935Smm if (ex >= 0x7fff0000) 74207753Smm return x * x; 75207753Smm 76207753Smm /* |x| in [0,0.5*ln2], return 1+expm1q(|x|)^2/(2*expq(|x|)) */ 77207753Smm if (ex < 0x3ffd62e4) /* 0.3465728759765625 */ 78207753Smm { 79207753Smm if (ex < 0x3fb80000) /* |x| < 2^-116 */ 80207753Smm return one; /* cosh(tiny) = 1 */ 81207753Smm t = expm1q (u.value); 82207753Smm w = one + t; 83207753Smm 84207753Smm return one + (t * t) / (w + w); 85207753Smm } 86207753Smm 87207753Smm /* |x| in [0.5*ln2,40], return (exp(|x|)+1/exp(|x|)/2; */ 88207753Smm if (ex < 0x40044000) 89207753Smm { 90207753Smm t = expq (u.value); 91207753Smm return half * t + half / t; 92207753Smm } 93207753Smm 94207753Smm /* |x| in [22, ln(maxdouble)] return half*exp(|x|) */ 95207753Smm if (ex <= 0x400c62e3) /* 11356.375 */ 96207753Smm return half * expq (u.value); 97207753Smm 98207753Smm /* |x| in [log(maxdouble), overflowthresold] */ 99207753Smm if (u.value <= ovf_thresh) 100207753Smm { 101207753Smm w = expq (half * u.value); 102207753Smm t = half * w; 103207753Smm return t * w; 104207753Smm } 105207753Smm 106207753Smm /* |x| > overflowthresold, cosh(x) overflow */ 107207753Smm return huge * huge; 108207753Smm} 109207753Smm