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 #pragma weak __coshl = coshl 31 32 #include "libm.h" 33 #include "longdouble.h" 34 35 /* 36 * COSH(X) 37 * RETURN THE HYPERBOLIC COSINE OF X 38 * 39 * Method : 40 * 1. Replace x by |x| (COSH(x) = COSH(-x)). 41 * 2. 42 * [ EXP(x) - 1 ]^2 43 * 0 <= x <= 0.3465 : COSH(x) := 1 + ------------------- 44 * 2*EXP(x) 45 * 46 * EXP(x) + 1/EXP(x) 47 * 0.3465 <= x <= thresh : COSH(x) := ------------------- 48 * 2 49 * thresh <= x <= lnovft : COSH(x) := EXP(x)/2 50 * lnovft <= x < INF : COSH(x) := SCALBN(EXP(x-MEP1*ln2),ME) 51 * 52 * 53 * here 54 * 0.3465 a number that is near one half of ln2. 55 * thresh a number such that 56 * EXP(thresh)+EXP(-thresh)=EXP(thresh) 57 * lnovft logarithm of the overflow threshold 58 * = MEP1*ln2 chopped to machine precision. 59 * ME maximum exponent 60 * MEP1 maximum exponent plus 1 61 * 62 * Special cases: 63 * COSH(x) is |x| if x is +INF, -INF, or NaN. 64 * only COSH(0)=1 is exact for finite x. 65 */ 66 67 static const long double C[] = { 68 0.5L, 69 1.0L, 70 0.3465L, 71 45.0L, 72 1.135652340629414394879149e+04L, 73 7.004447686242549087858985e-16L, 74 2.710505431213761085018632e-20L, /* 2^-65 */ 75 }; 76 77 #define half C[0] 78 #define one C[1] 79 #define thr1 C[2] 80 #define thr2 C[3] 81 #define lnovft C[4] 82 #define lnovlo C[5] 83 #define tinyl C[6] 84 85 long double 86 coshl(long double x) { 87 long double w, t; 88 89 w = fabsl(x); 90 if (!finitel(w)) 91 return (w + w); /* x is INF or NaN */ 92 if (w < thr1) { 93 if (w < tinyl) 94 return (one + w); /* inexact+directed rounding */ 95 t = expm1l(w); 96 w = one + t; 97 w = one + (t * t) / (w + w); 98 return (w); 99 } 100 if (w < thr2) { 101 t = expl(w); 102 return (half * (t + one / t)); 103 } 104 if (w <= lnovft) 105 return (half * expl(w)); 106 return (scalbnl(expl((w - lnovft) - lnovlo), 16383)); 107 } 108