lib/msun/src/catrig.c

   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(nan_mix(x, y), nan_mix(x, y)));
 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(nan_mix(x, y), nan_mix(x, y)));
 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(nan_mix(x, y), nan_mix(x, y)));
 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 case 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