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 __catan = catan 31 32 /* INDENT OFF */ 33 /* 34 * dcomplex catan(dcomplex z); 35 * 36 * If 37 * z = x + iy, 38 * 39 * then 40 * 1 ( 2x ) 1 2 2 41 * Re w = - arctan(-----------) = - ATAN2(2x, 1 - x - y ) 42 * 2 ( 2 2) 2 43 * (1 - x - y ) 44 * 45 * ( 2 2) 46 * 1 (x + (y+1) ) 1 4y 47 * Im w = - log(------------) .= --- log [ 1 + ------------- ] 48 * 4 ( 2 2) 4 2 2 49 * (x + (y-1) ) x + (y-1) 50 * 51 * 2 16 3 y 52 * = t - 2t + -- t - ..., where t = ----------------- 53 * 3 x*x + (y-1)*(y-1) 54 * 55 * Note that: if catan( x, y) = ( u, v), then 56 * catan(-x, y) = (-u, v) 57 * catan( x,-y) = ( u,-v) 58 * 59 * Also, catan(x,y) = -i*catanh(-y,x), or 60 * catanh(x,y) = i*catan(-y,x) 61 * So, if catanh(y,x) = (v,u), then catan(x,y) = -i*(-v,u) = (u,v), i.e., 62 * catan(x,y) = (u,v) 63 * 64 * EXCEPTION CASES (conform to ISO/IEC 9899:1999(E)): 65 * catan( 0 , 0 ) = (0 , 0 ) 66 * catan( NaN, 0 ) = (NaN , 0 ) 67 * catan( 0 , 1 ) = (0 , +inf) with divide-by-zero 68 * catan( inf, y ) = (pi/2 , 0 ) for finite +y 69 * catan( NaN, y ) = (NaN , NaN ) with invalid for finite y != 0 70 * catan( x , inf ) = (pi/2 , 0 ) for finite +x 71 * catan( inf, inf ) = (pi/2 , 0 ) 72 * catan( NaN, inf ) = (NaN , 0 ) 73 * catan( x , NaN ) = (NaN , NaN ) with invalid for finite x 74 * catan( inf, NaN ) = (pi/2 , +-0 ) 75 */ 76 /* INDENT ON */ 77 78 #include "libm.h" /* atan/atan2/fabs/log/log1p */ 79 #include "complex_wrapper.h" 80 81 /* INDENT OFF */ 82 static const double 83 pi_2 = 1.570796326794896558e+00, 84 zero = 0.0, 85 half = 0.5, 86 two = 2.0, 87 ln2 = 6.931471805599453094172321214581765680755e-0001, 88 one = 1.0; 89 /* INDENT ON */ 90 91 dcomplex 92 catan(dcomplex z) { 93 dcomplex ans; 94 double x, y, ax, ay, t; 95 int hx, hy, ix, iy; 96 unsigned lx, ly; 97 98 x = D_RE(z); 99 y = D_IM(z); 100 ax = fabs(x); 101 ay = fabs(y); 102 hx = HI_WORD(x); 103 lx = LO_WORD(x); 104 hy = HI_WORD(y); 105 ly = LO_WORD(y); 106 ix = hx & 0x7fffffff; 107 iy = hy & 0x7fffffff; 108 109 /* x is inf or NaN */ 110 if (ix >= 0x7ff00000) { 111 if (ISINF(ix, lx)) { 112 D_RE(ans) = pi_2; 113 D_IM(ans) = zero; 114 } else { 115 D_RE(ans) = x + x; 116 if ((iy | ly) == 0 || (ISINF(iy, ly))) 117 D_IM(ans) = zero; 118 else 119 D_IM(ans) = (fabs(y) - ay) / (fabs(y) - ay); 120 } 121 } else if (iy >= 0x7ff00000) { 122 /* y is inf or NaN */ 123 if (ISINF(iy, ly)) { 124 D_RE(ans) = pi_2; 125 D_IM(ans) = zero; 126 } else { 127 D_RE(ans) = (fabs(x) - ax) / (fabs(x) - ax); 128 D_IM(ans) = y; 129 } 130 } else if ((ix | lx) == 0) { 131 /* INDENT OFF */ 132 /* 133 * x = 0 134 * 1 1 135 * A = --- * atan2(2x, 1-x*x-y*y) = --- atan2(0,1-|y|) 136 * 2 2 137 * 138 * 1 [ (y+1)*(y+1) ] 1 2 1 2y 139 * B = - log [ ------------ ] = - log (1+ ---) or - log(1+ ----) 140 * 4 [ (y-1)*(y-1) ] 2 y-1 2 1-y 141 */ 142 /* INDENT ON */ 143 t = one - ay; 144 if (((iy - 0x3ff00000) | ly) == 0) { 145 /* y=1: catan(0,1)=(0,+inf) with 1/0 signal */ 146 D_IM(ans) = ay / ax; 147 D_RE(ans) = zero; 148 } else if (iy >= 0x3ff00000) { /* y>1 */ 149 D_IM(ans) = half * log1p(two / (-t)); 150 D_RE(ans) = pi_2; 151 } else { /* y<1 */ 152 D_IM(ans) = half * log1p((ay + ay) / t); 153 D_RE(ans) = zero; 154 } 155 } else if (iy < 0x3e200000 || ((ix - iy) >> 20) >= 30) { 156 /* INDENT OFF */ 157 /* 158 * Tiny y (relative to 1+|x|) 159 * |y| < E*(1+|x|) 160 * where E=2**-29, -35, -60 for double, double extended, quad precision 161 * 162 * 1 [ x<=1: atan(x) 163 * A = --- * atan2(2x, 1-x*x-y*y) ~ [ 1 1+x 164 * 2 [ x>=1: - atan2(2,(1-x)*(-----)) 165 * 2 x 166 * 167 * y/x 168 * B ~ t*(1-2t), where t = ----------------- is tiny 169 * x + (y-1)*(y-1)/x 170 */ 171 /* INDENT ON */ 172 if (ix < 0x3ff00000) 173 D_RE(ans) = atan(ax); 174 else 175 D_RE(ans) = half * atan2(two, (one - ax) * (one + 176 one / ax)); 177 if ((iy | ly) == 0) { 178 D_IM(ans) = ay; 179 } else { 180 if (ix < 0x3e200000) 181 t = ay / ((ay - one) * (ay - one)); 182 else if (ix > 0x41c00000) 183 t = (ay / ax) / ax; 184 else 185 t = ay / (ax * ax + (ay - one) * (ay - one)); 186 D_IM(ans) = t * (one - (t + t)); 187 } 188 } else if (iy >= 0x41c00000 && ((iy - ix) >> 20) >= 30) { 189 /* INDENT OFF */ 190 /* 191 * Huge y relative to 1+|x| 192 * |y| > Einv*(1+|x|), where Einv~2**(prec/2+3), 193 * 1 194 * A ~ --- * atan2(2x, -y*y) ~ pi/2 195 * 2 196 * y 197 * B ~ t*(1-2t), where t = --------------- is tiny 198 * (y-1)*(y-1) 199 */ 200 /* INDENT ON */ 201 D_RE(ans) = pi_2; 202 t = (ay / (ay - one)) / (ay - one); 203 D_IM(ans) = t * (one - (t + t)); 204 } else if (((iy - 0x3ff00000) | ly) == 0) { 205 /* INDENT OFF */ 206 /* 207 * y = 1 208 * 1 1 209 * A = --- * atan2(2x, -x*x) = --- atan2(2,-x) 210 * 2 2 211 * 212 * 1 [x*x + 4] 1 4 [ 0.5(log2-logx) if 213 * B = - log [-------] = - log (1+ ---) = [ |x|<E, else 0.25* 214 * 4 [ x*x ] 4 x*x [ log1p((2/x)*(2/x)) 215 */ 216 /* INDENT ON */ 217 D_RE(ans) = half * atan2(two, -ax); 218 if (ix < 0x3e200000) 219 D_IM(ans) = half * (ln2 - log(ax)); 220 else { 221 t = two / ax; 222 D_IM(ans) = 0.25 * log1p(t * t); 223 } 224 } else if (ix >= 0x43900000) { 225 /* INDENT OFF */ 226 /* 227 * Huge x: 228 * when |x| > 1/E^2, 229 * 1 pi 230 * A ~ --- * atan2(2x, -x*x-y*y) ~ --- 231 * 2 2 232 * y y/x 233 * B ~ t*(1-2t), where t = --------------- = (-------------- )/x 234 * x*x+(y-1)*(y-1) 1+((y-1)/x)^2 235 */ 236 /* INDENT ON */ 237 D_RE(ans) = pi_2; 238 t = ((ay / ax) / (one + ((ay - one) / ax) * ((ay - one) / 239 ax))) / ax; 240 D_IM(ans) = t * (one - (t + t)); 241 } else if (ix < 0x38b00000) { 242 /* INDENT OFF */ 243 /* 244 * Tiny x: 245 * when |x| < E^4, (note that y != 1) 246 * 1 1 247 * A = --- * atan2(2x, 1-x*x-y*y) ~ --- * atan2(2x,(1-y)*(1+y)) 248 * 2 2 249 * 250 * 1 [(y+1)*(y+1)] 1 2 1 2y 251 * B = - log [-----------] = - log (1+ ---) or - log(1+ ----) 252 * 4 [(y-1)*(y-1)] 2 y-1 2 1-y 253 */ 254 /* INDENT ON */ 255 D_RE(ans) = half * atan2(ax + ax, (one - ay) * (one + ay)); 256 if (iy >= 0x3ff00000) 257 D_IM(ans) = half * log1p(two / (ay - one)); 258 else 259 D_IM(ans) = half * log1p((ay + ay) / (one - ay)); 260 } else { 261 /* INDENT OFF */ 262 /* 263 * normal x,y 264 * 1 265 * A = --- * atan2(2x, 1-x*x-y*y) 266 * 2 267 * 268 * 1 [x*x+(y+1)*(y+1)] 1 4y 269 * B = - log [---------------] = - log (1+ -----------------) 270 * 4 [x*x+(y-1)*(y-1)] 4 x*x + (y-1)*(y-1) 271 */ 272 /* INDENT ON */ 273 t = one - ay; 274 if (iy >= 0x3fe00000 && iy < 0x40000000) { 275 /* y close to 1 */ 276 D_RE(ans) = half * (atan2((ax + ax), (t * (one + ay) - 277 ax * ax))); 278 } else if (ix >= 0x3fe00000 && ix < 0x40000000) { 279 /* x close to 1 */ 280 D_RE(ans) = half * atan2((ax + ax), ((one - ax) * 281 (one + ax) - ay * ay)); 282 } else 283 D_RE(ans) = half * atan2((ax + ax), ((one - ax * ax) - 284 ay * ay)); 285 D_IM(ans) = 0.25 * log1p((4.0 * ay) / (ax * ax + t * t)); 286 } 287 if (hx < 0) 288 D_RE(ans) = -D_RE(ans); 289 if (hy < 0) 290 D_IM(ans) = -D_IM(ans); 291 return (ans); 292 } 293