1 /*- 2 * Copyright (c) 2012 Stephen Montgomery-Smith <stephen@FreeBSD.ORG> 3 * All rights reserved. 4 * 5 * Redistribution and use in source and binary forms, with or without 6 * modification, are permitted provided that the following conditions 7 * are met: 8 * 1. Redistributions of source code must retain the above copyright 9 * notice, this list of conditions and the following disclaimer. 10 * 2. Redistributions in binary form must reproduce the above copyright 11 * notice, this list of conditions and the following disclaimer in the 12 * documentation and/or other materials provided with the distribution. 13 * 14 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 15 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 16 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 17 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 18 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 19 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 20 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 21 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 22 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 23 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 24 * SUCH DAMAGE. 25 */ 26 27 #include <sys/cdefs.h> 28 __FBSDID("$FreeBSD$"); 29 30 #include <complex.h> 31 #include <float.h> 32 33 #include "math.h" 34 #include "math_private.h" 35 36 #undef isinf 37 #define isinf(x) (fabs(x) == INFINITY) 38 #undef isnan 39 #define isnan(x) ((x) != (x)) 40 #define raise_inexact() do { volatile float junk = 1 + tiny; } while(0) 41 #undef signbit 42 #define signbit(x) (__builtin_signbit(x)) 43 44 /* We need that DBL_EPSILON^2/128 is larger than FOUR_SQRT_MIN. */ 45 static const double 46 A_crossover = 10, /* Hull et al suggest 1.5, but 10 works better */ 47 B_crossover = 0.6417, /* suggested by Hull et al */ 48 FOUR_SQRT_MIN = 0x1p-509, /* >= 4 * sqrt(DBL_MIN) */ 49 QUARTER_SQRT_MAX = 0x1p509, /* <= sqrt(DBL_MAX) / 4 */ 50 m_e = 2.7182818284590452e0, /* 0x15bf0a8b145769.0p-51 */ 51 m_ln2 = 6.9314718055994531e-1, /* 0x162e42fefa39ef.0p-53 */ 52 pio2_hi = 1.5707963267948966e0, /* 0x1921fb54442d18.0p-52 */ 53 RECIP_EPSILON = 1 / DBL_EPSILON, 54 SQRT_3_EPSILON = 2.5809568279517849e-8, /* 0x1bb67ae8584caa.0p-78 */ 55 SQRT_6_EPSILON = 3.6500241499888571e-8, /* 0x13988e1409212e.0p-77 */ 56 SQRT_MIN = 0x1p-511; /* >= sqrt(DBL_MIN) */ 57 58 static const volatile double 59 pio2_lo = 6.1232339957367659e-17; /* 0x11a62633145c07.0p-106 */ 60 static const volatile float 61 tiny = 0x1p-100; 62 63 static double complex clog_for_large_values(double complex z); 64 65 /* 66 * Testing indicates that all these functions are accurate up to 4 ULP. 67 * The functions casin(h) and cacos(h) are about 2.5 times slower than asinh. 68 * The functions catan(h) are a little under 2 times slower than atanh. 69 * 70 * The code for casinh, casin, cacos, and cacosh comes first. The code is 71 * rather complicated, and the four functions are highly interdependent. 72 * 73 * The code for catanh and catan comes at the end. It is much simpler than 74 * the other functions, and the code for these can be disconnected from the 75 * rest of the code. 76 */ 77 78 /* 79 * ================================ 80 * | casinh, casin, cacos, cacosh | 81 * ================================ 82 */ 83 84 /* 85 * The algorithm is very close to that in "Implementing the complex arcsine 86 * and arccosine functions using exception handling" by T. E. Hull, Thomas F. 87 * Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on 88 * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335, 89 * http://dl.acm.org/citation.cfm?id=275324. 90 * 91 * Throughout we use the convention z = x + I*y. 92 * 93 * casinh(z) = sign(x)*log(A+sqrt(A*A-1)) + I*asin(B) 94 * where 95 * A = (|z+I| + |z-I|) / 2 96 * B = (|z+I| - |z-I|) / 2 = y/A 97 * 98 * These formulas become numerically unstable: 99 * (a) for Re(casinh(z)) when z is close to the line segment [-I, I] (that 100 * is, Re(casinh(z)) is close to 0); 101 * (b) for Im(casinh(z)) when z is close to either of the intervals 102 * [I, I*infinity) or (-I*infinity, -I] (that is, |Im(casinh(z))| is 103 * close to PI/2). 104 * 105 * These numerical problems are overcome by defining 106 * f(a, b) = (hypot(a, b) - b) / 2 = a*a / (hypot(a, b) + b) / 2 107 * Then if A < A_crossover, we use 108 * log(A + sqrt(A*A-1)) = log1p((A-1) + sqrt((A-1)*(A+1))) 109 * A-1 = f(x, 1+y) + f(x, 1-y) 110 * and if B > B_crossover, we use 111 * asin(B) = atan2(y, sqrt(A*A - y*y)) = atan2(y, sqrt((A+y)*(A-y))) 112 * A-y = f(x, y+1) + f(x, y-1) 113 * where without loss of generality we have assumed that x and y are 114 * non-negative. 115 * 116 * Much of the difficulty comes because the intermediate computations may 117 * produce overflows or underflows. This is dealt with in the paper by Hull 118 * et al by using exception handling. We do this by detecting when 119 * computations risk underflow or overflow. The hardest part is handling the 120 * underflows when computing f(a, b). 121 * 122 * Note that the function f(a, b) does not appear explicitly in the paper by 123 * Hull et al, but the idea may be found on pages 308 and 309. Introducing the 124 * function f(a, b) allows us to concentrate many of the clever tricks in this 125 * paper into one function. 126 */ 127 128 /* 129 * Function f(a, b, hypot_a_b) = (hypot(a, b) - b) / 2. 130 * Pass hypot(a, b) as the third argument. 131 */ 132 static inline double 133 f(double a, double b, double hypot_a_b) 134 { 135 if (b < 0) 136 return ((hypot_a_b - b) / 2); 137 if (b == 0) 138 return (a / 2); 139 return (a * a / (hypot_a_b + b) / 2); 140 } 141 142 /* 143 * All the hard work is contained in this function. 144 * x and y are assumed positive or zero, and less than RECIP_EPSILON. 145 * Upon return: 146 * rx = Re(casinh(z)) = -Im(cacos(y + I*x)). 147 * B_is_usable is set to 1 if the value of B is usable. 148 * If B_is_usable is set to 0, sqrt_A2my2 = sqrt(A*A - y*y), and new_y = y. 149 * If returning sqrt_A2my2 has potential to result in an underflow, it is 150 * rescaled, and new_y is similarly rescaled. 151 */ 152 static inline void 153 do_hard_work(double x, double y, double *rx, int *B_is_usable, double *B, 154 double *sqrt_A2my2, double *new_y) 155 { 156 double R, S, A; /* A, B, R, and S are as in Hull et al. */ 157 double Am1, Amy; /* A-1, A-y. */ 158 159 R = hypot(x, y + 1); /* |z+I| */ 160 S = hypot(x, y - 1); /* |z-I| */ 161 162 /* A = (|z+I| + |z-I|) / 2 */ 163 A = (R + S) / 2; 164 /* 165 * Mathematically A >= 1. There is a small chance that this will not 166 * be so because of rounding errors. So we will make certain it is 167 * so. 168 */ 169 if (A < 1) 170 A = 1; 171 172 if (A < A_crossover) { 173 /* 174 * Am1 = fp + fm, where fp = f(x, 1+y), and fm = f(x, 1-y). 175 * rx = log1p(Am1 + sqrt(Am1*(A+1))) 176 */ 177 if (y == 1 && x < DBL_EPSILON * DBL_EPSILON / 128) { 178 /* 179 * fp is of order x^2, and fm = x/2. 180 * A = 1 (inexactly). 181 */ 182 *rx = sqrt(x); 183 } else if (x >= DBL_EPSILON * fabs(y - 1)) { 184 /* 185 * Underflow will not occur because 186 * x >= DBL_EPSILON^2/128 >= FOUR_SQRT_MIN 187 */ 188 Am1 = f(x, 1 + y, R) + f(x, 1 - y, S); 189 *rx = log1p(Am1 + sqrt(Am1 * (A + 1))); 190 } else if (y < 1) { 191 /* 192 * fp = x*x/(1+y)/4, fm = x*x/(1-y)/4, and 193 * A = 1 (inexactly). 194 */ 195 *rx = x / sqrt((1 - y) * (1 + y)); 196 } else { /* if (y > 1) */ 197 /* 198 * A-1 = y-1 (inexactly). 199 */ 200 *rx = log1p((y - 1) + sqrt((y - 1) * (y + 1))); 201 } 202 } else { 203 *rx = log(A + sqrt(A * A - 1)); 204 } 205 206 *new_y = y; 207 208 if (y < FOUR_SQRT_MIN) { 209 /* 210 * Avoid a possible underflow caused by y/A. For casinh this 211 * would be legitimate, but will be picked up by invoking atan2 212 * later on. For cacos this would not be legitimate. 213 */ 214 *B_is_usable = 0; 215 *sqrt_A2my2 = A * (2 / DBL_EPSILON); 216 *new_y = y * (2 / DBL_EPSILON); 217 return; 218 } 219 220 /* B = (|z+I| - |z-I|) / 2 = y/A */ 221 *B = y / A; 222 *B_is_usable = 1; 223 224 if (*B > B_crossover) { 225 *B_is_usable = 0; 226 /* 227 * Amy = fp + fm, where fp = f(x, y+1), and fm = f(x, y-1). 228 * sqrt_A2my2 = sqrt(Amy*(A+y)) 229 */ 230 if (y == 1 && x < DBL_EPSILON / 128) { 231 /* 232 * fp is of order x^2, and fm = x/2. 233 * A = 1 (inexactly). 234 */ 235 *sqrt_A2my2 = sqrt(x) * sqrt((A + y) / 2); 236 } else if (x >= DBL_EPSILON * fabs(y - 1)) { 237 /* 238 * Underflow will not occur because 239 * x >= DBL_EPSILON/128 >= FOUR_SQRT_MIN 240 * and 241 * x >= DBL_EPSILON^2 >= FOUR_SQRT_MIN 242 */ 243 Amy = f(x, y + 1, R) + f(x, y - 1, S); 244 *sqrt_A2my2 = sqrt(Amy * (A + y)); 245 } else if (y > 1) { 246 /* 247 * fp = x*x/(y+1)/4, fm = x*x/(y-1)/4, and 248 * A = y (inexactly). 249 * 250 * y < RECIP_EPSILON. So the following 251 * scaling should avoid any underflow problems. 252 */ 253 *sqrt_A2my2 = x * (4 / DBL_EPSILON / DBL_EPSILON) * y / 254 sqrt((y + 1) * (y - 1)); 255 *new_y = y * (4 / DBL_EPSILON / DBL_EPSILON); 256 } else { /* if (y < 1) */ 257 /* 258 * fm = 1-y >= DBL_EPSILON, fp is of order x^2, and 259 * A = 1 (inexactly). 260 */ 261 *sqrt_A2my2 = sqrt((1 - y) * (1 + y)); 262 } 263 } 264 } 265 266 /* 267 * casinh(z) = z + O(z^3) as z -> 0 268 * 269 * casinh(z) = sign(x)*clog(sign(x)*z) + O(1/z^2) as z -> infinity 270 * The above formula works for the imaginary part as well, because 271 * Im(casinh(z)) = sign(x)*atan2(sign(x)*y, fabs(x)) + O(y/z^3) 272 * as z -> infinity, uniformly in y 273 */ 274 double complex 275 casinh(double complex z) 276 { 277 double x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y; 278 int B_is_usable; 279 double complex w; 280 281 x = creal(z); 282 y = cimag(z); 283 ax = fabs(x); 284 ay = fabs(y); 285 286 if (isnan(x) || isnan(y)) { 287 /* casinh(+-Inf + I*NaN) = +-Inf + I*NaN */ 288 if (isinf(x)) 289 return (cpack(x, y + y)); 290 /* casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */ 291 if (isinf(y)) 292 return (cpack(y, x + x)); 293 /* casinh(NaN + I*0) = NaN + I*0 */ 294 if (y == 0) 295 return (cpack(x + x, y)); 296 /* 297 * All other cases involving NaN return NaN + I*NaN. 298 * C99 leaves it optional whether to raise invalid if one of 299 * the arguments is not NaN, so we opt not to raise it. 300 */ 301 return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); 302 } 303 304 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { 305 /* clog...() will raise inexact unless x or y is infinite. */ 306 if (signbit(x) == 0) 307 w = clog_for_large_values(z) + m_ln2; 308 else 309 w = clog_for_large_values(-z) + m_ln2; 310 return (cpack(copysign(creal(w), x), copysign(cimag(w), y))); 311 } 312 313 /* Avoid spuriously raising inexact for z = 0. */ 314 if (x == 0 && y == 0) 315 return (z); 316 317 /* All remaining cases are inexact. */ 318 raise_inexact(); 319 320 if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) 321 return (z); 322 323 do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y); 324 if (B_is_usable) 325 ry = asin(B); 326 else 327 ry = atan2(new_y, sqrt_A2my2); 328 return (cpack(copysign(rx, x), copysign(ry, y))); 329 } 330 331 /* 332 * casin(z) = reverse(casinh(reverse(z))) 333 * where reverse(x + I*y) = y + I*x = I*conj(z). 334 */ 335 double complex 336 casin(double complex z) 337 { 338 double complex w = casinh(cpack(cimag(z), creal(z))); 339 340 return (cpack(cimag(w), creal(w))); 341 } 342 343 /* 344 * cacos(z) = PI/2 - casin(z) 345 * but do the computation carefully so cacos(z) is accurate when z is 346 * close to 1. 347 * 348 * cacos(z) = PI/2 - z + O(z^3) as z -> 0 349 * 350 * cacos(z) = -sign(y)*I*clog(z) + O(1/z^2) as z -> infinity 351 * The above formula works for the real part as well, because 352 * Re(cacos(z)) = atan2(fabs(y), x) + O(y/z^3) 353 * as z -> infinity, uniformly in y 354 */ 355 double complex 356 cacos(double complex z) 357 { 358 double x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x; 359 int sx, sy; 360 int B_is_usable; 361 double complex w; 362 363 x = creal(z); 364 y = cimag(z); 365 sx = signbit(x); 366 sy = signbit(y); 367 ax = fabs(x); 368 ay = fabs(y); 369 370 if (isnan(x) || isnan(y)) { 371 /* cacos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */ 372 if (isinf(x)) 373 return (cpack(y + y, -INFINITY)); 374 /* cacos(NaN + I*+-Inf) = NaN + I*-+Inf */ 375 if (isinf(y)) 376 return (cpack(x + x, -y)); 377 /* cacos(0 + I*NaN) = PI/2 + I*NaN with inexact */ 378 if (x == 0) 379 return (cpack(pio2_hi + pio2_lo, y + y)); 380 /* 381 * All other cases involving NaN return NaN + I*NaN. 382 * C99 leaves it optional whether to raise invalid if one of 383 * the arguments is not NaN, so we opt not to raise it. 384 */ 385 return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); 386 } 387 388 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { 389 /* clog...() will raise inexact unless x or y is infinite. */ 390 w = clog_for_large_values(z); 391 rx = fabs(cimag(w)); 392 ry = creal(w) + m_ln2; 393 if (sy == 0) 394 ry = -ry; 395 return (cpack(rx, ry)); 396 } 397 398 /* Avoid spuriously raising inexact for z = 1. */ 399 if (x == 1 && y == 0) 400 return (cpack(0, -y)); 401 402 /* All remaining cases are inexact. */ 403 raise_inexact(); 404 405 if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) 406 return (cpack(pio2_hi - (x - pio2_lo), -y)); 407 408 do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x); 409 if (B_is_usable) { 410 if (sx == 0) 411 rx = acos(B); 412 else 413 rx = acos(-B); 414 } else { 415 if (sx == 0) 416 rx = atan2(sqrt_A2mx2, new_x); 417 else 418 rx = atan2(sqrt_A2mx2, -new_x); 419 } 420 if (sy == 0) 421 ry = -ry; 422 return (cpack(rx, ry)); 423 } 424 425 /* 426 * cacosh(z) = I*cacos(z) or -I*cacos(z) 427 * where the sign is chosen so Re(cacosh(z)) >= 0. 428 */ 429 double complex 430 cacosh(double complex z) 431 { 432 double complex w; 433 double rx, ry; 434 435 w = cacos(z); 436 rx = creal(w); 437 ry = cimag(w); 438 /* cacosh(NaN + I*NaN) = NaN + I*NaN */ 439 if (isnan(rx) && isnan(ry)) 440 return (cpack(ry, rx)); 441 /* cacosh(NaN + I*+-Inf) = +Inf + I*NaN */ 442 /* cacosh(+-Inf + I*NaN) = +Inf + I*NaN */ 443 if (isnan(rx)) 444 return (cpack(fabs(ry), rx)); 445 /* cacosh(0 + I*NaN) = NaN + I*NaN */ 446 if (isnan(ry)) 447 return (cpack(ry, ry)); 448 return (cpack(fabs(ry), copysign(rx, cimag(z)))); 449 } 450 451 /* 452 * Optimized version of clog() for |z| finite and larger than ~RECIP_EPSILON. 453 */ 454 static double complex 455 clog_for_large_values(double complex z) 456 { 457 double x, y; 458 double ax, ay, t; 459 460 x = creal(z); 461 y = cimag(z); 462 ax = fabs(x); 463 ay = fabs(y); 464 if (ax < ay) { 465 t = ax; 466 ax = ay; 467 ay = t; 468 } 469 470 /* 471 * Avoid overflow in hypot() when x and y are both very large. 472 * Divide x and y by E, and then add 1 to the logarithm. This depends 473 * on E being larger than sqrt(2). 474 * Dividing by E causes an insignificant loss of accuracy; however 475 * this method is still poor since it is uneccessarily slow. 476 */ 477 if (ax > DBL_MAX / 2) 478 return (cpack(log(hypot(x / m_e, y / m_e)) + 1, atan2(y, x))); 479 480 /* 481 * Avoid overflow when x or y is large. Avoid underflow when x or 482 * y is small. 483 */ 484 if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN) 485 return (cpack(log(hypot(x, y)), atan2(y, x))); 486 487 return (cpack(log(ax * ax + ay * ay) / 2, atan2(y, x))); 488 } 489 490 /* 491 * ================= 492 * | catanh, catan | 493 * ================= 494 */ 495 496 /* 497 * sum_squares(x,y) = x*x + y*y (or just x*x if y*y would underflow). 498 * Assumes x*x and y*y will not overflow. 499 * Assumes x and y are finite. 500 * Assumes y is non-negative. 501 * Assumes fabs(x) >= DBL_EPSILON. 502 */ 503 static inline double 504 sum_squares(double x, double y) 505 { 506 507 /* Avoid underflow when y is small. */ 508 if (y < SQRT_MIN) 509 return (x * x); 510 511 return (x * x + y * y); 512 } 513 514 /* 515 * real_part_reciprocal(x, y) = Re(1/(x+I*y)) = x/(x*x + y*y). 516 * Assumes x and y are not NaN, and one of x and y is larger than 517 * RECIP_EPSILON. We avoid unwarranted underflow. It is important to not use 518 * the code creal(1/z), because the imaginary part may produce an unwanted 519 * underflow. 520 * This is only called in a context where inexact is always raised before 521 * the call, so no effort is made to avoid or force inexact. 522 */ 523 static inline double 524 real_part_reciprocal(double x, double y) 525 { 526 double scale; 527 uint32_t hx, hy; 528 int32_t ix, iy; 529 530 /* 531 * This code is inspired by the C99 document n1124.pdf, Section G.5.1, 532 * example 2. 533 */ 534 GET_HIGH_WORD(hx, x); 535 ix = hx & 0x7ff00000; 536 GET_HIGH_WORD(hy, y); 537 iy = hy & 0x7ff00000; 538 #define BIAS (DBL_MAX_EXP - 1) 539 /* XXX more guard digits are useful iff there is extra precision. */ 540 #define CUTOFF (DBL_MANT_DIG / 2 + 1) /* just half or 1 guard digit */ 541 if (ix - iy >= CUTOFF << 20 || isinf(x)) 542 return (1 / x); /* +-Inf -> +-0 is special */ 543 if (iy - ix >= CUTOFF << 20) 544 return (x / y / y); /* should avoid double div, but hard */ 545 if (ix <= (BIAS + DBL_MAX_EXP / 2 - CUTOFF) << 20) 546 return (x / (x * x + y * y)); 547 scale = 1; 548 SET_HIGH_WORD(scale, 0x7ff00000 - ix); /* 2**(1-ilogb(x)) */ 549 x *= scale; 550 y *= scale; 551 return (x / (x * x + y * y) * scale); 552 } 553 554 /* 555 * catanh(z) = log((1+z)/(1-z)) / 2 556 * = log1p(4*x / |z-1|^2) / 4 557 * + I * atan2(2*y, (1-x)*(1+x)-y*y) / 2 558 * 559 * catanh(z) = z + O(z^3) as z -> 0 560 * 561 * catanh(z) = 1/z + sign(y)*I*PI/2 + O(1/z^3) as z -> infinity 562 * The above formula works for the real part as well, because 563 * Re(catanh(z)) = x/|z|^2 + O(x/z^4) 564 * as z -> infinity, uniformly in x 565 */ 566 double complex 567 catanh(double complex z) 568 { 569 double x, y, ax, ay, rx, ry; 570 571 x = creal(z); 572 y = cimag(z); 573 ax = fabs(x); 574 ay = fabs(y); 575 576 /* This helps handle many cases. */ 577 if (y == 0 && ax <= 1) 578 return (cpack(atanh(x), y)); 579 580 /* To ensure the same accuracy as atan(), and to filter out z = 0. */ 581 if (x == 0) 582 return (cpack(x, atan(y))); 583 584 if (isnan(x) || isnan(y)) { 585 /* catanh(+-Inf + I*NaN) = +-0 + I*NaN */ 586 if (isinf(x)) 587 return (cpack(copysign(0, x), y + y)); 588 /* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */ 589 if (isinf(y)) 590 return (cpack(copysign(0, x), 591 copysign(pio2_hi + pio2_lo, y))); 592 /* 593 * All other cases involving NaN return NaN + I*NaN. 594 * C99 leaves it optional whether to raise invalid if one of 595 * the arguments is not NaN, so we opt not to raise it. 596 */ 597 return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); 598 } 599 600 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) 601 return (cpack(real_part_reciprocal(x, y), 602 copysign(pio2_hi + pio2_lo, y))); 603 604 if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) { 605 /* 606 * z = 0 was filtered out above. All other cases must raise 607 * inexact, but this is the only only that needs to do it 608 * explicitly. 609 */ 610 raise_inexact(); 611 return (z); 612 } 613 614 if (ax == 1 && ay < DBL_EPSILON) 615 rx = (m_ln2 - log(ay)) / 2; 616 else 617 rx = log1p(4 * ax / sum_squares(ax - 1, ay)) / 4; 618 619 if (ax == 1) 620 ry = atan2(2, -ay) / 2; 621 else if (ay < DBL_EPSILON) 622 ry = atan2(2 * ay, (1 - ax) * (1 + ax)) / 2; 623 else 624 ry = atan2(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2; 625 626 return (cpack(copysign(rx, x), copysign(ry, y))); 627 } 628 629 /* 630 * catan(z) = reverse(catanh(reverse(z))) 631 * where reverse(x + I*y) = y + I*x = I*conj(z). 632 */ 633 double complex 634 catan(double complex z) 635 { 636 double complex w = catanh(cpack(cimag(z), creal(z))); 637 638 return (cpack(cimag(w), creal(w))); 639 } 640