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 __ccoshf = ccoshf
30
31#include "libm.h"
32#include "complex_wrapper.h"
33
34#if defined(__i386) && !defined(__amd64)
35extern int __swapRP(int);
36#endif
37
38static const float zero = 0.0F, half = 0.5F;
39
40fcomplex
41ccoshf(fcomplex z) {
42	float		t, x, y, S, C;
43	double		w;
44	int		hx, ix, hy, iy, n;
45	fcomplex	ans;
46
47	x = F_RE(z);
48	y = F_IM(z);
49	hx = THE_WORD(x);
50	ix = hx & 0x7fffffff;
51	hy = THE_WORD(y);
52	iy = hy & 0x7fffffff;
53	x = fabsf(x);
54	y = fabsf(y);
55
56	sincosf(y, &S, &C);
57	if (ix >= 0x41600000) {	/* |x| > 14 = prec/2 (14,28,34,60) */
58		if (ix >= 0x42B171AA) {	/* |x| > 88.722... ~ log(2**128) */
59			if (ix >= 0x7f800000) {	/* |x| is inf or NaN */
60				if (iy == 0) {
61					F_RE(ans) = x;
62					F_IM(ans) = y;
63				} else if (iy >= 0x7f800000) {
64					F_RE(ans) = x;
65					F_IM(ans) = x - y;
66				} else {
67					F_RE(ans) = C * x;
68					F_IM(ans) = S * x;
69				}
70			} else {
71#if defined(__i386) && !defined(__amd64)
72				int	rp = __swapRP(fp_extended);
73#endif
74				/* return (C, S) * exp(x) / 2 */
75				w = __k_cexp((double)x, &n);
76				F_RE(ans) = (float)scalbn(C * w, n - 1);
77				F_IM(ans) = (float)scalbn(S * w, n - 1);
78#if defined(__i386) && !defined(__amd64)
79				if (rp != fp_extended)
80					(void) __swapRP(rp);
81#endif
82			}
83		} else {
84			t = expf(x) * half;
85			F_RE(ans) = C * t;
86			F_IM(ans) = S * t;
87		}
88	} else {
89		if (ix == 0) {	/* x = 0, return (C,0) */
90			F_RE(ans) = C;
91			F_IM(ans) = zero;
92		} else {
93			F_RE(ans) = C * coshf(x);
94			F_IM(ans) = S * sinhf(x);
95		}
96	}
97	if ((hx ^ hy) < 0)
98		F_IM(ans) = -F_IM(ans);
99	return (ans);
100}
101