src/complex.c

   1 #if defined(__GNUC__)
   2 #if !defined _GNU_SOURCE
   3 #define _GNU_SOURCE 1
   4 #endif /* _GNU_SOURCE */
   5 #endif /* __GNUC__ */
   6 #include "a/a.h"
   7 #if A_PREREQ_GNUC(3, 0)
   8 #pragma GCC diagnostic ignored "-Wfloat-conversion"
   9 #endif /* -Wfloat-conversion */
  10 #if A_PREREQ_GNUC(3, 0) || __has_warning("-Wfloat-equal")
  11 #pragma GCC diagnostic ignored "-Wfloat-equal"
  12 #endif /* -Wfloat-equal */
  13 #define LIBA_COMPLEX_C
  14 #include "a/complex.h"
  15 /* compiler built-in complex number type */
  16 #if A_PREREQ_MSVC(18, 0)
  17 #include <complex.h> /* 12.0 */
  18 #if A_FLOAT_TYPE + 0 == A_FLOAT_SINGLE
  19 #define A_COMPLEX _Fcomplex
  20 #elif A_FLOAT_TYPE + 0 == A_FLOAT_DOUBLE
  21 #define A_COMPLEX _Dcomplex
  22 #elif A_FLOAT_TYPE + 0 == A_FLOAT_EXTEND
  23 #define A_COMPLEX _Lcomplex
  24 #endif /* A_COMPLEX */
  25 #elif defined(__GNUC__) || defined(__clang__)
  26 #include <complex.h>
  27 #define A_COMPLEX _Complex A_FLOAT
  28 #else /* !A_FLOAT_TYPE */
  29 #if A_PREREQ_GNUC(2, 95) || __has_warning("-Waggregate-return")
  30 #pragma GCC diagnostic ignored "-Waggregate-return"
  31 #endif /* -Waggregate-return */
  32 #endif /* A_FLOAT_TYPE */
  33 #include "a/math.h"
  34 #include <stdlib.h>
  35
  36 #if !defined __STDC_VERSION__ || (defined(_MSC_VER) && (_MSC_VER < 1800))
  37 #if !defined isinf
  38 #define isinf(x) ((x) + (x) == (x) && (x))
  39 #endif /* isinf */
  40 #endif /* __STDC_VERSION__ */
  41
  42 void a_complex_polar(a_complex *ctx, a_float rho, a_float theta)
  43 {
  44     ctx->real = rho * a_float_cos(theta);
  45     ctx->imag = rho * a_float_sin(theta);
  46 }
  47
  48 unsigned int a_complex_parse(a_complex *ctx, char const *str)
  49 {
  50 #if A_FLOAT_TYPE + 0 == A_FLOAT_SINGLE
  51 #define strtonum(string, endptr) strtof(string, endptr)
  52 #elif A_FLOAT_TYPE + 0 == A_FLOAT_DOUBLE
  53 #define strtonum(string, endptr) strtod(string, endptr)
  54 #elif A_FLOAT_TYPE + 0 == A_FLOAT_EXTEND
  55 #define strtonum(string, endptr) strtold(string, endptr)
  56 #endif /* A_FLOAT_TYPE */
  57     union
  58     {
  59         char const *s;
  60         char *p;
  61     } u;
  62     u.s = str;
  63     if (!str) { return 0; }
  64     for (ctx->real = 0; *u.s; ++u.s)
  65     {
  66         if (*u.s == '+' || *u.s == '-' || ('0' <= *u.s && *u.s <= '9') || *u.s == '.')
  67         {
  68             ctx->real = strtonum(u.s, &u.p);
  69             break;
  70         }
  71     }
  72     for (ctx->imag = 0; *u.s; ++u.s)
  73     {
  74         if (*u.s == '+' || *u.s == '-' || ('0' <= *u.s && *u.s <= '9') || *u.s == '.')
  75         {
  76             ctx->imag = strtonum(u.s, &u.p);
  77             break;
  78         }
  79     }
  80     return (unsigned int)(u.s - str);
  81 }
  82
  83 a_bool a_complex_eq(a_complex x, a_complex y)
  84 {
  85     return x.real == y.real && x.imag == y.imag;
  86 }
  87
  88 a_bool a_complex_ne(a_complex x, a_complex y)
  89 {
  90     return x.real != y.real || x.imag != y.imag;
  91 }
  92
  93 a_float a_complex_logabs(a_complex z)
  94 {
  95     a_float r, u;
  96     a_float const xabs = a_float_abs(z.real);
  97     a_float const yabs = a_float_abs(z.imag);
  98     if (xabs >= yabs)
  99     {
 100         r = xabs;
 101         u = yabs / xabs;
 102     }
 103     else
 104     {
 105         r = yabs;
 106         u = xabs / yabs;
 107     }
 108     return a_float_log(r) + A_FLOAT_C(0.5) * a_float_log1p(u * u);
 109 }
 110
 111 a_float a_complex_abs2(a_complex z)
 112 {
 113     return z.real * z.real + z.imag * z.imag;
 114 }
 115
 116 a_float a_complex_abs(a_complex z)
 117 {
 118     return a_float_hypot(z.real, z.imag);
 119 }
 120
 121 a_float a_complex_arg(a_complex z)
 122 {
 123     if (z.real != 0 || z.imag != 0)
 124     {
 125         return a_float_atan2(z.imag, z.real);
 126     }
 127     return 0;
 128 }
 129
 130 void a_complex_proj_(a_complex *ctx)
 131 {
 132     if (isinf(ctx->real) || isinf(ctx->imag))
 133     {
 134         ctx->real = A_FLOAT_INF;
 135         ctx->imag = ctx->imag < 0 ? A_FLOAT_C(-0.0) : A_FLOAT_C(0.0);
 136     }
 137 }
 138
 139 void a_complex_mul_(a_complex *ctx, a_complex z)
 140 {
 141     a_float const real = ctx->real;
 142     ctx->real = real * z.real - ctx->imag * z.imag;
 143     ctx->imag = real * z.imag + ctx->imag * z.real;
 144 }
 145
 146 void a_complex_div_(a_complex *ctx, a_complex z)
 147 {
 148     a_float const inv = 1 / a_complex_abs(z);
 149     a_float const xr = ctx->real * inv;
 150     a_float const xi = ctx->imag * inv;
 151     a_float const yr = z.real * inv;
 152     a_float const yi = z.imag * inv;
 153     ctx->real = xr * yr + xi * yi;
 154     ctx->imag = xi * yr - xr * yi;
 155 }
 156
 157 void a_complex_inv_(a_complex *ctx)
 158 {
 159     a_float const inv = 1 / a_complex_abs(*ctx);
 160     ctx->real = +inv * ctx->real * inv;
 161     ctx->imag = -inv * ctx->imag * inv;
 162 }
 163
 164 #if A_PREREQ_GNUC(2, 95) || __has_warning("-Wimplicit-function-declaration")
 165 #pragma GCC diagnostic ignored "-Wimplicit-function-declaration"
 166 #endif /* gcc 2.95+ */
 167 #if __has_warning("-Wincompatible-library-redeclaration")
 168 #pragma clang diagnostic ignored "-Wincompatible-library-redeclaration"
 169 #endif /* clang 4.0+ */
 170 #if A_PREREQ_GNUC(10, 0)
 171 #pragma GCC diagnostic ignored "-Wbuiltin-declaration-mismatch"
 172 #endif /* gcc 10.0+ */
 173
 174 #if defined(A_HAVE_CSQRT) && (A_HAVE_CSQRT + 0 < 1)
 175 #undef A_HAVE_CSQRT
 176 #endif /* A_HAVE_CSQRT */
 177 #if defined(A_HAVE_CSQRT) && !defined A_COMPLEX
 178 a_complex A_FLOAT_F(csqrt)(a_complex);
 179 #endif /* A_HAVE_CSQRT */
 180 void a_complex_sqrt_(a_complex *ctx)
 181 {
 182 #if defined(A_HAVE_CSQRT) && defined(A_COMPLEX)
 183     A_COMPLEX *const x = (A_COMPLEX *)ctx;
 184     *x = A_FLOAT_F(csqrt)(*x);
 185 #elif defined(A_HAVE_CSQRT)
 186     *ctx = A_FLOAT_F(csqrt)(*ctx);
 187 #else /* !A_HAVE_CSQRT */
 188     if (ctx->real != 0 || ctx->imag != 0)
 189     {
 190         a_float const x = a_float_abs(ctx->real);
 191         a_float const y = a_float_abs(ctx->imag);
 192         a_float w;
 193         if (x >= y)
 194         {
 195             a_float u = y / x;
 196             w = A_FLOAT_SQRT1_2 * a_float_sqrt(x) * a_float_sqrt(a_float_sqrt(u * u + 1) + 1);
 197         }
 198         else
 199         {
 200             a_float u = x / y;
 201             w = A_FLOAT_SQRT1_2 * a_float_sqrt(y) * a_float_sqrt(a_float_sqrt(u * u + 1) + u);
 202         }
 203         if (ctx->real >= 0)
 204         {
 205             ctx->real = w;
 206             ctx->imag = ctx->imag / (2 * w);
 207         }
 208         else
 209         {
 210             ctx->real = ctx->imag / (2 * w);
 211             ctx->imag = ctx->imag < 0 ? -w : w;
 212         }
 213     }
 214 #endif /* A_HAVE_CSQRT */
 215 }
 216
 217 void a_complex_sqrt_real(a_complex *ctx, a_float x)
 218 {
 219     if (x >= 0)
 220     {
 221         ctx->real = a_float_sqrt(x);
 222         ctx->imag = 0;
 223     }
 224     else
 225     {
 226         ctx->real = 0;
 227         ctx->imag = a_float_sqrt(-x);
 228     }
 229 }
 230
 231 #if defined(A_HAVE_CPOW) && (A_HAVE_CPOW + 0 < 1)
 232 #undef A_HAVE_CPOW
 233 #endif /* A_HAVE_CPOW */
 234 #if defined(A_HAVE_CPOW) && !defined A_COMPLEX
 235 a_complex A_FLOAT_F(cpow)(a_complex, a_complex);
 236 #endif /* A_HAVE_CPOW */
 237 void a_complex_pow_(a_complex *ctx, a_complex a)
 238 {
 239 #if defined(A_HAVE_CPOW) && defined(A_COMPLEX)
 240     A_COMPLEX *const x = (A_COMPLEX *)ctx;
 241     A_COMPLEX *const y = (A_COMPLEX *)&a;
 242     *x = A_FLOAT_F(cpow)(*x, *y);
 243 #elif defined(A_HAVE_CPOW)
 244     *ctx = A_FLOAT_F(cpow)(*ctx, a);
 245 #else /* !A_HAVE_CPOW */
 246     if (ctx->real != 0 || ctx->imag != 0)
 247     {
 248         a_float const logr = a_complex_logabs(*ctx);
 249         a_float const theta = a_complex_arg(*ctx);
 250         a_float const rho = a_float_exp(logr * a.real - theta * a.imag);
 251         a_float const beta = theta * a.real + logr * a.imag;
 252         a_complex_polar(ctx, rho, beta);
 253     }
 254     else
 255     {
 256         ctx->real = (a.real == 0 && a.imag == 0) ? 1 : 0;
 257     }
 258 #endif /* A_HAVE_CPOW */
 259 }
 260
 261 void a_complex_pow_real_(a_complex *ctx, a_float a)
 262 {
 263     if (ctx->real != 0 || ctx->imag != 0)
 264     {
 265         a_float const logr = a_complex_logabs(*ctx);
 266         a_float const theta = a_complex_arg(*ctx);
 267         a_float const rho = a_float_exp(logr * a);
 268         a_float const beta = theta * a;
 269         a_complex_polar(ctx, rho, beta);
 270     }
 271     else
 272     {
 273         ctx->real = (a == 0) ? 1 : 0;
 274     }
 275 }
 276
 277 #if defined(A_HAVE_CEXP) && (A_HAVE_CEXP + 0 < 1)
 278 #undef A_HAVE_CEXP
 279 #endif /* A_HAVE_CEXP */
 280 #if defined(A_HAVE_CEXP) && !defined A_COMPLEX
 281 a_complex A_FLOAT_F(cexp)(a_complex);
 282 #endif /* A_HAVE_CEXP */
 283 void a_complex_exp_(a_complex *ctx)
 284 {
 285 #if defined(A_HAVE_CEXP) && defined(A_COMPLEX)
 286     A_COMPLEX *const x = (A_COMPLEX *)ctx;
 287     *x = A_FLOAT_F(cexp)(*x);
 288 #elif defined(A_HAVE_CEXP)
 289     *ctx = A_FLOAT_F(cexp)(*ctx);
 290 #else /* !A_HAVE_CEXP */
 291     a_complex_polar(ctx, a_float_exp(ctx->real), ctx->imag);
 292 #endif /* A_HAVE_CEXP */
 293 }
 294
 295 #if defined(A_HAVE_CLOG) && (A_HAVE_CLOG + 0 < 1)
 296 #undef A_HAVE_CLOG
 297 #endif /* A_HAVE_CLOG */
 298 #if defined(A_HAVE_CLOG) && !defined A_COMPLEX
 299 a_complex A_FLOAT_F(clog)(a_complex);
 300 #endif /* A_HAVE_CLOG */
 301 void a_complex_log_(a_complex *ctx)
 302 {
 303 #if defined(A_HAVE_CLOG) && defined(A_COMPLEX)
 304     A_COMPLEX *const x = (A_COMPLEX *)ctx;
 305     *x = A_FLOAT_F(clog)(*x);
 306 #elif defined(A_HAVE_CLOG)
 307     *ctx = A_FLOAT_F(clog)(*ctx);
 308 #else /* !A_HAVE_CLOG */
 309     ctx->real = a_complex_logabs(*ctx);
 310     ctx->imag = a_complex_arg(*ctx);
 311 #endif /* A_HAVE_CLOG */
 312 }
 313
 314 void a_complex_log2_(a_complex *ctx)
 315 {
 316     a_complex_log_(ctx);
 317     a_complex_mul_real_(ctx, A_FLOAT_LN1_2);
 318 }
 319
 320 void a_complex_log10_(a_complex *ctx)
 321 {
 322     a_complex_log_(ctx);
 323     a_complex_mul_real_(ctx, A_FLOAT_LN1_10);
 324 }
 325
 326 void a_complex_logb_(a_complex *ctx, a_complex b)
 327 {
 328     a_complex_log_(ctx);
 329     a_complex_log_(&b);
 330     a_complex_div_(ctx, b);
 331 }
 332
 333 #if defined(A_HAVE_CSIN) && (A_HAVE_CSIN + 0 < 1)
 334 #undef A_HAVE_CSIN
 335 #endif /* A_HAVE_CSIN */
 336 #if defined(A_HAVE_CSIN) && !defined A_COMPLEX
 337 a_complex A_FLOAT_F(csin)(a_complex);
 338 #endif /* A_HAVE_CSIN */
 339 void a_complex_sin_(a_complex *ctx)
 340 {
 341 #if defined(A_HAVE_CSIN) && defined(A_COMPLEX)
 342     A_COMPLEX *const x = (A_COMPLEX *)ctx;
 343     *x = A_FLOAT_F(csin)(*x);
 344 #elif defined(A_HAVE_CSIN)
 345     *ctx = A_FLOAT_F(csin)(*ctx);
 346 #else /* !A_HAVE_CSIN */
 347     if (ctx->imag != 0)
 348     {
 349         a_float const real = ctx->real;
 350         ctx->real = a_float_sin(real) * a_float_cosh(ctx->imag);
 351         ctx->imag = a_float_cos(real) * a_float_sinh(ctx->imag);
 352     }
 353     else
 354     {
 355         ctx->real = a_float_sin(ctx->real);
 356     }
 357 #endif /* A_HAVE_CSIN */
 358 }
 359
 360 #if defined(A_HAVE_CCOS) && (A_HAVE_CCOS + 0 < 1)
 361 #undef A_HAVE_CCOS
 362 #endif /* A_HAVE_CCOS */
 363 #if defined(A_HAVE_CCOS) && !defined A_COMPLEX
 364 a_complex A_FLOAT_F(ccos)(a_complex);
 365 #endif /* A_HAVE_CCOS */
 366 void a_complex_cos_(a_complex *ctx)
 367 {
 368 #if defined(A_HAVE_CCOS) && defined(A_COMPLEX)
 369     A_COMPLEX *const x = (A_COMPLEX *)ctx;
 370     *x = A_FLOAT_F(ccos)(*x);
 371 #elif defined(A_HAVE_CCOS)
 372     *ctx = A_FLOAT_F(ccos)(*ctx);
 373 #else /* !A_HAVE_CCOS */
 374     if (ctx->imag != 0)
 375     {
 376         a_float const real = ctx->real;
 377         ctx->real = a_float_cos(real) * a_float_cosh(+ctx->imag);
 378         ctx->imag = a_float_sin(real) * a_float_sinh(-ctx->imag);
 379     }
 380     else
 381     {
 382         ctx->real = a_float_cos(ctx->real);
 383     }
 384 #endif /* A_HAVE_CCOS */
 385 }
 386
 387 #if defined(A_HAVE_CTAN) && (A_HAVE_CTAN + 0 < 1)
 388 #undef A_HAVE_CTAN
 389 #endif /* A_HAVE_CTAN */
 390 #if defined(A_HAVE_CTAN) && !defined A_COMPLEX
 391 a_complex A_FLOAT_F(ctan)(a_complex);
 392 #endif /* A_HAVE_CTAN */
 393 void a_complex_tan_(a_complex *ctx)
 394 {
 395 #if defined(A_HAVE_CTAN) && defined(A_COMPLEX)
 396     A_COMPLEX *const x = (A_COMPLEX *)ctx;
 397     *x = A_FLOAT_F(ctan)(*x);
 398 #elif defined(A_HAVE_CTAN)
 399     *ctx = A_FLOAT_F(ctan)(*ctx);
 400 #else /* !A_HAVE_CTAN */
 401     a_float const cr = a_float_cos(ctx->real);
 402     a_float const si = a_float_sinh(ctx->imag);
 403     a_float const den = cr * cr + si * si;
 404     ctx->real = A_FLOAT_C(0.5) * a_float_sin(2 * ctx->real) / den;
 405     if (a_float_abs(ctx->imag) < 1)
 406     {
 407         ctx->imag = A_FLOAT_C(0.5) * a_float_sinh(2 * ctx->imag) / den;
 408     }
 409     else
 410     {
 411         a_float const d = a_float_pow(cr / si, 2) + 1;
 412         ctx->imag = 1 / (a_float_tanh(ctx->imag) * d);
 413     }
 414 #endif /* A_HAVE_CTAN */
 415 }
 416
 417 void a_complex_sec_(a_complex *ctx)
 418 {
 419     a_complex_cos_(ctx);
 420     a_complex_inv_(ctx);
 421 }
 422
 423 void a_complex_csc_(a_complex *ctx)
 424 {
 425     a_complex_sin_(ctx);
 426     a_complex_inv_(ctx);
 427 }
 428
 429 void a_complex_cot_(a_complex *ctx)
 430 {
 431     a_complex_tan_(ctx);
 432     a_complex_inv_(ctx);
 433 }
 434
 435 #if defined(A_HAVE_CASIN) && (A_HAVE_CASIN + 0 < 1)
 436 #undef A_HAVE_CASIN
 437 #endif /* A_HAVE_CASIN */
 438 #if defined(A_HAVE_CASIN) && !defined A_COMPLEX
 439 a_complex A_FLOAT_F(casin)(a_complex);
 440 #endif /* A_HAVE_CASIN */
 441 void a_complex_asin_(a_complex *ctx)
 442 {
 443 #if defined(A_HAVE_CASIN) && defined(A_COMPLEX)
 444     A_COMPLEX *const x = (A_COMPLEX *)ctx;
 445     *x = A_FLOAT_F(casin)(*x);
 446 #elif defined(A_HAVE_CASIN)
 447     *ctx = A_FLOAT_F(casin)(*ctx);
 448 #else /* !A_HAVE_CASIN */
 449     if (ctx->imag != 0)
 450     {
 451         a_float const a_crossover = A_FLOAT_C(1.5);
 452         a_float const b_crossover = A_FLOAT_C(0.6417);
 453         a_float const x = a_float_abs(ctx->real);
 454         a_float const y = a_float_abs(ctx->imag);
 455         a_float const r = a_float_hypot(x + 1, y);
 456         a_float const s = a_float_hypot(x - 1, y);
 457         a_float const a = A_FLOAT_C(0.5) * (r + s);
 458         a_float const b = x / a;
 459         a_float const y2 = y * y;
 460         a_float const real = ctx->real;
 461         a_float const imag = ctx->imag;
 462         if (b <= b_crossover)
 463         {
 464             ctx->real = a_float_asin(b);
 465         }
 466         else if (x <= 1)
 467         {
 468             a_float const den = A_FLOAT_C(0.5) * (a + x) * (y2 / (r + x + 1) + (s + 1 - x));
 469             ctx->real = a_float_atan(x / a_float_sqrt(den));
 470         }
 471         else
 472         {
 473             a_float const apx = a + x;
 474             a_float const den = A_FLOAT_C(0.5) * (apx / (r + x + 1) + apx / (s + x - 1));
 475             ctx->real = a_float_atan(x / (a_float_sqrt(den) * y));
 476         }
 477         if (a <= a_crossover)
 478         {
 479             a_float am1;
 480             if (x < 1)
 481             {
 482                 am1 = A_FLOAT_C(0.5) * (y2 / (r + x + 1) + y2 / (s + 1 - x));
 483             }
 484             else
 485             {
 486                 am1 = A_FLOAT_C(0.5) * (y2 / (r + x + 1) + (s + x - 1));
 487             }
 488             ctx->imag = a_float_log1p(am1 + a_float_sqrt((a + 1) * am1));
 489         }
 490         else
 491         {
 492             ctx->imag = a_float_log(a + a_float_sqrt(a * a - 1));
 493         }
 494         if (real < 0) { ctx->real = -ctx->real; }
 495         if (imag < 0) { ctx->imag = -ctx->imag; }
 496     }
 497     else
 498     {
 499         a_complex_asin_real(ctx, ctx->real);
 500     }
 501 #endif /* A_HAVE_CASIN */
 502 }
 503
 504 void a_complex_asin_real(a_complex *ctx, a_float x)
 505 {
 506     if (a_float_abs(x) <= 1)
 507     {
 508         ctx->real = a_float_asin(x);
 509         ctx->imag = 0;
 510     }
 511     else
 512     {
 513         if (x > 0)
 514         {
 515             ctx->real = +A_FLOAT_PI_2;
 516             ctx->imag = -a_float_acosh(+x);
 517         }
 518         else
 519         {
 520             ctx->real = -A_FLOAT_PI_2;
 521             ctx->imag = +a_float_acosh(-x);
 522         }
 523     }
 524 }
 525
 526 #if defined(A_HAVE_CACOS) && (A_HAVE_CACOS + 0 < 1)
 527 #undef A_HAVE_CACOS
 528 #endif /* A_HAVE_CACOS */
 529 #if defined(A_HAVE_CACOS) && !defined A_COMPLEX
 530 a_complex A_FLOAT_F(cacos)(a_complex);
 531 #endif /* A_HAVE_CACOS */
 532 void a_complex_acos_(a_complex *ctx)
 533 {
 534 #if defined(A_HAVE_CACOS) && defined(A_COMPLEX)
 535     A_COMPLEX *const x = (A_COMPLEX *)ctx;
 536     *x = A_FLOAT_F(cacos)(*x);
 537 #elif defined(A_HAVE_CACOS)
 538     *ctx = A_FLOAT_F(cacos)(*ctx);
 539 #else /* !A_HAVE_CACOS */
 540     if (ctx->imag != 0)
 541     {
 542         a_float const a_crossover = A_FLOAT_C(1.5);
 543         a_float const b_crossover = A_FLOAT_C(0.6417);
 544         a_float const x = a_float_abs(ctx->real);
 545         a_float const y = a_float_abs(ctx->imag);
 546         a_float const r = a_float_hypot(x + 1, y);
 547         a_float const s = a_float_hypot(x - 1, y);
 548         a_float const a = A_FLOAT_C(0.5) * (r + s);
 549         a_float const b = x / a;
 550         a_float const y2 = y * y;
 551         a_float const real = ctx->real;
 552         a_float const imag = ctx->imag;
 553         if (b <= b_crossover)
 554         {
 555             ctx->real = a_float_acos(b);
 556         }
 557         else if (x <= 1)
 558         {
 559             a_float const den = A_FLOAT_C(0.5) * (a + x) * (y2 / (r + x + 1) + (s + 1 - x));
 560             ctx->real = a_float_atan(a_float_sqrt(den) / x);
 561         }
 562         else
 563         {
 564             a_float const apx = a + x;
 565             a_float const den = A_FLOAT_C(0.5) * (apx / (r + x + 1) + apx / (s + x - 1));
 566             ctx->real = a_float_atan(a_float_sqrt(den) * y / x);
 567         }
 568         if (a <= a_crossover)
 569         {
 570             a_float am1;
 571             if (x < 1)
 572             {
 573                 am1 = A_FLOAT_C(0.5) * (y2 / (r + x + 1) + y2 / (s + 1 - x));
 574             }
 575             else
 576             {
 577                 am1 = A_FLOAT_C(0.5) * (y2 / (r + x + 1) + (s + x - 1));
 578             }
 579             ctx->imag = a_float_log1p(am1 + a_float_sqrt((a + 1) * am1));
 580         }
 581         else
 582         {
 583             ctx->imag = a_float_log(a + a_float_sqrt(a * a - 1));
 584         }
 585         if (real < 0) { ctx->real = A_FLOAT_PI - ctx->real; }
 586         if (imag >= 0) { ctx->imag = -ctx->imag; }
 587     }
 588     else
 589     {
 590         a_complex_acos_real(ctx, ctx->real);
 591     }
 592 #endif /* A_HAVE_CACOS */
 593 }
 594
 595 void a_complex_acos_real(a_complex *ctx, a_float x)
 596 {
 597     if (a_float_abs(x) <= 1)
 598     {
 599         ctx->real = a_float_acos(x);
 600         ctx->imag = 0;
 601     }
 602     else
 603     {
 604         if (x > 0)
 605         {
 606             ctx->real = 0;
 607             ctx->imag = +a_float_acosh(+x);
 608         }
 609         else
 610         {
 611             ctx->real = A_FLOAT_PI;
 612             ctx->imag = -a_float_acosh(-x);
 613         }
 614     }
 615 }
 616
 617 #if defined(A_HAVE_CATAN) && (A_HAVE_CATAN + 0 < 1)
 618 #undef A_HAVE_CATAN
 619 #endif /* A_HAVE_CATAN */
 620 #if defined(A_HAVE_CATAN) && !defined A_COMPLEX
 621 a_complex A_FLOAT_F(catan)(a_complex);
 622 #endif /* A_HAVE_CATAN */
 623 void a_complex_atan_(a_complex *ctx)
 624 {
 625 #if defined(A_HAVE_CATAN) && defined(A_COMPLEX)
 626     A_COMPLEX *const x = (A_COMPLEX *)ctx;
 627     *x = A_FLOAT_F(catan)(*x);
 628 #elif defined(A_HAVE_CATAN)
 629     *ctx = A_FLOAT_F(catan)(*ctx);
 630 #else /* !A_HAVE_CATAN */
 631     if (ctx->imag != 0)
 632     {
 633         a_float const r = a_float_hypot(ctx->real, ctx->imag);
 634         a_float const u = 2 * ctx->imag / (r * r + 1);
 635         a_float const imag = ctx->imag;
 636         if (a_float_abs(u) < A_FLOAT_C(0.1))
 637         {
 638             ctx->imag = A_FLOAT_C(0.25) * (a_float_log1p(u) - a_float_log1p(-u));
 639         }
 640         else
 641         {
 642             a_float const a = a_float_hypot(ctx->real, ctx->imag + 1);
 643             a_float const b = a_float_hypot(ctx->real, ctx->imag - 1);
 644             ctx->imag = A_FLOAT_C(0.5) * a_float_log(a / b);
 645         }
 646         if (ctx->real != 0)
 647         {
 648             ctx->real = A_FLOAT_C(0.5) * a_float_atan2(2 * ctx->real, (1 + r) * (1 - r));
 649         }
 650         else
 651         {
 652             if (imag > 1)
 653             {
 654                 ctx->real = +A_FLOAT_PI;
 655             }
 656             else if (imag < -1)
 657             {
 658                 ctx->real = -A_FLOAT_PI;
 659             }
 660             else
 661             {
 662                 ctx->real = 0;
 663             }
 664         }
 665     }
 666     else
 667     {
 668         ctx->real = a_float_atan(ctx->real);
 669     }
 670 #endif /* A_HAVE_CATAN */
 671 }
 672
 673 void a_complex_asec_(a_complex *ctx)
 674 {
 675     a_complex_inv_(ctx);
 676     a_complex_acos_(ctx);
 677 }
 678
 679 void a_complex_asec_real(a_complex *ctx, a_float x)
 680 {
 681     if (x <= -1 || x >= 1)
 682     {
 683         ctx->real = a_float_acos(1 / x);
 684         ctx->imag = 0;
 685     }
 686     else
 687     {
 688         if (x >= 0)
 689         {
 690             ctx->real = 0;
 691             ctx->imag = +a_float_acosh(+1 / x);
 692         }
 693         else
 694         {
 695             ctx->real = A_FLOAT_PI;
 696             ctx->imag = -a_float_acosh(-1 / x);
 697         }
 698     }
 699 }
 700
 701 void a_complex_acsc_(a_complex *ctx)
 702 {
 703     a_complex_inv_(ctx);
 704     a_complex_asin_(ctx);
 705 }
 706
 707 void a_complex_acsc_real(a_complex *ctx, a_float x)
 708 {
 709     if (x <= -1 || x >= 1)
 710     {
 711         ctx->real = a_float_asin(1 / x);
 712         ctx->imag = 0;
 713     }
 714     else
 715     {
 716         if (x >= 0)
 717         {
 718             ctx->real = +A_FLOAT_PI_2;
 719             ctx->imag = -a_float_acosh(+1 / x);
 720         }
 721         else
 722         {
 723             ctx->real = -A_FLOAT_PI_2;
 724             ctx->imag = +a_float_acosh(-1 / x);
 725         }
 726     }
 727 }
 728
 729 void a_complex_acot_(a_complex *ctx)
 730 {
 731     if (ctx->real != 0 || ctx->imag != 0)
 732     {
 733         a_complex_inv_(ctx);
 734         a_complex_atan_(ctx);
 735     }
 736     else
 737     {
 738         ctx->real = A_FLOAT_PI_2;
 739     }
 740 }
 741
 742 #if defined(A_HAVE_CSINH) && (A_HAVE_CSINH + 0 < 1)
 743 #undef A_HAVE_CSINH
 744 #endif /* A_HAVE_CSINH */
 745 #if defined(A_HAVE_CSINH) && !defined A_COMPLEX
 746 a_complex A_FLOAT_F(csinh)(a_complex);
 747 #endif /* A_HAVE_CSINH */
 748 void a_complex_sinh_(a_complex *ctx)
 749 {
 750 #if defined(A_HAVE_CSINH) && defined(A_COMPLEX)
 751     A_COMPLEX *const x = (A_COMPLEX *)ctx;
 752     *x = A_FLOAT_F(csinh)(*x);
 753 #elif defined(A_HAVE_CSINH)
 754     *ctx = A_FLOAT_F(csinh)(*ctx);
 755 #else /* !A_HAVE_CSINH */
 756     a_float const real = ctx->real;
 757     ctx->real = a_float_sinh(real) * a_float_cos(ctx->imag);
 758     ctx->imag = a_float_cosh(real) * a_float_sin(ctx->imag);
 759 #endif /* A_HAVE_CSINH */
 760 }
 761
 762 #if defined(A_HAVE_CCOSH) && (A_HAVE_CCOSH + 0 < 1)
 763 #undef A_HAVE_CCOSH
 764 #endif /* A_HAVE_CCOSH */
 765 #if defined(A_HAVE_CCOSH) && !defined A_COMPLEX
 766 a_complex A_FLOAT_F(ccosh)(a_complex);
 767 #endif /* A_HAVE_CCOSH */
 768 void a_complex_cosh_(a_complex *ctx)
 769 {
 770 #if defined(A_HAVE_CCOSH) && defined(A_COMPLEX)
 771     A_COMPLEX *const x = (A_COMPLEX *)ctx;
 772     *x = A_FLOAT_F(ccosh)(*x);
 773 #elif defined(A_HAVE_CCOSH)
 774     *ctx = A_FLOAT_F(ccosh)(*ctx);
 775 #else /* !A_HAVE_CCOSH */
 776     a_float const real = ctx->real;
 777     ctx->real = a_float_cosh(real) * a_float_cos(ctx->imag);
 778     ctx->imag = a_float_sinh(real) * a_float_sin(ctx->imag);
 779 #endif /* A_HAVE_CCOSH */
 780 }
 781
 782 #if defined(A_HAVE_CTANH) && (A_HAVE_CTANH + 0 < 1)
 783 #undef A_HAVE_CTANH
 784 #endif /* A_HAVE_CTANH */
 785 #if defined(A_HAVE_CTANH) && !defined A_COMPLEX
 786 a_complex A_FLOAT_F(ctanh)(a_complex);
 787 #endif /* A_HAVE_CTANH */
 788 void a_complex_tanh_(a_complex *ctx)
 789 {
 790 #if defined(A_HAVE_CTANH) && defined(A_COMPLEX)
 791     A_COMPLEX *const x = (A_COMPLEX *)ctx;
 792     *x = A_FLOAT_F(ctanh)(*x);
 793 #elif defined(A_HAVE_CTANH)
 794     *ctx = A_FLOAT_F(ctanh)(*ctx);
 795 #else /* !A_HAVE_CTANH */
 796     a_float const ci = a_float_cos(ctx->imag);
 797     a_float const sr = a_float_sinh(ctx->real);
 798     a_float const den = ci * ci + sr * sr;
 799     ctx->imag = A_FLOAT_C(0.5) * a_float_sin(2 * ctx->imag) / den;
 800     if (a_float_abs(ctx->real) < 1)
 801     {
 802         ctx->real = a_float_sinh(ctx->real) * a_float_cosh(ctx->real) / den;
 803     }
 804     else
 805     {
 806         a_float const d = a_float_pow(ci / sr, 2) + 1;
 807         ctx->real = 1 / (a_float_tanh(ctx->real) * d);
 808     }
 809 #endif /* A_HAVE_CTANH */
 810 }
 811
 812 void a_complex_sech_(a_complex *ctx)
 813 {
 814     a_complex_cosh_(ctx);
 815     a_complex_inv_(ctx);
 816 }
 817
 818 void a_complex_csch_(a_complex *ctx)
 819 {
 820     a_complex_sinh_(ctx);
 821     a_complex_inv_(ctx);
 822 }
 823
 824 void a_complex_coth_(a_complex *ctx)
 825 {
 826     a_complex_tanh_(ctx);
 827     a_complex_inv_(ctx);
 828 }
 829
 830 #if defined(A_HAVE_CASINH) && (A_HAVE_CASINH + 0 < 1)
 831 #undef A_HAVE_CASINH
 832 #endif /* A_HAVE_CASINH */
 833 #if defined(A_HAVE_CASINH) && !defined A_COMPLEX
 834 a_complex A_FLOAT_F(casinh)(a_complex);
 835 #endif /* A_HAVE_CASINH */
 836 void a_complex_asinh_(a_complex *ctx)
 837 {
 838 #if defined(A_HAVE_CASINH) && defined(A_COMPLEX)
 839     A_COMPLEX *const x = (A_COMPLEX *)ctx;
 840     *x = A_FLOAT_F(casinh)(*x);
 841 #elif defined(A_HAVE_CASINH)
 842     *ctx = A_FLOAT_F(casinh)(*ctx);
 843 #else /* !A_HAVE_CASINH */
 844     a_complex_mul_imag_(ctx, +1);
 845     a_complex_asin_(ctx);
 846     a_complex_mul_imag_(ctx, -1);
 847 #endif /* A_HAVE_CASINH */
 848 }
 849
 850 #if defined(A_HAVE_CACOSH) && (A_HAVE_CACOSH + 0 < 1)
 851 #undef A_HAVE_CACOSH
 852 #endif /* A_HAVE_CACOSH */
 853 #if defined(A_HAVE_CACOSH) && !defined A_COMPLEX
 854 a_complex A_FLOAT_F(cacosh)(a_complex);
 855 #endif /* A_HAVE_CACOSH */
 856 void a_complex_acosh_(a_complex *ctx)
 857 {
 858 #if defined(A_HAVE_CACOSH) && defined(A_COMPLEX)
 859     A_COMPLEX *const x = (A_COMPLEX *)ctx;
 860     *x = A_FLOAT_F(cacosh)(*x);
 861 #elif defined(A_HAVE_CACOSH)
 862     *ctx = A_FLOAT_F(cacosh)(*ctx);
 863 #else /* !A_HAVE_CACOSH */
 864     a_complex_acsc_(ctx);
 865     a_complex_mul_imag_(ctx, ctx->imag > 0 ? -1 : +1);
 866 #endif /* A_HAVE_CACOSH */
 867 }
 868
 869 void a_complex_acosh_real(a_complex *ctx, a_float x)
 870 {
 871     if (x >= 1)
 872     {
 873         ctx->real = a_float_acosh(+x);
 874         ctx->imag = 0;
 875     }
 876     else if (x >= -1)
 877     {
 878         ctx->real = 0;
 879         ctx->imag = a_float_acos(x);
 880     }
 881     else
 882     {
 883         ctx->real = a_float_acosh(-x);
 884         ctx->imag = A_FLOAT_PI;
 885     }
 886 }
 887
 888 #undef A_HAVE_CATANH
 889 #if defined(A_HAVE_CATANH) && (A_HAVE_CATANH + 0 < 1)
 890 #undef A_HAVE_CATANH
 891 #endif /* A_HAVE_CATANH */
 892 #if defined(A_HAVE_CATANH) && !defined A_COMPLEX
 893 a_complex A_FLOAT_F(catanh)(a_complex);
 894 #endif /* A_HAVE_CATANH */
 895 void a_complex_atanh_(a_complex *ctx)
 896 {
 897 #if defined(A_HAVE_CATANH) && defined(A_COMPLEX)
 898     A_COMPLEX *const x = (A_COMPLEX *)ctx;
 899     *x = A_FLOAT_F(catanh)(*x);
 900 #elif defined(A_HAVE_CATANH)
 901     *ctx = A_FLOAT_F(catanh)(*ctx);
 902 #else /* !A_HAVE_CATANH */
 903     if (ctx->imag != 0)
 904     {
 905         a_complex_mul_imag_(ctx, +1);
 906         a_complex_atan_(ctx);
 907         a_complex_mul_imag_(ctx, -1);
 908     }
 909     else
 910     {
 911         a_complex_atanh_real(ctx, ctx->real);
 912     }
 913 #endif /* A_HAVE_CATANH */
 914 }
 915
 916 void a_complex_atanh_real(a_complex *ctx, a_float x)
 917 {
 918     if (x > -1 && x < 1)
 919     {
 920         ctx->real = a_float_atanh(x);
 921         ctx->imag = 0;
 922     }
 923     else
 924     {
 925         ctx->real = a_float_atanh(1 / x);
 926         ctx->imag = (x < 0) ? +A_FLOAT_PI_2 : -A_FLOAT_PI_2;
 927     }
 928 }
 929
 930 void a_complex_asech_(a_complex *ctx)
 931 {
 932     a_complex_inv_(ctx);
 933     a_complex_acosh_(ctx);
 934 }
 935
 936 void a_complex_acsch_(a_complex *ctx)
 937 {
 938     a_complex_inv_(ctx);
 939     a_complex_asinh_(ctx);
 940 }
 941
 942 void a_complex_acoth_(a_complex *ctx)
 943 {
 944     a_complex_inv_(ctx);
 945     a_complex_atanh_(ctx);
 946 }