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 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 23 */ 24 /* 25 * Copyright 2005 Sun Microsystems, Inc. All rights reserved. 26 * Use is subject to license terms. 27 */ 28 29 #pragma weak __cosh = cosh 30 31 /* INDENT OFF */ 32 /* 33 * cosh(x) 34 * Code originated from 4.3bsd. 35 * Modified by K.C. Ng for SUN 4.0 libm. 36 * Method : 37 * 1. Replace x by |x| (cosh(x) = cosh(-x)). 38 * 2. 39 * [ exp(x) - 1 ]^2 40 * 0 <= x <= 0.3465 : cosh(x) := 1 + ------------------- 41 * 2*exp(x) 42 * 43 * exp(x) + 1/exp(x) 44 * 0.3465 <= x <= 22 : cosh(x) := ------------------- 45 * 2 46 * 22 <= x <= lnovft : cosh(x) := exp(x)/2 47 * lnovft <= x < INF : cosh(x) := scalbn(exp(x-1024*ln2),1023) 48 * 49 * Note: .3465 is a number near one half of ln2. 50 * 51 * Special cases: 52 * cosh(x) is |x| if x is +INF, -INF, or NaN. 53 * only cosh(0)=1 is exact for finite x. 54 */ 55 /* INDENT ON */ 56 57 #include "libm.h" 58 59 static const double 60 ln2 = 6.93147180559945286227e-01, 61 ln2hi = 6.93147180369123816490e-01, 62 ln2lo = 1.90821492927058770002e-10, 63 lnovft = 7.09782712893383973096e+02; 64 65 double 66 cosh(double x) { 67 double t, w; 68 69 w = fabs(x); 70 if (!finite(w)) 71 return (w * w); 72 if (w < 0.3465) { 73 t = expm1(w); 74 w = 1.0 + t; 75 if (w != 1.0) 76 w = 1.0 + (t * t) / (w + w); 77 return (w); 78 } else if (w < 22.0) { 79 t = exp(w); 80 return (0.5 * (t + 1.0 / t)); 81 } else if (w <= lnovft) { 82 return (0.5 * exp(w)); 83 } else { 84 w = (w - 1024 * ln2hi) - 1024 * ln2lo; 85 if (w >= ln2) 86 return (_SVID_libm_err(x, x, 5)); 87 else 88 return (scalbn(exp(w), 1023)); 89 } 90 } 91