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 __casinl = casinl
31
32#include "libm.h"		/* asinl/atanl/fabsl/isinfl/log1pl/logl/sqrtl */
33#include "complex_wrapper.h"
34#include "longdouble.h"
35
36/* INDENT OFF */
37static const long double
38zero = 0.0L,
39one = 1.0L,
40Acrossover = 1.5L,
41Bcrossover = 0.6417L,
42half = 0.5L,
43ln2 = 6.931471805599453094172321214581765680755e-0001L,
44Foursqrtu = 7.3344154702193886624856495681939326638255e-2466L,	/* 2**-8189 */
45#if defined(__x86)
46E = 5.4210108624275221700372640043497085571289e-20L,	/* 2**-64 */
47pi_4 = 0.7853981633974483095739921312272713294078130L,
48pi_4_l = 4.1668714592604391641479322342670193036704898e-20L,
49pi_2 = 1.5707963267948966191479842624545426588156260L,
50pi_2_l = 8.3337429185208783282958644685340386073409796e-20L;
51
52#else
53E = 9.6296497219361792652798897129246365926905e-35L,	/* 2**-113 */
54pi_4 = 0.7853981633974483096156608458198756993697670L,
55pi_4_l = 2.1679525325309452561992610065108379921905808e-35L,
56pi_2 = 1.5707963267948966192313216916397513987395340L,
57pi_2_l = 4.3359050650618905123985220130216759843811616e-35L;
58
59#endif
60/* INDENT ON */
61
62#if defined(__x86)
63static const int ip1 = 0x40400000;	/* 2**65 */
64#else
65static const int ip1 = 0x40710000;	/* 2**114 */
66#endif
67
68ldcomplex
69casinl(ldcomplex z) {
70	long double x, y, t, R, S, A, Am1, B, y2, xm1, xp1, Apx;
71	int ix, iy, hx, hy;
72	ldcomplex ans;
73
74	x = LD_RE(z);
75	y = LD_IM(z);
76	hx = HI_XWORD(x);
77	hy = HI_XWORD(y);
78	ix = hx & 0x7fffffff;
79	iy = hy & 0x7fffffff;
80	x = fabsl(x);
81	y = fabsl(y);
82
83	/* special cases */
84
85	/* x is inf or NaN */
86	if (ix >= 0x7fff0000) {	/* x is inf or NaN */
87		if (isinfl(x)) {	/* x is INF */
88			LD_IM(ans) = x;
89			if (iy >= 0x7fff0000) {
90				if (isinfl(y))
91					/* casin(inf + i inf) = pi/4 + i inf */
92					LD_RE(ans) = pi_4 + pi_4_l;
93				else	/* casin(inf + i NaN) = NaN + i inf */
94					LD_RE(ans) = y + y;
95			} else	/* casin(inf + iy) = pi/2 + i inf */
96				LD_RE(ans) = pi_2 + pi_2_l;
97		} else {		/* x is NaN */
98			if (iy >= 0x7fff0000) {
99				/* INDENT OFF */
100				/*
101				 * casin(NaN + i inf) = NaN  + i inf
102				 * casin(NaN + i NaN) = NaN  + i NaN
103				 */
104				/* INDENT ON */
105				LD_IM(ans) = y + y;
106				LD_RE(ans) = x + x;
107			} else {
108				/* INDENT OFF */
109				/* casin(NaN + i y ) = NaN  + i NaN */
110				/* INDENT ON */
111				LD_IM(ans) = LD_RE(ans) = x + y;
112			}
113		}
114		if (hx < 0)
115			LD_RE(ans) = -LD_RE(ans);
116		if (hy < 0)
117			LD_IM(ans) = -LD_IM(ans);
118		return (ans);
119	}
120
121	/* casin(+0 + i 0) = 0 + i 0. */
122	if (x == zero && y == zero)
123		return (z);
124
125	if (iy >= 0x7fff0000) {	/* y is inf or NaN */
126		if (isinfl(y)) {	/* casin(x + i inf) = 0 + i inf */
127			LD_IM(ans) = y;
128			LD_RE(ans) = zero;
129		} else {		/* casin(x + i NaN) = NaN + i NaN */
130			LD_IM(ans) = x + y;
131			if (x == zero)
132				LD_RE(ans) = x;
133			else
134				LD_RE(ans) = y;
135		}
136		if (hx < 0)
137			LD_RE(ans) = -LD_RE(ans);
138		if (hy < 0)
139			LD_IM(ans) = -LD_IM(ans);
140		return (ans);
141	}
142
143	if (y == zero) {	/* region 1: y=0 */
144		if (ix < 0x3fff0000) {	/* |x| < 1 */
145			LD_RE(ans) = asinl(x);
146			LD_IM(ans) = zero;
147		} else {
148			LD_RE(ans) = pi_2 + pi_2_l;
149			if (ix >= ip1)	/* |x| >= i386 ? 2**65 : 2**114 */
150				LD_IM(ans) = ln2 + logl(x);
151			else if (ix >= 0x3fff8000)	/* x > Acrossover */
152				LD_IM(ans) = logl(x + sqrtl((x - one) * (x +
153					one)));
154			else {
155				xm1 = x - one;
156				LD_IM(ans) = log1pl(xm1 + sqrtl(xm1 * (x +
157					one)));
158			}
159		}
160	} else if (y <= E * fabsl(x - one)) {	/* region 2: y < tiny*|x-1| */
161		if (ix < 0x3fff0000) {	/* x < 1 */
162			LD_RE(ans) = asinl(x);
163			LD_IM(ans) = y / sqrtl((one + x) * (one - x));
164		} else {
165			LD_RE(ans) = pi_2 + pi_2_l;
166			if (ix >= ip1)	/* i386 ? 2**65 : 2**114 */
167				LD_IM(ans) = ln2 + logl(x);
168			else if (ix >= 0x3fff8000)	/* x > Acrossover */
169				LD_IM(ans) = logl(x + sqrtl((x - one) * (x +
170					one)));
171			else
172				LD_IM(ans) = log1pl((x - one) + sqrtl((x -
173					one) * (x + one)));
174		}
175	} else if (y < Foursqrtu) {	/* region 3 */
176		t = sqrtl(y);
177		LD_RE(ans) = pi_2 - (t - pi_2_l);
178		LD_IM(ans) = t;
179	} else if (E * y - one >= x) {	/* region 4 */
180		LD_RE(ans) = x / y;	/* need to fix underflow cases */
181		LD_IM(ans) = ln2 + logl(y);
182	} else if (ix >= 0x5ffb0000 || iy >= 0x5ffb0000) {
183		/* region 5: x+1 and y are both (>= sqrt(max)/8) i.e. 2**8188 */
184		t = x / y;
185		LD_RE(ans) = atanl(t);
186		LD_IM(ans) = ln2 + logl(y) + half * log1pl(t * t);
187	} else if (x < Foursqrtu) {
188		/* region 6: x is very small, < 4sqrt(min) */
189		A = sqrtl(one + y * y);
190		LD_RE(ans) = x / A;	/* may underflow */
191		if (iy >= 0x3fff8000)	/* if y > Acrossover */
192			LD_IM(ans) = logl(y + A);
193		else
194			LD_IM(ans) = half * log1pl((y + y) * (y + A));
195	} else {	/* safe region */
196		y2 = y * y;
197		xp1 = x + one;
198		xm1 = x - one;
199		R = sqrtl(xp1 * xp1 + y2);
200		S = sqrtl(xm1 * xm1 + y2);
201		A = half * (R + S);
202		B = x / A;
203		if (B <= Bcrossover)
204			LD_RE(ans) = asinl(B);
205		else {		/* use atan and an accurate approx to a-x */
206			Apx = A + x;
207			if (x <= one)
208				LD_RE(ans) = atanl(x / sqrtl(half * Apx * (y2 /
209					(R + xp1) + (S - xm1))));
210			else
211				LD_RE(ans) = atanl(x / (y * sqrtl(half * (Apx /
212					(R + xp1) + Apx / (S + xm1)))));
213		}
214		if (A <= Acrossover) {
215			/* use log1p and an accurate approx to A-1 */
216			if (x < one)
217				Am1 = half * (y2 / (R + xp1) + y2 / (S - xm1));
218			else
219				Am1 = half * (y2 / (R + xp1) + (S + xm1));
220			LD_IM(ans) = log1pl(Am1 + sqrtl(Am1 * (A + one)));
221		} else {
222			LD_IM(ans) = logl(A + sqrtl(A * A - one));
223		}
224	}
225
226	if (hx < 0)
227		LD_RE(ans) = -LD_RE(ans);
228	if (hy < 0)
229		LD_IM(ans) = -LD_IM(ans);
230
231	return (ans);
232}
233