src/complex.c

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