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 __catanl = catanl
31
32/* INDENT OFF */
33/*
34 * ldcomplex catanl(ldcomplex z);
35 *
36 * Atan(z) return A + Bi where,
37 *            1
38 *	A = --- * atan2(2x, 1-x*x-y*y)
39 *            2
40 *
41 *            1      [ x*x + (y+1)*(y+1) ]   1               4y
42 *       B = --- log [ ----------------- ] = - log (1+ -----------------)
43 *            4      [ x*x + (y-1)*(y-1) ]   4         x*x + (y-1)*(y-1)
44 *
45 *                 2    16  3                         y
46 *         = t - 2t   + -- t  - ..., where t = -----------------
47 *                      3                      x*x + (y-1)*(y-1)
48 * Proof:
49 * Let w = atan(z=x+yi) = A + B i. Then tan(w) = z.
50 * Since sin(w) = (exp(iw)-exp(-iw))/(2i), cos(w)=(exp(iw)+exp(-iw))/(2),
51 * Let p = exp(iw), then z = tan(w) = ((p-1/p)/(p+1/p))/i, or
52 * iz = (p*p-1)/(p*p+1), or, after simplification,
53 *	p*p = (1+iz)/(1-iz)			            ... (1)
54 * LHS of (1) = exp(2iw) = exp(2i(A+Bi)) = exp(-2B)*exp(2iA)
55 *            = exp(-2B)*(cos(2A)+i*sin(2A))	            ... (2)
56 *              1-y+ix   (1-y+ix)*(1+y+ix)   1-x*x-y*y + 2xi
57 * RHS of (1) = ------ = ----------------- = --------------- ... (3)
58 *              1+y-ix    (1+y)**2 + x**2    (1+y)**2 + x**2
59 *
60 * Comparing the real and imaginary parts of (2) and (3), we have:
61 * 	cos(2A) : 1-x*x-y*y = sin(2A) : 2x
62 * and hence
63 *	tan(2A) = 2x/(1-x*x-y*y), or
64 *	A = 0.5 * atan2(2x, 1-x*x-y*y)	                    ... (4)
65 *
66 * For the imaginary part B, Note that |p*p| = exp(-2B), and
67 *	|1+iz|   |i-z|   hypot(x,(y-1))
68 *       |----| = |---| = --------------
69 *	|1-iz|   |i+z|   hypot(x,(y+1))
70 * Thus
71 *                 x*x + (y+1)*(y+1)
72 *	exp(4B) = -----------------, or
73 *                 x*x + (y-1)*(y-1)
74 *
75 *            1     [x^2+(y+1)^2]   1             4y
76 *       B =  - log [-----------] = - log(1+ -------------)  ... (5)
77 *            4     [x^2+(y-1)^2]   4         x^2+(y-1)^2
78 *
79 * QED.
80 *
81 * Note that: if catan( x, y) = ( u, v), then
82 *               catan(-x, y) = (-u, v)
83 *               catan( x,-y) = ( u,-v)
84 *
85 * Also,   catan(x,y) = -i*catanh(-y,x), or
86 *        catanh(x,y) =  i*catan(-y,x)
87 * So, if catanh(y,x) = (v,u), then catan(x,y) = -i*(-v,u) = (u,v), i.e.,
88 *         catan(x,y) = (u,v)
89 *
90 * EXCEPTION CASES (conform to ISO/IEC 9899:1999(E)):
91 *    catan( 0  , 0   ) =  (0    ,  0   )
92 *    catan( NaN, 0   ) =  (NaN  ,  0   )
93 *    catan( 0  , 1   ) =  (0    ,  +inf) with divide-by-zero
94 *    catan( inf, y   ) =  (pi/2 ,  0   ) for finite +y
95 *    catan( NaN, y   ) =  (NaN  ,  NaN ) with invalid for finite y != 0
96 *    catan( x  , inf ) =  (pi/2 ,  0   ) for finite +x
97 *    catan( inf, inf ) =  (pi/2 ,  0   )
98 *    catan( NaN, inf ) =  (NaN  ,  0   )
99 *    catan( x  , NaN ) =  (NaN  ,  NaN ) with invalid for finite x
100 *    catan( inf, NaN ) =  (pi/2 ,  +-0 )
101 */
102/* INDENT ON */
103
104#include "libm.h"	/* atan2l/atanl/fabsl/isinfl/iszerol/log1pl/logl */
105#include "complex_wrapper.h"
106#include "longdouble.h"
107
108/* INDENT OFF */
109static const long double
110zero = 0.0L,
111one = 1.0L,
112two = 2.0L,
113half = 0.5L,
114ln2 = 6.931471805599453094172321214581765680755e-0001L,
115pi_2 = 1.570796326794896619231321691639751442098584699687552910487472L,
116#if defined(__x86)
117E = 2.910383045673370361328125000000000000000e-11L,	/* 2**-35 */
118Einv = 3.435973836800000000000000000000000000000e+10L;	/* 2**+35 */
119#else
120E = 8.673617379884035472059622406959533691406e-19L,	/* 2**-60 */
121Einv = 1.152921504606846976000000000000000000000e18L;	/* 2**+60 */
122#endif
123/* INDENT ON */
124
125ldcomplex
126catanl(ldcomplex z) {
127	ldcomplex ans;
128	long double x, y, t1, ax, ay, t;
129	int hx, hy, ix, iy;
130
131	x = LD_RE(z);
132	y = LD_IM(z);
133	ax = fabsl(x);
134	ay = fabsl(y);
135	hx = HI_XWORD(x);
136	hy = HI_XWORD(y);
137	ix = hx & 0x7fffffff;
138	iy = hy & 0x7fffffff;
139
140	/* x is inf or NaN */
141	if (ix >= 0x7fff0000) {
142		if (isinfl(x)) {
143			LD_RE(ans) = pi_2;
144			LD_IM(ans) = zero;
145		} else {
146			LD_RE(ans) = x + x;
147			if (iszerol(y) || (isinfl(y)))
148				LD_IM(ans) = zero;
149			else
150				LD_IM(ans) = (fabsl(y) - ay) / (fabsl(y) - ay);
151		}
152	} else if (iy >= 0x7fff0000) {
153		/* y is inf or NaN */
154		if (isinfl(y)) {
155			LD_RE(ans) = pi_2;
156			LD_IM(ans) = zero;
157		} else {
158			LD_RE(ans) = (fabsl(x) - ax) / (fabsl(x) - ax);
159			LD_IM(ans) = y;
160		}
161	} else if (iszerol(x)) {
162		/* INDENT OFF */
163		/*
164		 * x = 0
165		 *      1                            1
166		 * A = --- * atan2(2x, 1-x*x-y*y) = --- atan2(0,1-|y|)
167		 *      2                            2
168		 *
169		 *     1     [ (y+1)*(y+1) ]   1          2      1         2y
170		 * B = - log [ ----------- ] = - log (1+ ---) or - log(1+ ----)
171		 *     4     [ (y-1)*(y-1) ]   2         y-1     2         1-y
172		 */
173		/* INDENT ON */
174		t = one - ay;
175		if (ay == one) {
176			/* y=1: catan(0,1)=(0,+inf) with 1/0 signal */
177			LD_IM(ans) = ay / ax;
178			LD_RE(ans) = zero;
179		} else if (ay > one) {	/* y>1 */
180			LD_IM(ans) = half * log1pl(two / (-t));
181			LD_RE(ans) = pi_2;
182		} else {		/* y<1 */
183			LD_IM(ans) = half * log1pl((ay + ay) / t);
184			LD_RE(ans) = zero;
185		}
186	} else if (ay < E * (one + ax)) {
187		/* INDENT OFF */
188		/*
189		 * Tiny y (relative to 1+|x|)
190		 *     |y| < E*(1+|x|)
191		 * where E=2**-29, -35, -60 for double, extended, quad precision
192		 *
193		 *     1                         [x<=1:   atan(x)
194		 * A = - * atan2(2x,1-x*x-y*y) ~ [      1                 1+x
195		 *     2                         [x>=1: - atan2(2,(1-x)*(-----))
196		 *                                      2                  x
197		 *
198		 *                               y/x
199		 * B ~ t*(1-2t), where t = ----------------- is tiny
200		 *                         x + (y-1)*(y-1)/x
201		 *
202		 *                           y
203		 * (when x < 2**-60, t = ----------- )
204		 *                       (y-1)*(y-1)
205		 */
206		/* INDENT ON */
207		if (ay == zero)
208			LD_IM(ans) = ay;
209		else {
210			t1 = ay - one;
211			if (ix < 0x3fc30000)
212				t = ay / (t1 * t1);
213			else if (ix > 0x403b0000)
214				t = (ay / ax) / ax;
215			else
216				t = ay / (ax * ax + t1 * t1);
217			LD_IM(ans) = t * (one - two * t);
218		}
219		if (ix < 0x3fff0000)
220			LD_RE(ans) = atanl(ax);
221		else
222			LD_RE(ans) = half * atan2l(two, (one - ax) * (one +
223				one / ax));
224
225	} else if (ay > Einv * (one + ax)) {
226		/* INDENT OFF */
227		/*
228		 * Huge y relative to 1+|x|
229		 *     |y| > Einv*(1+|x|), where Einv~2**(prec/2+3),
230		 *      1
231		 * A ~ --- * atan2(2x, -y*y) ~ pi/2
232		 *      2
233		 *                               y
234		 * B ~ t*(1-2t), where t = --------------- is tiny
235		 *                          (y-1)*(y-1)
236		 */
237		/* INDENT ON */
238		LD_RE(ans) = pi_2;
239		t = (ay / (ay - one)) / (ay - one);
240		LD_IM(ans) = t * (one - (t + t));
241	} else if (ay == one) {
242		/* INDENT OFF */
243		/*
244		 * y=1
245		 *     1                      1
246		 * A = - * atan2(2x, -x*x) = --- atan2(2,-x)
247		 *     2                      2
248		 *
249		 *     1     [ x*x+4]   1          4     [ 0.5(log2-logx) if
250		 * B = - log [ -----] = - log (1+ ---) = [ |x|<E, else 0.25*
251		 *     4     [  x*x ]   4         x*x    [ log1p((2/x)*(2/x))
252		 */
253		/* INDENT ON */
254		LD_RE(ans) = half * atan2l(two, -ax);
255		if (ax < E)
256			LD_IM(ans) = half * (ln2 - logl(ax));
257		else {
258			t = two / ax;
259			LD_IM(ans) = 0.25L * log1pl(t * t);
260		}
261	} else if (ax > Einv * Einv) {
262		/* INDENT OFF */
263		/*
264		 * Huge x:
265		 * when |x| > 1/E^2,
266		 *      1                           pi
267		 * A ~ --- * atan2(2x, -x*x-y*y) ~ ---
268		 *      2                           2
269		 *                               y                 y/x
270		 * B ~ t*(1-2t), where t = --------------- = (-------------- )/x
271		 *                         x*x+(y-1)*(y-1)     1+((y-1)/x)^2
272		 */
273		/* INDENT ON */
274		LD_RE(ans) = pi_2;
275		t = ((ay / ax) / (one + ((ay - one) / ax) * ((ay - one) /
276			ax))) / ax;
277		LD_IM(ans) = t * (one - (t + t));
278	} else if (ax < E * E * E * E) {
279		/* INDENT OFF */
280		/*
281		 * Tiny x:
282		 * when |x| < E^4,  (note that y != 1)
283		 *      1                            1
284		 * A = --- * atan2(2x, 1-x*x-y*y) ~ --- * atan2(2x,1-y*y)
285		 *      2                            2
286		 *
287		 *     1     [ (y+1)*(y+1) ]   1          2      1         2y
288		 * B = - log [ ----------- ] = - log (1+ ---) or - log(1+ ----)
289		 *     4     [ (y-1)*(y-1) ]   2         y-1     2         1-y
290		 */
291		/* INDENT ON */
292		LD_RE(ans) = half * atan2l(ax + ax, (one - ay) * (one + ay));
293		if (ay > one)	/* y>1 */
294			LD_IM(ans) = half * log1pl(two / (ay - one));
295		else		/* y<1 */
296			LD_IM(ans) = half * log1pl((ay + ay) / (one - ay));
297	} else {
298		/* INDENT OFF */
299		/*
300		 * normal x,y
301		 *      1
302		 * A = --- * atan2(2x, 1-x*x-y*y)
303		 *      2
304		 *
305		 *     1     [ x*x+(y+1)*(y+1) ]   1               4y
306		 * B = - log [ --------------- ] = - log (1+ -----------------)
307		 *     4     [ x*x+(y-1)*(y-1) ]   4         x*x + (y-1)*(y-1)
308		 */
309		/* INDENT ON */
310		t = one - ay;
311		if (iy >= 0x3ffe0000 && iy < 0x40000000) {
312			/* y close to 1 */
313			LD_RE(ans) = half * (atan2l((ax + ax), (t * (one +
314				ay) - ax * ax)));
315		} else if (ix >= 0x3ffe0000 && ix < 0x40000000) {
316			/* x close to 1 */
317			LD_RE(ans) = half * atan2l((ax + ax), ((one - ax) *
318				(one + ax) - ay * ay));
319		} else
320			LD_RE(ans) = half * atan2l((ax + ax), ((one - ax *
321				ax) - ay * ay));
322		LD_IM(ans) = 0.25L * log1pl((4.0L * ay) / (ax * ax + t * t));
323	}
324	if (hx < 0)
325		LD_RE(ans) = -LD_RE(ans);
326	if (hy < 0)
327		LD_IM(ans) = -LD_IM(ans);
328	return (ans);
329}
330