1 /* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 5 * Common Development and Distribution License (the "License"). 6 * You may not use this file except in compliance with the License. 7 * 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9 * or http://www.opensolaris.org/os/licensing. 10 * See the License for the specific language governing permissions 11 * and limitations under the License. 12 * 13 * When distributing Covered Code, include this CDDL HEADER in each 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15 * If applicable, add the following below this CDDL HEADER, with the 16 * fields enclosed by brackets "[]" replaced with your own identifying 17 * information: Portions Copyright [yyyy] [name of copyright owner] 18 * 19 * CDDL HEADER END 20 */ 21 22 /* 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 24 */ 25 /* 26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 27 * Use is subject to license terms. 28 */ 29 30 #if defined(ELFOBJ) 31 #pragma weak coshl = __coshl 32 #endif 33 34 #include "libm.h" 35 #include "longdouble.h" 36 37 /* 38 * COSH(X) 39 * RETURN THE HYPERBOLIC COSINE OF X 40 * 41 * Method : 42 * 1. Replace x by |x| (COSH(x) = COSH(-x)). 43 * 2. 44 * [ EXP(x) - 1 ]^2 45 * 0 <= x <= 0.3465 : COSH(x) := 1 + ------------------- 46 * 2*EXP(x) 47 * 48 * EXP(x) + 1/EXP(x) 49 * 0.3465 <= x <= thresh : COSH(x) := ------------------- 50 * 2 51 * thresh <= x <= lnovft : COSH(x) := EXP(x)/2 52 * lnovft <= x < INF : COSH(x) := SCALBN(EXP(x-MEP1*ln2),ME) 53 * 54 * 55 * here 56 * 0.3465 a number that is near one half of ln2. 57 * thresh a number such that 58 * EXP(thresh)+EXP(-thresh)=EXP(thresh) 59 * lnovft logarithm of the overflow threshold 60 * = MEP1*ln2 chopped to machine precision. 61 * ME maximum exponent 62 * MEP1 maximum exponent plus 1 63 * 64 * Special cases: 65 * COSH(x) is |x| if x is +INF, -INF, or NaN. 66 * only COSH(0)=1 is exact for finite x. 67 */ 68 69 static const long double C[] = { 70 0.5L, 71 1.0L, 72 0.3465L, 73 45.0L, 74 1.135652340629414394879149e+04L, 75 7.004447686242549087858985e-16L, 76 2.710505431213761085018632e-20L, /* 2^-65 */ 77 }; 78 79 #define half C[0] 80 #define one C[1] 81 #define thr1 C[2] 82 #define thr2 C[3] 83 #define lnovft C[4] 84 #define lnovlo C[5] 85 #define tinyl C[6] 86 87 long double 88 coshl(long double x) { 89 long double w, t; 90 91 w = fabsl(x); 92 if (!finitel(w)) 93 return (w + w); /* x is INF or NaN */ 94 if (w < thr1) { 95 if (w < tinyl) 96 return (one + w); /* inexact+directed rounding */ 97 t = expm1l(w); 98 w = one + t; 99 w = one + (t * t) / (w + w); 100 return (w); 101 } 102 if (w < thr2) { 103 t = expl(w); 104 return (half * (t + one / t)); 105 } 106 if (w <= lnovft) 107 return (half * expl(w)); 108 return (scalbnl(expl((w - lnovft) - lnovlo), 16383)); 109 } 110