25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER START
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * The contents of this file are subject to the terms of the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Common Development and Distribution License (the "License").
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You may not use this file except in compliance with the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * See the License for the specific language governing permissions
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * and limitations under the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER END
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Use is subject to license terms.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* INDENT OFF */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * dcomplex catan(dcomplex z);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 1 ( 2x ) 1 2 2
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Re w = - arctan(-----------) = - ATAN2(2x, 1 - x - y )
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (1 - x - y )
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 1 (x + (y+1) ) 1 4y
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Im w = - log(------------) .= --- log [ 1 + ------------- ]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 4 ( 2 2) 4 2 2
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (x + (y-1) ) x + (y-1)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * = t - 2t + -- t - ..., where t = -----------------
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 3 x*x + (y-1)*(y-1)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Note that: if catan( x, y) = ( u, v), then
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * catan(-x, y) = (-u, v)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * catan( x,-y) = ( u,-v)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Also, catan(x,y) = -i*catanh(-y,x), or
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * catanh(x,y) = i*catan(-y,x)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * So, if catanh(y,x) = (v,u), then catan(x,y) = -i*(-v,u) = (u,v), i.e.,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * catan(x,y) = (u,v)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * EXCEPTION CASES (conform to ISO/IEC 9899:1999(E)):
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * catan( 0 , 0 ) = (0 , 0 )
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * catan( NaN, 0 ) = (NaN , 0 )
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * catan( 0 , 1 ) = (0 , +inf) with divide-by-zero
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * catan( inf, y ) = (pi/2 , 0 ) for finite +y
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * catan( NaN, y ) = (NaN , NaN ) with invalid for finite y != 0
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * catan( x , inf ) = (pi/2 , 0 ) for finite +x
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * catan( inf, inf ) = (pi/2 , 0 )
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * catan( NaN, inf ) = (NaN , 0 )
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * catan( x , NaN ) = (NaN , NaN ) with invalid for finite x
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * catan( inf, NaN ) = (pi/2 , +-0 )
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* INDENT ON */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#include "libm.h" /* atan/atan2/fabs/log/log1p */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* INDENT OFF */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ln2 = 6.931471805599453094172321214581765680755e-0001,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* INDENT ON */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* x is inf or NaN */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis D_IM(ans) = (fabs(y) - ay) / (fabs(y) - ay);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* y is inf or NaN */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis D_RE(ans) = (fabs(x) - ax) / (fabs(x) - ax);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* INDENT OFF */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * A = --- * atan2(2x, 1-x*x-y*y) = --- atan2(0,1-|y|)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 1 [ (y+1)*(y+1) ] 1 2 1 2y
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * B = - log [ ------------ ] = - log (1+ ---) or - log(1+ ----)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 4 [ (y-1)*(y-1) ] 2 y-1 2 1-y
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* INDENT ON */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* y=1: catan(0,1)=(0,+inf) with 1/0 signal */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else { /* y<1 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (iy < 0x3e200000 || ((ix - iy) >> 20) >= 30) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* INDENT OFF */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Tiny y (relative to 1+|x|)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * |y| < E*(1+|x|)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * where E=2**-29, -35, -60 for double, double extended, quad precision
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 1 [ x<=1: atan(x)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * A = --- * atan2(2x, 1-x*x-y*y) ~ [ 1 1+x
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 2 [ x>=1: - atan2(2,(1-x)*(-----))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * B ~ t*(1-2t), where t = ----------------- is tiny
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * x + (y-1)*(y-1)/x
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* INDENT ON */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis D_RE(ans) = half * atan2(two, (one - ax) * (one +
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t = ay / (ax * ax + (ay - one) * (ay - one));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (iy >= 0x41c00000 && ((iy - ix) >> 20) >= 30) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* INDENT OFF */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Huge y relative to 1+|x|
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * |y| > Einv*(1+|x|), where Einv~2**(prec/2+3),
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * A ~ --- * atan2(2x, -y*y) ~ pi/2
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * B ~ t*(1-2t), where t = --------------- is tiny
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* INDENT ON */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (((iy - 0x3ff00000) | ly) == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* INDENT OFF */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * A = --- * atan2(2x, -x*x) = --- atan2(2,-x)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 1 [x*x + 4] 1 4 [ 0.5(log2-logx) if
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * B = - log [-------] = - log (1+ ---) = [ |x|<E, else 0.25*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 4 [ x*x ] 4 x*x [ log1p((2/x)*(2/x))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* INDENT ON */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* INDENT OFF */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * when |x| > 1/E^2,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * A ~ --- * atan2(2x, -x*x-y*y) ~ ---
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * B ~ t*(1-2t), where t = --------------- = (-------------- )/x
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * x*x+(y-1)*(y-1) 1+((y-1)/x)^2
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* INDENT ON */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t = ((ay / ax) / (one + ((ay - one) / ax) * ((ay - one) /
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* INDENT OFF */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * when |x| < E^4, (note that y != 1)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * A = --- * atan2(2x, 1-x*x-y*y) ~ --- * atan2(2x,(1-y)*(1+y))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 1 [(y+1)*(y+1)] 1 2 1 2y
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * B = - log [-----------] = - log (1+ ---) or - log(1+ ----)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 4 [(y-1)*(y-1)] 2 y-1 2 1-y
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* INDENT ON */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis D_RE(ans) = half * atan2(ax + ax, (one - ay) * (one + ay));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis D_IM(ans) = half * log1p(two / (ay - one));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis D_IM(ans) = half * log1p((ay + ay) / (one - ay));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* INDENT OFF */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * A = --- * atan2(2x, 1-x*x-y*y)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 1 [x*x+(y+1)*(y+1)] 1 4y
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * B = - log [---------------] = - log (1+ -----------------)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 4 [x*x+(y-1)*(y-1)] 4 x*x + (y-1)*(y-1)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* INDENT ON */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (iy >= 0x3fe00000 && iy < 0x40000000) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* y close to 1 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis D_RE(ans) = half * (atan2((ax + ax), (t * (one + ay) -
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (ix >= 0x3fe00000 && ix < 0x40000000) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* x close to 1 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis D_RE(ans) = half * atan2((ax + ax), ((one - ax) *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis D_RE(ans) = half * atan2((ax + ax), ((one - ax * ax) -
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis D_IM(ans) = 0.25 * log1p((4.0 * ay) / (ax * ax + t * t));