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