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 __ccosh = ccosh
31
32/* INDENT OFF */
33/*
34 * dcomplex ccosh(dcomplex z);
35 *
36 *             z      -z       x                      -x
37 *            e   +  e        e  (cos(y)+i*sin(y)) + e  (cos(-y)+i*sin(-y))
38 * cosh 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) cosh(x)  + i sin(y) sinh(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 *     cosh z = -----  if |x| >= (P/2 + 2)*ln2
60 *                2
61 *
62 * EXCEPTION (conform to ISO/IEC 9899:1999(E)):
63 *      ccosh(0,0)=(1,0)
64 *      ccosh(0,inf)=(NaN,+-0)
65 *      ccosh(0,NaN)=(NaN,+-0)
66 *      ccosh(x,inf) = (NaN,NaN) for finite non-zero x
67 *      ccosh(x,NaN) = (NaN,NaN) for finite non-zero x
68 *      ccosh(inf,0) = (inf, 0)
69 *      ccosh(inf,y) = (inf*cos(y),inf*sin(y)) for finite non-zero y
70 *      ccosh(inf,inf) = (+-inf,NaN)
71 *      ccosh(inf,NaN) = (+inf,NaN)
72 *      ccosh(NaN,0) = (NaN,+-0)
73 *      ccosh(NaN,y) = (NaN,NaN) for non-zero y
74 *      ccosh(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
81dcomplex
82ccosh(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				t = __k_cexp(x, &n);
114						/* return exp(x)=t*2**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 (C,0) */
125			D_RE(ans) = C;
126			D_IM(ans) = 0.0;
127		} else {
128			D_RE(ans) = C * cosh(x);
129			D_IM(ans) = S * sinh(x);
130		}
131	}
132	if ((hx ^ hy) < 0)
133		D_IM(ans) = -D_IM(ans);
134	return (ans);
135}
136