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 __casin = casin 31 32 /* INDENT OFF */ 33 /* 34 * dcomplex casin(dcomplex z); 35 * 36 * Alogrithm 37 * (based on T.E.Hull, Thomas F. Fairgrieve and Ping Tak Peter Tang's 38 * paper "Implementing the Complex Arcsine and Arccosine Functins Using 39 * Exception Handling", ACM TOMS, Vol 23, pp 299-335) 40 * 41 * The principal value of complex inverse sine function casin(z), 42 * where z = x+iy, can be defined by 43 * 44 * casin(z) = asin(B) + i sign(y) log (A + sqrt(A*A-1)), 45 * 46 * where the log function is the natural log, and 47 * ____________ ____________ 48 * 1 / 2 2 1 / 2 2 49 * A = --- / (x+1) + y + --- / (x-1) + y 50 * 2 \/ 2 \/ 51 * ____________ ____________ 52 * 1 / 2 2 1 / 2 2 53 * B = --- / (x+1) + y - --- / (x-1) + y . 54 * 2 \/ 2 \/ 55 * 56 * The Branch cuts are on the real line from -inf to -1 and from 1 to inf. 57 * The real and imaginary parts are based on Abramowitz and Stegun 58 * [Handbook of Mathematic Functions, 1972]. The sign of the imaginary 59 * part is chosen to be the generally considered the principal value of 60 * this function. 61 * 62 * Notes:1. A is the average of the distances from z to the points (1,0) 63 * and (-1,0) in the complex z-plane, and in particular A>=1. 64 * 2. B is in [-1,1], and A*B = x. 65 * 66 * Special notes: if casin( x, y) = ( u, v), then 67 * casin(-x, y) = (-u, v), 68 * casin( x,-y) = ( u,-v), 69 * in general, we have casin(conj(z)) = conj(casin(z)) 70 * casin(-z) = -casin(z) 71 * casin(z) = pi/2 - cacos(z) 72 * 73 * EXCEPTION CASES (conform to ISO/IEC 9899:1999(E)): 74 * casin( 0 + i 0 ) = 0 + i 0 75 * casin( 0 + i NaN ) = 0 + i NaN 76 * casin( x + i inf ) = 0 + i inf for finite x 77 * casin( x + i NaN ) = NaN + i NaN with invalid for finite x != 0 78 * casin(inf + iy ) = pi/2 + i inf finite y 79 * casin(inf + i inf) = pi/4 + i inf 80 * casin(inf + i NaN) = NaN + i inf 81 * casin(NaN + i y ) = NaN + i NaN for finite y 82 * casin(NaN + i inf) = NaN + i inf 83 * casin(NaN + i NaN) = NaN + i NaN 84 * 85 * Special Regions (better formula for accuracy and for avoiding spurious 86 * overflow or underflow) (all x and y are assumed nonnegative): 87 * case 1: y = 0 88 * case 2: tiny y relative to x-1: y <= ulp(0.5)*|x-1| 89 * case 3: tiny y: y < 4 sqrt(u), where u = minimum normal number 90 * case 4: huge y relative to x+1: y >= (1+x)/ulp(0.5) 91 * case 5: huge x and y: x and y >= sqrt(M)/8, where M = maximum normal number 92 * case 6: tiny x: x < 4 sqrt(u) 93 * -------- 94 * case 1 & 2. y=0 or y/|x-1| is tiny. We have 95 * ____________ _____________ 96 * / 2 2 / y 2 97 * / (x+-1) + y = |x+-1| / 1 + (------) 98 * \/ \/ |x+-1| 99 * 100 * 1 y 2 101 * ~ |x+-1| ( 1 + --- (------) ) 102 * 2 |x+-1| 103 * 104 * 2 105 * y 106 * = |x+-1| + --------. 107 * 2|x+-1| 108 * 109 * Consequently, it is not difficult to see that 110 * 2 111 * y 112 * [ 1 + ------------ , if x < 1, 113 * [ 2(1+x)(1-x) 114 * [ 115 * [ 116 * [ x, if x = 1 (y = 0), 117 * [ 118 * A ~= [ 2 119 * [ x * y 120 * [ x + ------------ , if x > 1 121 * [ 2(1+x)(x-1) 122 * 123 * and hence 124 * ______ 2 125 * / 2 y y 126 * A + \/ A - 1 ~ 1 + ---------------- + -----------, if x < 1, 127 * sqrt((x+1)(1-x)) 2(x+1)(1-x) 128 * 129 * 130 * ~ x + sqrt((x-1)*(x+1)), if x >= 1. 131 * 132 * 2 133 * y 134 * [ x(1 - ------------), if x < 1, 135 * [ 2(1+x)(1-x) 136 * B = x/A ~ [ 137 * [ 1, if x = 1, 138 * [ 139 * [ 2 140 * [ y 141 * [ 1 - ------------ , if x > 1, 142 * [ 2(1+x)(1-x) 143 * Thus 144 * [ asin(x) + i y/sqrt((x-1)*(x+1)), if x < 1 145 * casin(x+i*y)=[ 146 * [ pi/2 + i log(x+sqrt(x*x-1)), if x >= 1 147 * 148 * case 3. y < 4 sqrt(u), where u = minimum normal x. 149 * After case 1 and 2, this will only occurs when x=1. When x=1, we have 150 * A = (sqrt(4+y*y)+y)/2 ~ 1 + y/2 + y^2/8 + ... 151 * and 152 * B = 1/A = 1 - y/2 + y^2/8 + ... 153 * Since 154 * asin(x) = pi/2-2*asin(sqrt((1-x)/2)) 155 * asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ... 156 * we have, for the real part asin(B), 157 * asin(1-y/2) ~ pi/2 - 2 asin(sqrt(y/4)) 158 * ~ pi/2 - sqrt(y) 159 * For the imaginary part, 160 * log(A+sqrt(A*A-1)) ~ log(1+y/2+sqrt(2*y/2)) 161 * = log(1+y/2+sqrt(y)) 162 * = (y/2+sqrt(y)) - (y/2+sqrt(y))^2/2 + ... 163 * ~ sqrt(y) - y*(sqrt(y)+y/2)/2 164 * ~ sqrt(y) 165 * 166 * case 4. y >= (x+1)ulp(0.5). In this case, A ~ y and B ~ x/y. Thus 167 * real part = asin(B) ~ x/y (be careful, x/y may underflow) 168 * and 169 * imag part = log(y+sqrt(y*y-one)) 170 * 171 * 172 * case 5. Both x and y are large: x and y > sqrt(M)/8, where M = maximum x 173 * In this case, 174 * A ~ sqrt(x*x+y*y) 175 * B ~ x/sqrt(x*x+y*y). 176 * Thus 177 * real part = asin(B) = atan(x/y), 178 * imag part = log(A+sqrt(A*A-1)) ~ log(2A) 179 * = log(2) + 0.5*log(x*x+y*y) 180 * = log(2) + log(y) + 0.5*log(1+(x/y)^2) 181 * 182 * case 6. x < 4 sqrt(u). In this case, we have 183 * A ~ sqrt(1+y*y), B = x/sqrt(1+y*y). 184 * Since B is tiny, we have 185 * real part = asin(B) ~ B = x/sqrt(1+y*y) 186 * imag part = log(A+sqrt(A*A-1)) = log (A+sqrt(y*y)) 187 * = log(y+sqrt(1+y*y)) 188 * = 0.5*log(y^2+2ysqrt(1+y^2)+1+y^2) 189 * = 0.5*log(1+2y(y+sqrt(1+y^2))); 190 * = 0.5*log1p(2y(y+A)); 191 * 192 * casin(z) = asin(B) + i sign(y) log (A + sqrt(A*A-1)), 193 */ 194 /* INDENT ON */ 195 196 #include "libm.h" /* asin/atan/fabs/log/log1p/sqrt */ 197 #include "complex_wrapper.h" 198 199 /* INDENT OFF */ 200 static const double 201 zero = 0.0, 202 one = 1.0, 203 E = 1.11022302462515654042e-16, /* 2**-53 */ 204 ln2 = 6.93147180559945286227e-01, 205 pi_2 = 1.570796326794896558e+00, 206 pi_2_l = 6.123233995736765886e-17, 207 pi_4 = 7.85398163397448278999e-01, 208 Foursqrtu = 5.96667258496016539463e-154, /* 2**(-509) */ 209 Acrossover = 1.5, 210 Bcrossover = 0.6417, 211 half = 0.5; 212 /* INDENT ON */ 213 214 dcomplex 215 casin(dcomplex z) { 216 double x, y, t, R, S, A, Am1, B, y2, xm1, xp1, Apx; 217 int ix, iy, hx, hy; 218 unsigned lx, ly; 219 dcomplex ans; 220 221 x = D_RE(z); 222 y = D_IM(z); 223 hx = HI_WORD(x); 224 lx = LO_WORD(x); 225 hy = HI_WORD(y); 226 ly = LO_WORD(y); 227 ix = hx & 0x7fffffff; 228 iy = hy & 0x7fffffff; 229 x = fabs(x); 230 y = fabs(y); 231 232 /* special cases */ 233 234 /* x is inf or NaN */ 235 if (ix >= 0x7ff00000) { /* x is inf or NaN */ 236 if (ISINF(ix, lx)) { /* x is INF */ 237 D_IM(ans) = x; 238 if (iy >= 0x7ff00000) { 239 if (ISINF(iy, ly)) 240 /* casin(inf + i inf) = pi/4 + i inf */ 241 D_RE(ans) = pi_4; 242 else /* casin(inf + i NaN) = NaN + i inf */ 243 D_RE(ans) = y + y; 244 } else /* casin(inf + iy) = pi/2 + i inf */ 245 D_RE(ans) = pi_2; 246 } else { /* x is NaN */ 247 if (iy >= 0x7ff00000) { 248 /* INDENT OFF */ 249 /* 250 * casin(NaN + i inf) = NaN + i inf 251 * casin(NaN + i NaN) = NaN + i NaN 252 */ 253 /* INDENT ON */ 254 D_IM(ans) = y + y; 255 D_RE(ans) = x + x; 256 } else { 257 /* casin(NaN + i y ) = NaN + i NaN */ 258 D_IM(ans) = D_RE(ans) = x + y; 259 } 260 } 261 if (hx < 0) 262 D_RE(ans) = -D_RE(ans); 263 if (hy < 0) 264 D_IM(ans) = -D_IM(ans); 265 return (ans); 266 } 267 268 /* casin(+0 + i 0 ) = 0 + i 0. */ 269 if ((ix | lx | iy | ly) == 0) 270 return (z); 271 272 if (iy >= 0x7ff00000) { /* y is inf or NaN */ 273 if (ISINF(iy, ly)) { /* casin(x + i inf) = 0 + i inf */ 274 D_IM(ans) = y; 275 D_RE(ans) = zero; 276 } else { /* casin(x + i NaN) = NaN + i NaN */ 277 D_IM(ans) = x + y; 278 if ((ix | lx) == 0) 279 D_RE(ans) = x; 280 else 281 D_RE(ans) = y; 282 } 283 if (hx < 0) 284 D_RE(ans) = -D_RE(ans); 285 if (hy < 0) 286 D_IM(ans) = -D_IM(ans); 287 return (ans); 288 } 289 290 if ((iy | ly) == 0) { /* region 1: y=0 */ 291 if (ix < 0x3ff00000) { /* |x| < 1 */ 292 D_RE(ans) = asin(x); 293 D_IM(ans) = zero; 294 } else { 295 D_RE(ans) = pi_2; 296 if (ix >= 0x43500000) /* |x| >= 2**54 */ 297 D_IM(ans) = ln2 + log(x); 298 else if (ix >= 0x3ff80000) /* x > Acrossover */ 299 D_IM(ans) = log(x + sqrt((x - one) * (x + 300 one))); 301 else { 302 xm1 = x - one; 303 D_IM(ans) = log1p(xm1 + sqrt(xm1 * (x + one))); 304 } 305 } 306 } else if (y <= E * fabs(x - one)) { /* region 2: y < tiny*|x-1| */ 307 if (ix < 0x3ff00000) { /* x < 1 */ 308 D_RE(ans) = asin(x); 309 D_IM(ans) = y / sqrt((one + x) * (one - x)); 310 } else { 311 D_RE(ans) = pi_2; 312 if (ix >= 0x43500000) { /* |x| >= 2**54 */ 313 D_IM(ans) = ln2 + log(x); 314 } else if (ix >= 0x3ff80000) /* x > Acrossover */ 315 D_IM(ans) = log(x + sqrt((x - one) * (x + 316 one))); 317 else 318 D_IM(ans) = log1p((x - one) + sqrt((x - one) * 319 (x + one))); 320 } 321 } else if (y < Foursqrtu) { /* region 3 */ 322 t = sqrt(y); 323 D_RE(ans) = pi_2 - (t - pi_2_l); 324 D_IM(ans) = t; 325 } else if (E * y - one >= x) { /* region 4 */ 326 D_RE(ans) = x / y; /* need to fix underflow cases */ 327 D_IM(ans) = ln2 + log(y); 328 } else if (ix >= 0x5fc00000 || iy >= 0x5fc00000) { /* x,y>2**509 */ 329 /* region 5: x+1 or y is very large (>= sqrt(max)/8) */ 330 t = x / y; 331 D_RE(ans) = atan(t); 332 D_IM(ans) = ln2 + log(y) + half * log1p(t * t); 333 } else if (x < Foursqrtu) { 334 /* region 6: x is very small, < 4sqrt(min) */ 335 A = sqrt(one + y * y); 336 D_RE(ans) = x / A; /* may underflow */ 337 if (iy >= 0x3ff80000) /* if y > Acrossover */ 338 D_IM(ans) = log(y + A); 339 else 340 D_IM(ans) = half * log1p((y + y) * (y + A)); 341 } else { /* safe region */ 342 y2 = y * y; 343 xp1 = x + one; 344 xm1 = x - one; 345 R = sqrt(xp1 * xp1 + y2); 346 S = sqrt(xm1 * xm1 + y2); 347 A = half * (R + S); 348 B = x / A; 349 350 if (B <= Bcrossover) 351 D_RE(ans) = asin(B); 352 else { /* use atan and an accurate approx to a-x */ 353 Apx = A + x; 354 if (x <= one) 355 D_RE(ans) = atan(x / sqrt(half * Apx * (y2 / 356 (R + xp1) + (S - xm1)))); 357 else 358 D_RE(ans) = atan(x / (y * sqrt(half * (Apx / 359 (R + xp1) + Apx / (S + xm1))))); 360 } 361 if (A <= Acrossover) { 362 /* use log1p and an accurate approx to A-1 */ 363 if (x < one) 364 Am1 = half * (y2 / (R + xp1) + y2 / (S - xm1)); 365 else 366 Am1 = half * (y2 / (R + xp1) + (S + xm1)); 367 D_IM(ans) = log1p(Am1 + sqrt(Am1 * (A + one))); 368 } else { 369 D_IM(ans) = log(A + sqrt(A * A - one)); 370 } 371 } 372 373 if (hx < 0) 374 D_RE(ans) = -D_RE(ans); 375 if (hy < 0) 376 D_IM(ans) = -D_IM(ans); 377 378 return (ans); 379 } 380