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