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