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 */ 109 static const long double 110 zero = 0.0L, 111 one = 1.0L, 112 two = 2.0L, 113 half = 0.5L, 114 ln2 = 6.931471805599453094172321214581765680755e-0001L, 115 pi_2 = 1.570796326794896619231321691639751442098584699687552910487472L, 116 #if defined(__x86) 117 E = 2.910383045673370361328125000000000000000e-11L, /* 2**-35 */ 118 Einv = 3.435973836800000000000000000000000000000e+10L; /* 2**+35 */ 119 #else 120 E = 8.673617379884035472059622406959533691406e-19L, /* 2**-60 */ 121 Einv = 1.152921504606846976000000000000000000000e18L; /* 2**+60 */ 122 #endif 123 /* INDENT ON */ 124 125 ldcomplex 126 catanl(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