/*
* CDDL HEADER START
*
* The contents of this file are subject to the terms of the
* Common Development and Distribution License (the "License").
* You may not use this file except in compliance with the License.
*
* You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
* See the License for the specific language governing permissions
* and limitations under the License.
*
* When distributing Covered Code, include this CDDL HEADER in each
* file and include the License file at usr/src/OPENSOLARIS.LICENSE.
* If applicable, add the following below this CDDL HEADER, with the
* fields enclosed by brackets "[]" replaced with your own identifying
* information: Portions Copyright [yyyy] [name of copyright owner]
*
* CDDL HEADER END
*/
/*
* Copyright 2011 Nexenta Systems, Inc. All rights reserved.
*/
/*
* Copyright 2006 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
/* INDENT OFF */
/*
* ldcomplex catanl(ldcomplex z);
*
* Atan(z) return A + Bi where,
* 1
* A = --- * atan2(2x, 1-x*x-y*y)
* 2
*
* 1 [ x*x + (y+1)*(y+1) ] 1 4y
* B = --- log [ ----------------- ] = - log (1+ -----------------)
* 4 [ x*x + (y-1)*(y-1) ] 4 x*x + (y-1)*(y-1)
*
* 2 16 3 y
* = t - 2t + -- t - ..., where t = -----------------
* 3 x*x + (y-1)*(y-1)
* Proof:
* Let w = atan(z=x+yi) = A + B i. Then tan(w) = z.
* Since sin(w) = (exp(iw)-exp(-iw))/(2i), cos(w)=(exp(iw)+exp(-iw))/(2),
* Let p = exp(iw), then z = tan(w) = ((p-1/p)/(p+1/p))/i, or
* iz = (p*p-1)/(p*p+1), or, after simplification,
* p*p = (1+iz)/(1-iz) ... (1)
* LHS of (1) = exp(2iw) = exp(2i(A+Bi)) = exp(-2B)*exp(2iA)
* = exp(-2B)*(cos(2A)+i*sin(2A)) ... (2)
* 1-y+ix (1-y+ix)*(1+y+ix) 1-x*x-y*y + 2xi
* RHS of (1) = ------ = ----------------- = --------------- ... (3)
* 1+y-ix (1+y)**2 + x**2 (1+y)**2 + x**2
*
* Comparing the real and imaginary parts of (2) and (3), we have:
* cos(2A) : 1-x*x-y*y = sin(2A) : 2x
* and hence
* tan(2A) = 2x/(1-x*x-y*y), or
* A = 0.5 * atan2(2x, 1-x*x-y*y) ... (4)
*
* For the imaginary part B, Note that |p*p| = exp(-2B), and
* |1+iz| |i-z| hypot(x,(y-1))
* |----| = |---| = --------------
* |1-iz| |i+z| hypot(x,(y+1))
* Thus
* x*x + (y+1)*(y+1)
* exp(4B) = -----------------, or
* x*x + (y-1)*(y-1)
*
* 1 [x^2+(y+1)^2] 1 4y
* B = - log [-----------] = - log(1+ -------------) ... (5)
* 4 [x^2+(y-1)^2] 4 x^2+(y-1)^2
*
* QED.
*
* Note that: if catan( x, y) = ( u, v), then
* catan(-x, y) = (-u, v)
* catan( x,-y) = ( u,-v)
*
* Also, catan(x,y) = -i*catanh(-y,x), or
* catanh(x,y) = i*catan(-y,x)
* So, if catanh(y,x) = (v,u), then catan(x,y) = -i*(-v,u) = (u,v), i.e.,
* catan(x,y) = (u,v)
*
* catan( 0 , 0 ) = (0 , 0 )
* catan( NaN, 0 ) = (NaN , 0 )
* catan( 0 , 1 ) = (0 , +inf) with divide-by-zero
* catan( inf, y ) = (pi/2 , 0 ) for finite +y
* catan( NaN, y ) = (NaN , NaN ) with invalid for finite y != 0
* catan( x , inf ) = (pi/2 , 0 ) for finite +x
* catan( inf, inf ) = (pi/2 , 0 )
* catan( NaN, inf ) = (NaN , 0 )
* catan( x , NaN ) = (NaN , NaN ) with invalid for finite x
* catan( inf, NaN ) = (pi/2 , +-0 )
*/
/* INDENT ON */
#include "complex_wrapper.h"
#include "longdouble.h"
/* INDENT OFF */
static const long double
#if defined(__x86)
E = 2.910383045673370361328125000000000000000e-11L, /* 2**-35 */
#else
E = 8.673617379884035472059622406959533691406e-19L, /* 2**-60 */
#endif
/* INDENT ON */
x = LD_RE(z);
y = LD_IM(z);
/* x is inf or NaN */
if (ix >= 0x7fff0000) {
if (isinfl(x)) {
} else {
else
}
} else if (iy >= 0x7fff0000) {
/* y is inf or NaN */
if (isinfl(y)) {
} else {
}
} else if (iszerol(x)) {
/* INDENT OFF */
/*
* x = 0
* 1 1
* A = --- * atan2(2x, 1-x*x-y*y) = --- atan2(0,1-|y|)
* 2 2
*
* 1 [ (y+1)*(y+1) ] 1 2 1 2y
* B = - log [ ----------- ] = - log (1+ ---) or - log(1+ ----)
* 4 [ (y-1)*(y-1) ] 2 y-1 2 1-y
*/
/* INDENT ON */
/* y=1: catan(0,1)=(0,+inf) with 1/0 signal */
} else { /* y<1 */
}
/* INDENT OFF */
/*
* Tiny y (relative to 1+|x|)
* |y| < E*(1+|x|)
* where E=2**-29, -35, -60 for double, extended, quad precision
*
* 1 [x<=1: atan(x)
* A = - * atan2(2x,1-x*x-y*y) ~ [ 1 1+x
* 2 [x>=1: - atan2(2,(1-x)*(-----))
* 2 x
*
* y/x
* B ~ t*(1-2t), where t = ----------------- is tiny
* x + (y-1)*(y-1)/x
*
* y
* (when x < 2**-60, t = ----------- )
* (y-1)*(y-1)
*/
/* INDENT ON */
else {
if (ix < 0x3fc30000)
else if (ix > 0x403b0000)
else
}
if (ix < 0x3fff0000)
else
/* INDENT OFF */
/*
* Huge y relative to 1+|x|
* |y| > Einv*(1+|x|), where Einv~2**(prec/2+3),
* 1
* A ~ --- * atan2(2x, -y*y) ~ pi/2
* 2
* y
* B ~ t*(1-2t), where t = --------------- is tiny
* (y-1)*(y-1)
*/
/* INDENT ON */
/* INDENT OFF */
/*
* y=1
* 1 1
* A = - * atan2(2x, -x*x) = --- atan2(2,-x)
* 2 2
*
* 1 [ x*x+4] 1 4 [ 0.5(log2-logx) if
* B = - log [ -----] = - log (1+ ---) = [ |x|<E, else 0.25*
* 4 [ x*x ] 4 x*x [ log1p((2/x)*(2/x))
*/
/* INDENT ON */
if (ax < E)
else {
}
/* INDENT OFF */
/*
* Huge x:
* when |x| > 1/E^2,
* 1 pi
* A ~ --- * atan2(2x, -x*x-y*y) ~ ---
* 2 2
* y y/x
* B ~ t*(1-2t), where t = --------------- = (-------------- )/x
* x*x+(y-1)*(y-1) 1+((y-1)/x)^2
*/
/* INDENT ON */
} else if (ax < E * E * E * E) {
/* INDENT OFF */
/*
* Tiny x:
* when |x| < E^4, (note that y != 1)
* 1 1
* A = --- * atan2(2x, 1-x*x-y*y) ~ --- * atan2(2x,1-y*y)
* 2 2
*
* 1 [ (y+1)*(y+1) ] 1 2 1 2y
* B = - log [ ----------- ] = - log (1+ ---) or - log(1+ ----)
* 4 [ (y-1)*(y-1) ] 2 y-1 2 1-y
*/
/* INDENT ON */
else /* y<1 */
} else {
/* INDENT OFF */
/*
* normal x,y
* 1
* A = --- * atan2(2x, 1-x*x-y*y)
* 2
*
* 1 [ x*x+(y+1)*(y+1) ] 1 4y
* B = - log [ --------------- ] = - log (1+ -----------------)
* 4 [ x*x+(y-1)*(y-1) ] 4 x*x + (y-1)*(y-1)
*/
/* INDENT ON */
/* y close to 1 */
/* x close to 1 */
} else
}
if (hx < 0)
if (hy < 0)
return (ans);
}