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 csinh = __csinh 31 32 /* INDENT OFF */ 33 /* 34 * dcomplex csinh(dcomplex z); 35 * 36 * z -z x -x 37 * e - e e (cos(y)+i*sin(y)) - e (cos(-y)+i*sin(-y)) 38 * sinh z = -------------- = --------------------------------------------- 39 * 2 2 40 * x -x x -x 41 * cos(y) ( e - e ) + i*sin(y) (e + e ) 42 * = -------------------------------------------- 43 * 2 44 * 45 * = cos(y) sinh(x) + i sin(y) cosh(x) 46 * 47 * Implementation Note 48 * ------------------- 49 * 50 * |x| -|x| |x| -2|x| -2|x| -P-4 51 * Note that e +- e = e ( 1 +- e ). If e < 2 , where 52 * 53 * P stands for the number of significant bits of the machine precision, 54 * |x| 55 * then the result will be rounded to e . Therefore, we have 56 * 57 * z 58 * e 59 * sinh z = ----- if |x| >= (P/2 + 2)*ln2 60 * 2 61 * 62 * EXCEPTION (conform to ISO/IEC 9899:1999(E)): 63 * csinh(0,0)=(0,0) 64 * csinh(0,inf)=(+-0,NaN) 65 * csinh(0,NaN)=(+-0,NaN) 66 * csinh(x,inf) = (NaN,NaN) for finite positive x 67 * csinh(x,NaN) = (NaN,NaN) for finite non-zero x 68 * csinh(inf,0) = (inf, 0) 69 * csinh(inf,y) = (inf*cos(y),inf*sin(y)) for positive finite y 70 * csinh(inf,inf) = (+-inf,NaN) 71 * csinh(inf,NaN) = (+-inf,NaN) 72 * csinh(NaN,0) = (NaN,0) 73 * csinh(NaN,y) = (NaN,NaN) for non-zero y 74 * csinh(NaN,NaN) = (NaN,NaN) 75 */ 76 /* INDENT ON */ 77 78 #include "libm.h" /* cosh/exp/fabs/scalbn/sinh/sincos/__k_cexp */ 79 #include "complex_wrapper.h" 80 81 dcomplex 82 csinh(dcomplex z) { 83 double t, x, y, S, C; 84 int hx, ix, lx, hy, iy, ly, n; 85 dcomplex ans; 86 87 x = D_RE(z); 88 y = D_IM(z); 89 hx = HI_WORD(x); 90 lx = LO_WORD(x); 91 ix = hx & 0x7fffffff; 92 hy = HI_WORD(y); 93 ly = LO_WORD(y); 94 iy = hy & 0x7fffffff; 95 x = fabs(x); 96 y = fabs(y); 97 98 (void) sincos(y, &S, &C); 99 if (ix >= 0x403c0000) { /* |x| > 28 = prec/2 (14,28,34,60) */ 100 if (ix >= 0x40862E42) { /* |x| > 709.78... ~ log(2**1024) */ 101 if (ix >= 0x7ff00000) { /* |x| is inf or NaN */ 102 if ((iy | ly) == 0) { 103 D_RE(ans) = x; 104 D_IM(ans) = y; 105 } else if (iy >= 0x7ff00000) { 106 D_RE(ans) = x; 107 D_IM(ans) = x - y; 108 } else { 109 D_RE(ans) = C * x; 110 D_IM(ans) = S * x; 111 } 112 } else { 113 /* return exp(x)=t*2**n */ 114 t = __k_cexp(x, &n); 115 D_RE(ans) = scalbn(C * t, n - 1); 116 D_IM(ans) = scalbn(S * t, n - 1); 117 } 118 } else { 119 t = exp(x) * 0.5; 120 D_RE(ans) = C * t; 121 D_IM(ans) = S * t; 122 } 123 } else { 124 if ((ix | lx) == 0) { /* x = 0, return (0,S) */ 125 D_RE(ans) = 0.0; 126 D_IM(ans) = S; 127 } else { 128 D_RE(ans) = C * sinh(x); 129 D_IM(ans) = S * cosh(x); 130 } 131 } 132 if (hx < 0) 133 D_RE(ans) = -D_RE(ans); 134 if (hy < 0) 135 D_IM(ans) = -D_IM(ans); 136 return (ans); 137 } 138