src/complex.c

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