125c28e8Piotr Jasiukajtis/*
225c28e8Piotr Jasiukajtis * CDDL HEADER START
325c28e8Piotr Jasiukajtis *
425c28e8Piotr Jasiukajtis * The contents of this file are subject to the terms of the
525c28e8Piotr Jasiukajtis * Common Development and Distribution License (the "License").
625c28e8Piotr Jasiukajtis * You may not use this file except in compliance with the License.
725c28e8Piotr Jasiukajtis *
825c28e8Piotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
925c28e8Piotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
1025c28e8Piotr Jasiukajtis * See the License for the specific language governing permissions
1125c28e8Piotr Jasiukajtis * and limitations under the License.
1225c28e8Piotr Jasiukajtis *
1325c28e8Piotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
1425c28e8Piotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
1525c28e8Piotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
1625c28e8Piotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
1725c28e8Piotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
1825c28e8Piotr Jasiukajtis *
1925c28e8Piotr Jasiukajtis * CDDL HEADER END
2025c28e8Piotr Jasiukajtis */
2125c28e8Piotr Jasiukajtis
2225c28e8Piotr Jasiukajtis/*
2325c28e8Piotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
2425c28e8Piotr Jasiukajtis */
2525c28e8Piotr Jasiukajtis/*
2625c28e8Piotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
2725c28e8Piotr Jasiukajtis * Use is subject to license terms.
2825c28e8Piotr Jasiukajtis */
2925c28e8Piotr Jasiukajtis
30ddc0e0bRichard Lowe#pragma weak __catanl = catanl
3125c28e8Piotr Jasiukajtis
3225c28e8Piotr Jasiukajtis/* INDENT OFF */
3325c28e8Piotr Jasiukajtis/*
3425c28e8Piotr Jasiukajtis * ldcomplex catanl(ldcomplex z);
3525c28e8Piotr Jasiukajtis *
3625c28e8Piotr Jasiukajtis * Atan(z) return A + Bi where,
3725c28e8Piotr Jasiukajtis *            1
3825c28e8Piotr Jasiukajtis *	A = --- * atan2(2x, 1-x*x-y*y)
3925c28e8Piotr Jasiukajtis *            2
4025c28e8Piotr Jasiukajtis *
4125c28e8Piotr Jasiukajtis *            1      [ x*x + (y+1)*(y+1) ]   1               4y
4225c28e8Piotr Jasiukajtis *       B = --- log [ ----------------- ] = - log (1+ -----------------)
4325c28e8Piotr Jasiukajtis *            4      [ x*x + (y-1)*(y-1) ]   4         x*x + (y-1)*(y-1)
4425c28e8Piotr Jasiukajtis *
4525c28e8Piotr Jasiukajtis *                 2    16  3                         y
4625c28e8Piotr Jasiukajtis *         = t - 2t   + -- t  - ..., where t = -----------------
4725c28e8Piotr Jasiukajtis *                      3                      x*x + (y-1)*(y-1)
4825c28e8Piotr Jasiukajtis * Proof:
4925c28e8Piotr Jasiukajtis * Let w = atan(z=x+yi) = A + B i. Then tan(w) = z.
5025c28e8Piotr Jasiukajtis * Since sin(w) = (exp(iw)-exp(-iw))/(2i), cos(w)=(exp(iw)+exp(-iw))/(2),
5125c28e8Piotr Jasiukajtis * Let p = exp(iw), then z = tan(w) = ((p-1/p)/(p+1/p))/i, or
5225c28e8Piotr Jasiukajtis * iz = (p*p-1)/(p*p+1), or, after simplification,
5325c28e8Piotr Jasiukajtis *	p*p = (1+iz)/(1-iz)			            ... (1)
5425c28e8Piotr Jasiukajtis * LHS of (1) = exp(2iw) = exp(2i(A+Bi)) = exp(-2B)*exp(2iA)
5525c28e8Piotr Jasiukajtis *            = exp(-2B)*(cos(2A)+i*sin(2A))	            ... (2)
5625c28e8Piotr Jasiukajtis *              1-y+ix   (1-y+ix)*(1+y+ix)   1-x*x-y*y + 2xi
5725c28e8Piotr Jasiukajtis * RHS of (1) = ------ = ----------------- = --------------- ... (3)
5825c28e8Piotr Jasiukajtis *              1+y-ix    (1+y)**2 + x**2    (1+y)**2 + x**2
5925c28e8Piotr Jasiukajtis *
6025c28e8Piotr Jasiukajtis * Comparing the real and imaginary parts of (2) and (3), we have:
6125c28e8Piotr Jasiukajtis * 	cos(2A) : 1-x*x-y*y = sin(2A) : 2x
6225c28e8Piotr Jasiukajtis * and hence
6325c28e8Piotr Jasiukajtis *	tan(2A) = 2x/(1-x*x-y*y), or
6425c28e8Piotr Jasiukajtis *	A = 0.5 * atan2(2x, 1-x*x-y*y)	                    ... (4)
6525c28e8Piotr Jasiukajtis *
6625c28e8Piotr Jasiukajtis * For the imaginary part B, Note that |p*p| = exp(-2B), and
6725c28e8Piotr Jasiukajtis *	|1+iz|   |i-z|   hypot(x,(y-1))
6825c28e8Piotr Jasiukajtis *       |----| = |---| = --------------
6925c28e8Piotr Jasiukajtis *	|1-iz|   |i+z|   hypot(x,(y+1))
7025c28e8Piotr Jasiukajtis * Thus
7125c28e8Piotr Jasiukajtis *                 x*x + (y+1)*(y+1)
7225c28e8Piotr Jasiukajtis *	exp(4B) = -----------------, or
7325c28e8Piotr Jasiukajtis *                 x*x + (y-1)*(y-1)
7425c28e8Piotr Jasiukajtis *
7525c28e8Piotr Jasiukajtis *            1     [x^2+(y+1)^2]   1             4y
7625c28e8Piotr Jasiukajtis *       B =  - log [-----------] = - log(1+ -------------)  ... (5)
7725c28e8Piotr Jasiukajtis *            4     [x^2+(y-1)^2]   4         x^2+(y-1)^2
7825c28e8Piotr Jasiukajtis *
7925c28e8Piotr Jasiukajtis * QED.
8025c28e8Piotr Jasiukajtis *
8125c28e8Piotr Jasiukajtis * Note that: if catan( x, y) = ( u, v), then
8225c28e8Piotr Jasiukajtis *               catan(-x, y) = (-u, v)
8325c28e8Piotr Jasiukajtis *               catan( x,-y) = ( u,-v)
8425c28e8Piotr Jasiukajtis *
8525c28e8Piotr Jasiukajtis * Also,   catan(x,y) = -i*catanh(-y,x), or
8625c28e8Piotr Jasiukajtis *        catanh(x,y) =  i*catan(-y,x)
8725c28e8Piotr Jasiukajtis * So, if catanh(y,x) = (v,u), then catan(x,y) = -i*(-v,u) = (u,v), i.e.,
8825c28e8Piotr Jasiukajtis *         catan(x,y) = (u,v)
8925c28e8Piotr Jasiukajtis *
9025c28e8Piotr Jasiukajtis * EXCEPTION CASES (conform to ISO/IEC 9899:1999(E)):
9125c28e8Piotr Jasiukajtis *    catan( 0  , 0   ) =  (0    ,  0   )
9225c28e8Piotr Jasiukajtis *    catan( NaN, 0   ) =  (NaN  ,  0   )
9325c28e8Piotr Jasiukajtis *    catan( 0  , 1   ) =  (0    ,  +inf) with divide-by-zero
9425c28e8Piotr Jasiukajtis *    catan( inf, y   ) =  (pi/2 ,  0   ) for finite +y
9525c28e8Piotr Jasiukajtis *    catan( NaN, y   ) =  (NaN  ,  NaN ) with invalid for finite y != 0
9625c28e8Piotr Jasiukajtis *    catan( x  , inf ) =  (pi/2 ,  0   ) for finite +x
9725c28e8Piotr Jasiukajtis *    catan( inf, inf ) =  (pi/2 ,  0   )
9825c28e8Piotr Jasiukajtis *    catan( NaN, inf ) =  (NaN  ,  0   )
9925c28e8Piotr Jasiukajtis *    catan( x  , NaN ) =  (NaN  ,  NaN ) with invalid for finite x
10025c28e8Piotr Jasiukajtis *    catan( inf, NaN ) =  (pi/2 ,  +-0 )
10125c28e8Piotr Jasiukajtis */
10225c28e8Piotr Jasiukajtis/* INDENT ON */
10325c28e8Piotr Jasiukajtis
10425c28e8Piotr Jasiukajtis#include "libm.h"	/* atan2l/atanl/fabsl/isinfl/iszerol/log1pl/logl */
10525c28e8Piotr Jasiukajtis#include "complex_wrapper.h"
10625c28e8Piotr Jasiukajtis#include "longdouble.h"
10725c28e8Piotr Jasiukajtis
10825c28e8Piotr Jasiukajtis/* INDENT OFF */
10925c28e8Piotr Jasiukajtisstatic const long double
11025c28e8Piotr Jasiukajtiszero = 0.0L,
11125c28e8Piotr Jasiukajtisone = 1.0L,
11225c28e8Piotr Jasiukajtistwo = 2.0L,
11325c28e8Piotr Jasiukajtishalf = 0.5L,
11425c28e8Piotr Jasiukajtisln2 = 6.931471805599453094172321214581765680755e-0001L,
11525c28e8Piotr Jasiukajtispi_2 = 1.570796326794896619231321691639751442098584699687552910487472L,
11625c28e8Piotr Jasiukajtis#if defined(__x86)
11725c28e8Piotr JasiukajtisE = 2.910383045673370361328125000000000000000e-11L,	/* 2**-35 */
11825c28e8Piotr JasiukajtisEinv = 3.435973836800000000000000000000000000000e+10L;	/* 2**+35 */
11925c28e8Piotr Jasiukajtis#else
12025c28e8Piotr JasiukajtisE = 8.673617379884035472059622406959533691406e-19L,	/* 2**-60 */
12125c28e8Piotr JasiukajtisEinv = 1.152921504606846976000000000000000000000e18L;	/* 2**+60 */
12225c28e8Piotr Jasiukajtis#endif
12325c28e8Piotr Jasiukajtis/* INDENT ON */
12425c28e8Piotr Jasiukajtis
12525c28e8Piotr Jasiukajtisldcomplex
12625c28e8Piotr Jasiukajtiscatanl(ldcomplex z) {
12725c28e8Piotr Jasiukajtis	ldcomplex ans;
12825c28e8Piotr Jasiukajtis	long double x, y, t1, ax, ay, t;
12925c28e8Piotr Jasiukajtis	int hx, hy, ix, iy;
13025c28e8Piotr Jasiukajtis
13125c28e8Piotr Jasiukajtis	x = LD_RE(z);
13225c28e8Piotr Jasiukajtis	y = LD_IM(z);
13325c28e8Piotr Jasiukajtis	ax = fabsl(x);
13425c28e8Piotr Jasiukajtis	ay = fabsl(y);
13525c28e8Piotr Jasiukajtis	hx = HI_XWORD(x);
13625c28e8Piotr Jasiukajtis	hy = HI_XWORD(y);
13725c28e8Piotr Jasiukajtis	ix = hx & 0x7fffffff;
13825c28e8Piotr Jasiukajtis	iy = hy & 0x7fffffff;
13925c28e8Piotr Jasiukajtis
14025c28e8Piotr Jasiukajtis	/* x is inf or NaN */
14125c28e8Piotr Jasiukajtis	if (ix >= 0x7fff0000) {
14225c28e8Piotr Jasiukajtis		if (isinfl(x)) {
14325c28e8Piotr Jasiukajtis			LD_RE(ans) = pi_2;
14425c28e8Piotr Jasiukajtis			LD_IM(ans) = zero;
14525c28e8Piotr Jasiukajtis		} else {
14625c28e8Piotr Jasiukajtis			LD_RE(ans) = x + x;
14725c28e8Piotr Jasiukajtis			if (iszerol(y) || (isinfl(y)))
14825c28e8Piotr Jasiukajtis				LD_IM(ans) = zero;
14925c28e8Piotr Jasiukajtis			else
15025c28e8Piotr Jasiukajtis				LD_IM(ans) = (fabsl(y) - ay) / (fabsl(y) - ay);
15125c28e8Piotr Jasiukajtis		}
15225c28e8Piotr Jasiukajtis	} else if (iy >= 0x7fff0000) {
15325c28e8Piotr Jasiukajtis		/* y is inf or NaN */
15425c28e8Piotr Jasiukajtis		if (isinfl(y)) {
15525c28e8Piotr Jasiukajtis			LD_RE(ans) = pi_2;
15625c28e8Piotr Jasiukajtis			LD_IM(ans) = zero;
15725c28e8Piotr Jasiukajtis		} else {
15825c28e8Piotr Jasiukajtis			LD_RE(ans) = (fabsl(x) - ax) / (fabsl(x) - ax);
15925c28e8Piotr Jasiukajtis			LD_IM(ans) = y;
16025c28e8Piotr Jasiukajtis		}
16125c28e8Piotr Jasiukajtis	} else if (iszerol(x)) {
16225c28e8Piotr Jasiukajtis		/* INDENT OFF */
16325c28e8Piotr Jasiukajtis		/*
16425c28e8Piotr Jasiukajtis		 * x = 0
16525c28e8Piotr Jasiukajtis		 *      1                            1
16625c28e8Piotr Jasiukajtis		 * A = --- * atan2(2x, 1-x*x-y*y) = --- atan2(0,1-|y|)
16725c28e8Piotr Jasiukajtis		 *      2                            2
16825c28e8Piotr Jasiukajtis		 *
16925c28e8Piotr Jasiukajtis		 *     1     [ (y+1)*(y+1) ]   1          2      1         2y
17025c28e8Piotr Jasiukajtis		 * B = - log [ ----------- ] = - log (1+ ---) or - log(1+ ----)
17125c28e8Piotr Jasiukajtis		 *     4     [ (y-1)*(y-1) ]   2         y-1     2         1-y
17225c28e8Piotr Jasiukajtis		 */
17325c28e8Piotr Jasiukajtis		/* INDENT ON */
17425c28e8Piotr Jasiukajtis		t = one - ay;
17525c28e8Piotr Jasiukajtis		if (ay == one) {
17625c28e8Piotr Jasiukajtis			/* y=1: catan(0,1)=(0,+inf) with 1/0 signal */
17725c28e8Piotr Jasiukajtis			LD_IM(ans) = ay / ax;
17825c28e8Piotr Jasiukajtis			LD_RE(ans) = zero;
17925c28e8Piotr Jasiukajtis		} else if (ay > one) {	/* y>1 */
18025c28e8Piotr Jasiukajtis			LD_IM(ans) = half * log1pl(two / (-t));
18125c28e8Piotr Jasiukajtis			LD_RE(ans) = pi_2;
18225c28e8Piotr Jasiukajtis		} else {		/* y<1 */
18325c28e8Piotr Jasiukajtis			LD_IM(ans) = half * log1pl((ay + ay) / t);
18425c28e8Piotr Jasiukajtis			LD_RE(ans) = zero;
18525c28e8Piotr Jasiukajtis		}
18625c28e8Piotr Jasiukajtis	} else if (ay < E * (one + ax)) {
18725c28e8Piotr Jasiukajtis		/* INDENT OFF */
18825c28e8Piotr Jasiukajtis		/*
18925c28e8Piotr Jasiukajtis		 * Tiny y (relative to 1+|x|)
19025c28e8Piotr Jasiukajtis		 *     |y| < E*(1+|x|)
19125c28e8Piotr Jasiukajtis		 * where E=2**-29, -35, -60 for double, extended, quad precision
19225c28e8Piotr Jasiukajtis		 *
19325c28e8Piotr Jasiukajtis		 *     1                         [x<=1:   atan(x)
19425c28e8Piotr Jasiukajtis		 * A = - * atan2(2x,1-x*x-y*y) ~ [      1                 1+x
19525c28e8Piotr Jasiukajtis		 *     2                         [x>=1: - atan2(2,(1-x)*(-----))
19625c28e8Piotr Jasiukajtis		 *                                      2                  x
19725c28e8Piotr Jasiukajtis		 *
19825c28e8Piotr Jasiukajtis		 *                               y/x
19925c28e8Piotr Jasiukajtis		 * B ~ t*(1-2t), where t = ----------------- is tiny
20025c28e8Piotr Jasiukajtis		 *                         x + (y-1)*(y-1)/x
20125c28e8Piotr Jasiukajtis		 *
20225c28e8Piotr Jasiukajtis		 *                           y
20325c28e8Piotr Jasiukajtis		 * (when x < 2**-60, t = ----------- )
20425c28e8Piotr Jasiukajtis		 *                       (y-1)*(y-1)
20525c28e8Piotr Jasiukajtis		 */
20625c28e8Piotr Jasiukajtis		/* INDENT ON */
20725c28e8Piotr Jasiukajtis		if (ay == zero)
20825c28e8Piotr Jasiukajtis			LD_IM(ans) = ay;
20925c28e8Piotr Jasiukajtis		else {
21025c28e8Piotr Jasiukajtis			t1 = ay - one;
21125c28e8Piotr Jasiukajtis			if (ix < 0x3fc30000)
21225c28e8Piotr Jasiukajtis				t = ay / (t1 * t1);
21325c28e8Piotr Jasiukajtis			else if (ix > 0x403b0000)
21425c28e8Piotr Jasiukajtis				t = (ay / ax) / ax;
21525c28e8Piotr Jasiukajtis			else
21625c28e8Piotr Jasiukajtis				t = ay / (ax * ax + t1 * t1);
21725c28e8Piotr Jasiukajtis			LD_IM(ans) = t * (one - two * t);
21825c28e8Piotr Jasiukajtis		}
21925c28e8Piotr Jasiukajtis		if (ix < 0x3fff0000)
22025c28e8Piotr Jasiukajtis			LD_RE(ans) = atanl(ax);
22125c28e8Piotr Jasiukajtis		else
22225c28e8Piotr Jasiukajtis			LD_RE(ans) = half * atan2l(two, (one - ax) * (one +
22325c28e8Piotr Jasiukajtis				one / ax));
22425c28e8Piotr Jasiukajtis
22525c28e8Piotr Jasiukajtis	} else if (ay > Einv * (one + ax)) {
22625c28e8Piotr Jasiukajtis		/* INDENT OFF */
22725c28e8Piotr Jasiukajtis		/*
22825c28e8Piotr Jasiukajtis		 * Huge y relative to 1+|x|
22925c28e8Piotr Jasiukajtis		 *     |y| > Einv*(1+|x|), where Einv~2**(prec/2+3),
23025c28e8Piotr Jasiukajtis		 *      1
23125c28e8Piotr Jasiukajtis		 * A ~ --- * atan2(2x, -y*y) ~ pi/2
23225c28e8Piotr Jasiukajtis		 *      2
23325c28e8Piotr Jasiukajtis		 *                               y
23425c28e8Piotr Jasiukajtis		 * B ~ t*(1-2t), where t = --------------- is tiny
23525c28e8Piotr Jasiukajtis		 *                          (y-1)*(y-1)
23625c28e8Piotr Jasiukajtis		 */
23725c28e8Piotr Jasiukajtis		/* INDENT ON */
23825c28e8Piotr Jasiukajtis		LD_RE(ans) = pi_2;
23925c28e8Piotr Jasiukajtis		t = (ay / (ay - one)) / (ay - one);
24025c28e8Piotr Jasiukajtis		LD_IM(ans) = t * (one - (t + t));
24125c28e8Piotr Jasiukajtis	} else if (ay == one) {
24225c28e8Piotr Jasiukajtis		/* INDENT OFF */
24325c28e8Piotr Jasiukajtis		/*
24425c28e8Piotr Jasiukajtis		 * y=1
24525c28e8Piotr Jasiukajtis		 *     1                      1
24625c28e8Piotr Jasiukajtis		 * A = - * atan2(2x, -x*x) = --- atan2(2,-x)
24725c28e8Piotr Jasiukajtis		 *     2                      2
24825c28e8Piotr Jasiukajtis		 *
24925c28e8Piotr Jasiukajtis		 *     1     [ x*x+4]   1          4     [ 0.5(log2-logx) if
25025c28e8Piotr Jasiukajtis		 * B = - log [ -----] = - log (1+ ---) = [ |x|<E, else 0.25*
25125c28e8Piotr Jasiukajtis		 *     4     [  x*x ]   4         x*x    [ log1p((2/x)*(2/x))
25225c28e8Piotr Jasiukajtis		 */
25325c28e8Piotr Jasiukajtis		/* INDENT ON */
25425c28e8Piotr Jasiukajtis		LD_RE(ans) = half * atan2l(two, -ax);
25525c28e8Piotr Jasiukajtis		if (ax < E)
25625c28e8Piotr Jasiukajtis			LD_IM(ans) = half * (ln2 - logl(ax));
25725c28e8Piotr Jasiukajtis		else {
25825c28e8Piotr Jasiukajtis			t = two / ax;
25925c28e8Piotr Jasiukajtis			LD_IM(ans) = 0.25L * log1pl(t * t);
26025c28e8Piotr Jasiukajtis		}
26125c28e8Piotr Jasiukajtis	} else if (ax > Einv * Einv) {
26225c28e8Piotr Jasiukajtis		/* INDENT OFF */
26325c28e8Piotr Jasiukajtis		/*
26425c28e8Piotr Jasiukajtis		 * Huge x:
26525c28e8Piotr Jasiukajtis		 * when |x| > 1/E^2,
26625c28e8Piotr Jasiukajtis		 *      1                           pi
26725c28e8Piotr Jasiukajtis		 * A ~ --- * atan2(2x, -x*x-y*y) ~ ---
26825c28e8Piotr Jasiukajtis		 *      2                           2
26925c28e8Piotr Jasiukajtis		 *                               y                 y/x
27025c28e8Piotr Jasiukajtis		 * B ~ t*(1-2t), where t = --------------- = (-------------- )/x
27125c28e8Piotr Jasiukajtis		 *                         x*x+(y-1)*(y-1)     1+((y-1)/x)^2
27225c28e8Piotr Jasiukajtis		 */
27325c28e8Piotr Jasiukajtis		/* INDENT ON */
27425c28e8Piotr Jasiukajtis		LD_RE(ans) = pi_2;
27525c28e8Piotr Jasiukajtis		t = ((ay / ax) / (one + ((ay - one) / ax) * ((ay - one) /
27625c28e8Piotr Jasiukajtis			ax))) / ax;
27725c28e8Piotr Jasiukajtis		LD_IM(ans) = t * (one - (t + t));
27825c28e8Piotr Jasiukajtis	} else if (ax < E * E * E * E) {
27925c28e8Piotr Jasiukajtis		/* INDENT OFF */
28025c28e8Piotr Jasiukajtis		/*
28125c28e8Piotr Jasiukajtis		 * Tiny x:
28225c28e8Piotr Jasiukajtis		 * when |x| < E^4,  (note that y != 1)
28325c28e8Piotr Jasiukajtis		 *      1                            1
28425c28e8Piotr Jasiukajtis		 * A = --- * atan2(2x, 1-x*x-y*y) ~ --- * atan2(2x,1-y*y)
28525c28e8Piotr Jasiukajtis		 *      2                            2
28625c28e8Piotr Jasiukajtis		 *
28725c28e8Piotr Jasiukajtis		 *     1     [ (y+1)*(y+1) ]   1          2      1         2y
28825c28e8Piotr Jasiukajtis		 * B = - log [ ----------- ] = - log (1+ ---) or - log(1+ ----)
28925c28e8Piotr Jasiukajtis		 *     4     [ (y-1)*(y-1) ]   2         y-1     2         1-y
29025c28e8Piotr Jasiukajtis		 */
29125c28e8Piotr Jasiukajtis		/* INDENT ON */
29225c28e8Piotr Jasiukajtis		LD_RE(ans) = half * atan2l(ax + ax, (one - ay) * (one + ay));
29325c28e8Piotr Jasiukajtis		if (ay > one)	/* y>1 */
29425c28e8Piotr Jasiukajtis			LD_IM(ans) = half * log1pl(two / (ay - one));
29525c28e8Piotr Jasiukajtis		else		/* y<1 */
29625c28e8Piotr Jasiukajtis			LD_IM(ans) = half * log1pl((ay + ay) / (one - ay));
29725c28e8Piotr Jasiukajtis	} else {
29825c28e8Piotr Jasiukajtis		/* INDENT OFF */
29925c28e8Piotr Jasiukajtis		/*
30025c28e8Piotr Jasiukajtis		 * normal x,y
30125c28e8Piotr Jasiukajtis		 *      1
30225c28e8Piotr Jasiukajtis		 * A = --- * atan2(2x, 1-x*x-y*y)
30325c28e8Piotr Jasiukajtis		 *      2
30425c28e8Piotr Jasiukajtis		 *
30525c28e8Piotr Jasiukajtis		 *     1     [ x*x+(y+1)*(y+1) ]   1               4y
30625c28e8Piotr Jasiukajtis		 * B = - log [ --------------- ] = - log (1+ -----------------)
30725c28e8Piotr Jasiukajtis		 *     4     [ x*x+(y-1)*(y-1) ]   4         x*x + (y-1)*(y-1)
30825c28e8Piotr Jasiukajtis		 */
30925c28e8Piotr Jasiukajtis		/* INDENT ON */
31025c28e8Piotr Jasiukajtis		t = one - ay;
31125c28e8Piotr Jasiukajtis		if (iy >= 0x3ffe0000 && iy < 0x40000000) {
31225c28e8Piotr Jasiukajtis			/* y close to 1 */
31325c28e8Piotr Jasiukajtis			LD_RE(ans) = half * (atan2l((ax + ax), (t * (one +
31425c28e8Piotr Jasiukajtis				ay) - ax * ax)));
31525c28e8Piotr Jasiukajtis		} else if (ix >= 0x3ffe0000 && ix < 0x40000000) {
31625c28e8Piotr Jasiukajtis			/* x close to 1 */
31725c28e8Piotr Jasiukajtis			LD_RE(ans) = half * atan2l((ax + ax), ((one - ax) *
31825c28e8Piotr Jasiukajtis				(one + ax) - ay * ay));
31925c28e8Piotr Jasiukajtis		} else
32025c28e8Piotr Jasiukajtis			LD_RE(ans) = half * atan2l((ax + ax), ((one - ax *
32125c28e8Piotr Jasiukajtis				ax) - ay * ay));
32225c28e8Piotr Jasiukajtis		LD_IM(ans) = 0.25L * log1pl((4.0L * ay) / (ax * ax + t * t));
32325c28e8Piotr Jasiukajtis	}
32425c28e8Piotr Jasiukajtis	if (hx < 0)
32525c28e8Piotr Jasiukajtis		LD_RE(ans) = -LD_RE(ans);
32625c28e8Piotr Jasiukajtis	if (hy < 0)
32725c28e8Piotr Jasiukajtis		LD_IM(ans) = -LD_IM(ans);
32825c28e8Piotr Jasiukajtis	return (ans);
32925c28e8Piotr Jasiukajtis}
330