1 /* @(#)e_sinh.c 5.1 93/09/24 */ 2 /* 3 * ==================================================== 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5 * 6 * Developed at SunPro, a Sun Microsystems, Inc. business. 7 * Permission to use, copy, modify, and distribute this 8 * software is freely granted, provided that this notice 9 * is preserved. 10 * ==================================================== 11 */ 12 13 /* 14 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> 15 * 16 * Permission to use, copy, modify, and distribute this software for any 17 * purpose with or without fee is hereby granted, provided that the above 18 * copyright notice and this permission notice appear in all copies. 19 * 20 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES 21 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF 22 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR 23 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 24 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN 25 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 26 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 27 */ 28 29 /* sinhl(x) 30 * Method : 31 * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2 32 * 1. Replace x by |x| (sinhl(-x) = -sinhl(x)). 33 * 2. 34 * E + E/(E+1) 35 * 0 <= x <= 25 : sinhl(x) := --------------, E=expm1l(x) 36 * 2 37 * 38 * 25 <= x <= lnovft : sinhl(x) := expl(x)/2 39 * lnovft <= x <= ln2ovft: sinhl(x) := expl(x/2)/2 * expl(x/2) 40 * ln2ovft < x : sinhl(x) := x*shuge (overflow) 41 * 42 * Special cases: 43 * sinhl(x) is |x| if x is +INF, -INF, or NaN. 44 * only sinhl(0)=0 is exact for finite x. 45 */ 46 47 #include <openlibm_math.h> 48 49 #include "math_private.h" 50 51 static const long double one = 1.0, shuge = 1.0e4931L, 52 ovf_thresh = 1.1357216553474703894801348310092223067821E4L; 53 54 long double 55 sinhl(long double x) 56 { 57 long double t, w, h; 58 u_int32_t jx, ix; 59 ieee_quad_shape_type u; 60 61 /* Words of |x|. */ 62 u.value = x; 63 jx = u.parts32.mswhi; 64 ix = jx & 0x7fffffff; 65 66 /* x is INF or NaN */ 67 if (ix >= 0x7fff0000) 68 return x + x; 69 70 h = 0.5; 71 if (jx & 0x80000000) 72 h = -h; 73 74 /* Absolute value of x. */ 75 u.parts32.mswhi = ix; 76 77 /* |x| in [0,40], return sign(x)*0.5*(E+E/(E+1))) */ 78 if (ix <= 0x40044000) 79 { 80 if (ix < 0x3fc60000) /* |x| < 2^-57 */ 81 if (shuge + x > one) 82 return x; /* sinh(tiny) = tiny with inexact */ 83 t = expm1l (u.value); 84 if (ix < 0x3fff0000) 85 return h * (2.0 * t - t * t / (t + one)); 86 return h * (t + t / (t + one)); 87 } 88 89 /* |x| in [40, log(maxdouble)] return 0.5*exp(|x|) */ 90 if (ix <= 0x400c62e3) /* 11356.375 */ 91 return h * expl (u.value); 92 93 /* |x| in [log(maxdouble), overflowthreshold] 94 Overflow threshold is log(2 * maxdouble). */ 95 if (u.value <= ovf_thresh) 96 { 97 w = expl (0.5 * u.value); 98 t = h * w; 99 return t * w; 100 } 101 102 /* |x| > overflowthreshold, sinhl(x) overflow */ 103 return x * shuge; 104 } 105