libgfortran/intrinsics/c99_functions.c

   1 /* Implementation of various C99 functions
   2    Copyright (C) 2004-2014 Free Software Foundation, Inc.
   3
   4 This file is part of the GNU Fortran 95 runtime library (libgfortran).
   5
   6 Libgfortran is free software; you can redistribute it and/or
   7 modify it under the terms of the GNU General Public
   8 License as published by the Free Software Foundation; either
   9 version 3 of the License, or (at your option) any later version.
  10
  11 Libgfortran is distributed in the hope that it will be useful,
  12 but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 GNU General Public License for more details.
  15
  16 Under Section 7 of GPL version 3, you are granted additional
  17 permissions described in the GCC Runtime Library Exception, version
  18 3.1, as published by the Free Software Foundation.
  19
  20 You should have received a copy of the GNU General Public License and
  21 a copy of the GCC Runtime Library Exception along with this program;
  22 see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  23 <http://www.gnu.org/licenses/>.  */
  24
  25 #include "config.h"
  26
  27 #define C99_PROTOS_H WE_DONT_WANT_PROTOS_NOW
  28 #include "libgfortran.h"
  29
  30 /* On a C99 system "I" (with I*I = -1) should be defined in complex.h;
  31    if not, we define a fallback version here.  */
  32 #ifndef I
  33 # if defined(_Imaginary_I)
  34 #   define I _Imaginary_I
  35 # elif defined(_Complex_I)
  36 #   define I _Complex_I
  37 # else
  38 #   define I (1.0fi)
  39 # endif
  40 #endif
  41
  42 /* Macros to get real and imaginary parts of a complex, and set
  43    a complex value.  */
  44 #define REALPART(z) (__real__(z))
  45 #define IMAGPART(z) (__imag__(z))
  46 #define COMPLEX_ASSIGN(z_, r_, i_) {__real__(z_) = (r_); __imag__(z_) = (i_);}
  47
  48
  49 /* Prototypes are included to silence -Wstrict-prototypes
  50    -Wmissing-prototypes.  */
  51
  52
  53 /* Wrappers for systems without the various C99 single precision Bessel
  54    functions.  */
  55
  56 #if defined(HAVE_J0) && ! defined(HAVE_J0F)
  57 #define HAVE_J0F 1
  58 float j0f (float);
  59
  60 float
  61 j0f (float x)
  62 {
  63   return (float) j0 ((double) x);
  64 }
  65 #endif
  66
  67 #if defined(HAVE_J1) && !defined(HAVE_J1F)
  68 #define HAVE_J1F 1
  69 float j1f (float);
  70
  71 float j1f (float x)
  72 {
  73   return (float) j1 ((double) x);
  74 }
  75 #endif
  76
  77 #if defined(HAVE_JN) && !defined(HAVE_JNF)
  78 #define HAVE_JNF 1
  79 float jnf (int, float);
  80
  81 float
  82 jnf (int n, float x)
  83 {
  84   return (float) jn (n, (double) x);
  85 }
  86 #endif
  87
  88 #if defined(HAVE_Y0) && !defined(HAVE_Y0F)
  89 #define HAVE_Y0F 1
  90 float y0f (float);
  91
  92 float
  93 y0f (float x)
  94 {
  95   return (float) y0 ((double) x);
  96 }
  97 #endif
  98
  99 #if defined(HAVE_Y1) && !defined(HAVE_Y1F)
 100 #define HAVE_Y1F 1
 101 float y1f (float);
 102
 103 float
 104 y1f (float x)
 105 {
 106   return (float) y1 ((double) x);
 107 }
 108 #endif
 109
 110 #if defined(HAVE_YN) && !defined(HAVE_YNF)
 111 #define HAVE_YNF 1
 112 float ynf (int, float);
 113
 114 float
 115 ynf (int n, float x)
 116 {
 117   return (float) yn (n, (double) x);
 118 }
 119 #endif
 120
 121
 122 /* Wrappers for systems without the C99 erff() and erfcf() functions.  */
 123
 124 #if defined(HAVE_ERF) && !defined(HAVE_ERFF)
 125 #define HAVE_ERFF 1
 126 float erff (float);
 127
 128 float
 129 erff (float x)
 130 {
 131   return (float) erf ((double) x);
 132 }
 133 #endif
 134
 135 #if defined(HAVE_ERFC) && !defined(HAVE_ERFCF)
 136 #define HAVE_ERFCF 1
 137 float erfcf (float);
 138
 139 float
 140 erfcf (float x)
 141 {
 142   return (float) erfc ((double) x);
 143 }
 144 #endif
 145
 146
 147 #ifndef HAVE_ACOSF
 148 #define HAVE_ACOSF 1
 149 float acosf (float x);
 150
 151 float
 152 acosf (float x)
 153 {
 154   return (float) acos (x);
 155 }
 156 #endif
 157
 158 #if HAVE_ACOSH && !HAVE_ACOSHF
 159 float acoshf (float x);
 160
 161 float
 162 acoshf (float x)
 163 {
 164   return (float) acosh ((double) x);
 165 }
 166 #endif
 167
 168 #ifndef HAVE_ASINF
 169 #define HAVE_ASINF 1
 170 float asinf (float x);
 171
 172 float
 173 asinf (float x)
 174 {
 175   return (float) asin (x);
 176 }
 177 #endif
 178
 179 #if HAVE_ASINH && !HAVE_ASINHF
 180 float asinhf (float x);
 181
 182 float
 183 asinhf (float x)
 184 {
 185   return (float) asinh ((double) x);
 186 }
 187 #endif
 188
 189 #ifndef HAVE_ATAN2F
 190 #define HAVE_ATAN2F 1
 191 float atan2f (float y, float x);
 192
 193 float
 194 atan2f (float y, float x)
 195 {
 196   return (float) atan2 (y, x);
 197 }
 198 #endif
 199
 200 #ifndef HAVE_ATANF
 201 #define HAVE_ATANF 1
 202 float atanf (float x);
 203
 204 float
 205 atanf (float x)
 206 {
 207   return (float) atan (x);
 208 }
 209 #endif
 210
 211 #if HAVE_ATANH && !HAVE_ATANHF
 212 float atanhf (float x);
 213
 214 float
 215 atanhf (float x)
 216 {
 217   return (float) atanh ((double) x);
 218 }
 219 #endif
 220
 221 #ifndef HAVE_CEILF
 222 #define HAVE_CEILF 1
 223 float ceilf (float x);
 224
 225 float
 226 ceilf (float x)
 227 {
 228   return (float) ceil (x);
 229 }
 230 #endif
 231
 232 #ifndef HAVE_COPYSIGNF
 233 #define HAVE_COPYSIGNF 1
 234 float copysignf (float x, float y);
 235
 236 float
 237 copysignf (float x, float y)
 238 {
 239   return (float) copysign (x, y);
 240 }
 241 #endif
 242
 243 #ifndef HAVE_COSF
 244 #define HAVE_COSF 1
 245 float cosf (float x);
 246
 247 float
 248 cosf (float x)
 249 {
 250   return (float) cos (x);
 251 }
 252 #endif
 253
 254 #ifndef HAVE_COSHF
 255 #define HAVE_COSHF 1
 256 float coshf (float x);
 257
 258 float
 259 coshf (float x)
 260 {
 261   return (float) cosh (x);
 262 }
 263 #endif
 264
 265 #ifndef HAVE_EXPF
 266 #define HAVE_EXPF 1
 267 float expf (float x);
 268
 269 float
 270 expf (float x)
 271 {
 272   return (float) exp (x);
 273 }
 274 #endif
 275
 276 #ifndef HAVE_FABSF
 277 #define HAVE_FABSF 1
 278 float fabsf (float x);
 279
 280 float
 281 fabsf (float x)
 282 {
 283   return (float) fabs (x);
 284 }
 285 #endif
 286
 287 #ifndef HAVE_FLOORF
 288 #define HAVE_FLOORF 1
 289 float floorf (float x);
 290
 291 float
 292 floorf (float x)
 293 {
 294   return (float) floor (x);
 295 }
 296 #endif
 297
 298 #ifndef HAVE_FMODF
 299 #define HAVE_FMODF 1
 300 float fmodf (float x, float y);
 301
 302 float
 303 fmodf (float x, float y)
 304 {
 305   return (float) fmod (x, y);
 306 }
 307 #endif
 308
 309 #ifndef HAVE_FREXPF
 310 #define HAVE_FREXPF 1
 311 float frexpf (float x, int *exp);
 312
 313 float
 314 frexpf (float x, int *exp)
 315 {
 316   return (float) frexp (x, exp);
 317 }
 318 #endif
 319
 320 #ifndef HAVE_HYPOTF
 321 #define HAVE_HYPOTF 1
 322 float hypotf (float x, float y);
 323
 324 float
 325 hypotf (float x, float y)
 326 {
 327   return (float) hypot (x, y);
 328 }
 329 #endif
 330
 331 #ifndef HAVE_LOGF
 332 #define HAVE_LOGF 1
 333 float logf (float x);
 334
 335 float
 336 logf (float x)
 337 {
 338   return (float) log (x);
 339 }
 340 #endif
 341
 342 #ifndef HAVE_LOG10F
 343 #define HAVE_LOG10F 1
 344 float log10f (float x);
 345
 346 float
 347 log10f (float x)
 348 {
 349   return (float) log10 (x);
 350 }
 351 #endif
 352
 353 #ifndef HAVE_SCALBN
 354 #define HAVE_SCALBN 1
 355 double scalbn (double x, int y);
 356
 357 double
 358 scalbn (double x, int y)
 359 {
 360 #if (FLT_RADIX == 2) && defined(HAVE_LDEXP)
 361   return ldexp (x, y);
 362 #else
 363   return x * pow (FLT_RADIX, y);
 364 #endif
 365 }
 366 #endif
 367
 368 #ifndef HAVE_SCALBNF
 369 #define HAVE_SCALBNF 1
 370 float scalbnf (float x, int y);
 371
 372 float
 373 scalbnf (float x, int y)
 374 {
 375   return (float) scalbn (x, y);
 376 }
 377 #endif
 378
 379 #ifndef HAVE_SINF
 380 #define HAVE_SINF 1
 381 float sinf (float x);
 382
 383 float
 384 sinf (float x)
 385 {
 386   return (float) sin (x);
 387 }
 388 #endif
 389
 390 #ifndef HAVE_SINHF
 391 #define HAVE_SINHF 1
 392 float sinhf (float x);
 393
 394 float
 395 sinhf (float x)
 396 {
 397   return (float) sinh (x);
 398 }
 399 #endif
 400
 401 #ifndef HAVE_SQRTF
 402 #define HAVE_SQRTF 1
 403 float sqrtf (float x);
 404
 405 float
 406 sqrtf (float x)
 407 {
 408   return (float) sqrt (x);
 409 }
 410 #endif
 411
 412 #ifndef HAVE_TANF
 413 #define HAVE_TANF 1
 414 float tanf (float x);
 415
 416 float
 417 tanf (float x)
 418 {
 419   return (float) tan (x);
 420 }
 421 #endif
 422
 423 #ifndef HAVE_TANHF
 424 #define HAVE_TANHF 1
 425 float tanhf (float x);
 426
 427 float
 428 tanhf (float x)
 429 {
 430   return (float) tanh (x);
 431 }
 432 #endif
 433
 434 #ifndef HAVE_TRUNC
 435 #define HAVE_TRUNC 1
 436 double trunc (double x);
 437
 438 double
 439 trunc (double x)
 440 {
 441   if (!isfinite (x))
 442     return x;
 443
 444   if (x < 0.0)
 445     return - floor (-x);
 446   else
 447     return floor (x);
 448 }
 449 #endif
 450
 451 #ifndef HAVE_TRUNCF
 452 #define HAVE_TRUNCF 1
 453 float truncf (float x);
 454
 455 float
 456 truncf (float x)
 457 {
 458   return (float) trunc (x);
 459 }
 460 #endif
 461
 462 #ifndef HAVE_NEXTAFTERF
 463 #define HAVE_NEXTAFTERF 1
 464 /* This is a portable implementation of nextafterf that is intended to be
 465    independent of the floating point format or its in memory representation.
 466    This implementation works correctly with denormalized values.  */
 467 float nextafterf (float x, float y);
 468
 469 float
 470 nextafterf (float x, float y)
 471 {
 472   /* This variable is marked volatile to avoid excess precision problems
 473      on some platforms, including IA-32.  */
 474   volatile float delta;
 475   float absx, denorm_min;
 476
 477   if (isnan (x) || isnan (y))
 478     return x + y;
 479   if (x == y)
 480     return x;
 481   if (!isfinite (x))
 482     return x > 0 ? __FLT_MAX__ : - __FLT_MAX__;
 483
 484   /* absx = fabsf (x);  */
 485   absx = (x < 0.0) ? -x : x;
 486
 487   /* __FLT_DENORM_MIN__ is non-zero iff the target supports denormals.  */
 488   if (__FLT_DENORM_MIN__ == 0.0f)
 489     denorm_min = __FLT_MIN__;
 490   else
 491     denorm_min = __FLT_DENORM_MIN__;
 492
 493   if (absx < __FLT_MIN__)
 494     delta = denorm_min;
 495   else
 496     {
 497       float frac;
 498       int exp;
 499
 500       /* Discard the fraction from x.  */
 501       frac = frexpf (absx, &exp);
 502       delta = scalbnf (0.5f, exp);
 503
 504       /* Scale x by the epsilon of the representation.  By rights we should
 505          have been able to combine this with scalbnf, but some targets don't
 506          get that correct with denormals.  */
 507       delta *= __FLT_EPSILON__;
 508
 509       /* If we're going to be reducing the absolute value of X, and doing so
 510          would reduce the exponent of X, then the delta to be applied is
 511          one exponent smaller.  */
 512       if (frac == 0.5f && (y < x) == (x > 0))
 513         delta *= 0.5f;
 514
 515       /* If that underflows to zero, then we're back to the minimum.  */
 516       if (delta == 0.0f)
 517         delta = denorm_min;
 518     }
 519
 520   if (y < x)
 521     delta = -delta;
 522
 523   return x + delta;
 524 }
 525 #endif
 526
 527
 528 #ifndef HAVE_POWF
 529 #define HAVE_POWF 1
 530 float powf (float x, float y);
 531
 532 float
 533 powf (float x, float y)
 534 {
 535   return (float) pow (x, y);
 536 }
 537 #endif
 538
 539
 540 #ifndef HAVE_ROUND
 541 #define HAVE_ROUND 1
 542 /* Round to nearest integral value.  If the argument is halfway between two
 543    integral values then round away from zero.  */
 544 double round (double x);
 545
 546 double
 547 round (double x)
 548 {
 549    double t;
 550    if (!isfinite (x))
 551      return (x);
 552
 553    if (x >= 0.0)
 554     {
 555       t = floor (x);
 556       if (t - x <= -0.5)
 557         t += 1.0;
 558       return (t);
 559     }
 560    else
 561     {
 562       t = floor (-x);
 563       if (t + x <= -0.5)
 564         t += 1.0;
 565       return (-t);
 566     }
 567 }
 568 #endif
 569
 570
 571 /* Algorithm by Steven G. Kargl.  */
 572
 573 #if !defined(HAVE_ROUNDL)
 574 #define HAVE_ROUNDL 1
 575 long double roundl (long double x);
 576
 577 #if defined(HAVE_CEILL)
 578 /* Round to nearest integral value.  If the argument is halfway between two
 579    integral values then round away from zero.  */
 580
 581 long double
 582 roundl (long double x)
 583 {
 584    long double t;
 585    if (!isfinite (x))
 586      return (x);
 587
 588    if (x >= 0.0)
 589     {
 590       t = ceill (x);
 591       if (t - x > 0.5)
 592         t -= 1.0;
 593       return (t);
 594     }
 595    else
 596     {
 597       t = ceill (-x);
 598       if (t + x > 0.5)
 599         t -= 1.0;
 600       return (-t);
 601     }
 602 }
 603 #else
 604
 605 /* Poor version of roundl for system that don't have ceill.  */
 606 long double
 607 roundl (long double x)
 608 {
 609   if (x > DBL_MAX || x < -DBL_MAX)
 610     {
 611 #ifdef HAVE_NEXTAFTERL
 612       long double prechalf = nextafterl (0.5L, LDBL_MAX);
 613 #else
 614       static long double prechalf = 0.5L;
 615 #endif
 616       return (GFC_INTEGER_LARGEST) (x + (x > 0 ? prechalf : -prechalf));
 617     }
 618   else
 619     /* Use round().  */
 620     return round ((double) x);
 621 }
 622
 623 #endif
 624 #endif
 625
 626 #ifndef HAVE_ROUNDF
 627 #define HAVE_ROUNDF 1
 628 /* Round to nearest integral value.  If the argument is halfway between two
 629    integral values then round away from zero.  */
 630 float roundf (float x);
 631
 632 float
 633 roundf (float x)
 634 {
 635    float t;
 636    if (!isfinite (x))
 637      return (x);
 638
 639    if (x >= 0.0)
 640     {
 641       t = floorf (x);
 642       if (t - x <= -0.5)
 643         t += 1.0;
 644       return (t);
 645     }
 646    else
 647     {
 648       t = floorf (-x);
 649       if (t + x <= -0.5)
 650         t += 1.0;
 651       return (-t);
 652     }
 653 }
 654 #endif
 655
 656
 657 /* lround{f,,l} and llround{f,,l} functions.  */
 658
 659 #if !defined(HAVE_LROUNDF) && defined(HAVE_ROUNDF)
 660 #define HAVE_LROUNDF 1
 661 long int lroundf (float x);
 662
 663 long int
 664 lroundf (float x)
 665 {
 666   return (long int) roundf (x);
 667 }
 668 #endif
 669
 670 #if !defined(HAVE_LROUND) && defined(HAVE_ROUND)
 671 #define HAVE_LROUND 1
 672 long int lround (double x);
 673
 674 long int
 675 lround (double x)
 676 {
 677   return (long int) round (x);
 678 }
 679 #endif
 680
 681 #if !defined(HAVE_LROUNDL) && defined(HAVE_ROUNDL)
 682 #define HAVE_LROUNDL 1
 683 long int lroundl (long double x);
 684
 685 long int
 686 lroundl (long double x)
 687 {
 688   return (long long int) roundl (x);
 689 }
 690 #endif
 691
 692 #if !defined(HAVE_LLROUNDF) && defined(HAVE_ROUNDF)
 693 #define HAVE_LLROUNDF 1
 694 long long int llroundf (float x);
 695
 696 long long int
 697 llroundf (float x)
 698 {
 699   return (long long int) roundf (x);
 700 }
 701 #endif
 702
 703 #if !defined(HAVE_LLROUND) && defined(HAVE_ROUND)
 704 #define HAVE_LLROUND 1
 705 long long int llround (double x);
 706
 707 long long int
 708 llround (double x)
 709 {
 710   return (long long int) round (x);
 711 }
 712 #endif
 713
 714 #if !defined(HAVE_LLROUNDL) && defined(HAVE_ROUNDL)
 715 #define HAVE_LLROUNDL 1
 716 long long int llroundl (long double x);
 717
 718 long long int
 719 llroundl (long double x)
 720 {
 721   return (long long int) roundl (x);
 722 }
 723 #endif
 724
 725
 726 #ifndef HAVE_LOG10L
 727 #define HAVE_LOG10L 1
 728 /* log10 function for long double variables. The version provided here
 729    reduces the argument until it fits into a double, then use log10.  */
 730 long double log10l (long double x);
 731
 732 long double
 733 log10l (long double x)
 734 {
 735 #if LDBL_MAX_EXP > DBL_MAX_EXP
 736   if (x > DBL_MAX)
 737     {
 738       double val;
 739       int p2_result = 0;
 740       if (x > 0x1p16383L) { p2_result += 16383; x /= 0x1p16383L; }
 741       if (x > 0x1p8191L) { p2_result += 8191; x /= 0x1p8191L; }
 742       if (x > 0x1p4095L) { p2_result += 4095; x /= 0x1p4095L; }
 743       if (x > 0x1p2047L) { p2_result += 2047; x /= 0x1p2047L; }
 744       if (x > 0x1p1023L) { p2_result += 1023; x /= 0x1p1023L; }
 745       val = log10 ((double) x);
 746       return (val + p2_result * .30102999566398119521373889472449302L);
 747     }
 748 #endif
 749 #if LDBL_MIN_EXP < DBL_MIN_EXP
 750   if (x < DBL_MIN)
 751     {
 752       double val;
 753       int p2_result = 0;
 754       if (x < 0x1p-16380L) { p2_result += 16380; x /= 0x1p-16380L; }
 755       if (x < 0x1p-8189L) { p2_result += 8189; x /= 0x1p-8189L; }
 756       if (x < 0x1p-4093L) { p2_result += 4093; x /= 0x1p-4093L; }
 757       if (x < 0x1p-2045L) { p2_result += 2045; x /= 0x1p-2045L; }
 758       if (x < 0x1p-1021L) { p2_result += 1021; x /= 0x1p-1021L; }
 759       val = fabs (log10 ((double) x));
 760       return (- val - p2_result * .30102999566398119521373889472449302L);
 761     }
 762 #endif
 763     return log10 (x);
 764 }
 765 #endif
 766
 767
 768 #ifndef HAVE_FLOORL
 769 #define HAVE_FLOORL 1
 770 long double floorl (long double x);
 771
 772 long double
 773 floorl (long double x)
 774 {
 775   /* Zero, possibly signed.  */
 776   if (x == 0)
 777     return x;
 778
 779   /* Large magnitude.  */
 780   if (x > DBL_MAX || x < (-DBL_MAX))
 781     return x;
 782
 783   /* Small positive values.  */
 784   if (x >= 0 && x < DBL_MIN)
 785     return 0;
 786
 787   /* Small negative values.  */
 788   if (x < 0 && x > (-DBL_MIN))
 789     return -1;
 790
 791   return floor (x);
 792 }
 793 #endif
 794
 795
 796 #ifndef HAVE_FMODL
 797 #define HAVE_FMODL 1
 798 long double fmodl (long double x, long double y);
 799
 800 long double
 801 fmodl (long double x, long double y)
 802 {
 803   if (y == 0.0L)
 804     return 0.0L;
 805
 806   /* Need to check that the result has the same sign as x and magnitude
 807      less than the magnitude of y.  */
 808   return x - floorl (x / y) * y;
 809 }
 810 #endif
 811
 812
 813 #if !defined(HAVE_CABSF)
 814 #define HAVE_CABSF 1
 815 float cabsf (float complex z);
 816
 817 float
 818 cabsf (float complex z)
 819 {
 820   return hypotf (REALPART (z), IMAGPART (z));
 821 }
 822 #endif
 823
 824 #if !defined(HAVE_CABS)
 825 #define HAVE_CABS 1
 826 double cabs (double complex z);
 827
 828 double
 829 cabs (double complex z)
 830 {
 831   return hypot (REALPART (z), IMAGPART (z));
 832 }
 833 #endif
 834
 835 #if !defined(HAVE_CABSL) && defined(HAVE_HYPOTL)
 836 #define HAVE_CABSL 1
 837 long double cabsl (long double complex z);
 838
 839 long double
 840 cabsl (long double complex z)
 841 {
 842   return hypotl (REALPART (z), IMAGPART (z));
 843 }
 844 #endif
 845
 846
 847 #if !defined(HAVE_CARGF)
 848 #define HAVE_CARGF 1
 849 float cargf (float complex z);
 850
 851 float
 852 cargf (float complex z)
 853 {
 854   return atan2f (IMAGPART (z), REALPART (z));
 855 }
 856 #endif
 857
 858 #if !defined(HAVE_CARG)
 859 #define HAVE_CARG 1
 860 double carg (double complex z);
 861
 862 double
 863 carg (double complex z)
 864 {
 865   return atan2 (IMAGPART (z), REALPART (z));
 866 }
 867 #endif
 868
 869 #if !defined(HAVE_CARGL) && defined(HAVE_ATAN2L)
 870 #define HAVE_CARGL 1
 871 long double cargl (long double complex z);
 872
 873 long double
 874 cargl (long double complex z)
 875 {
 876   return atan2l (IMAGPART (z), REALPART (z));
 877 }
 878 #endif
 879
 880
 881 /* exp(z) = exp(a)*(cos(b) + i sin(b))  */
 882 #if !defined(HAVE_CEXPF)
 883 #define HAVE_CEXPF 1
 884 float complex cexpf (float complex z);
 885
 886 float complex
 887 cexpf (float complex z)
 888 {
 889   float a, b;
 890   float complex v;
 891
 892   a = REALPART (z);
 893   b = IMAGPART (z);
 894   COMPLEX_ASSIGN (v, cosf (b), sinf (b));
 895   return expf (a) * v;
 896 }
 897 #endif
 898
 899 #if !defined(HAVE_CEXP)
 900 #define HAVE_CEXP 1
 901 double complex cexp (double complex z);
 902
 903 double complex
 904 cexp (double complex z)
 905 {
 906   double a, b;
 907   double complex v;
 908
 909   a = REALPART (z);
 910   b = IMAGPART (z);
 911   COMPLEX_ASSIGN (v, cos (b), sin (b));
 912   return exp (a) * v;
 913 }
 914 #endif
 915
 916 #if !defined(HAVE_CEXPL) && defined(HAVE_COSL) && defined(HAVE_SINL) && defined(EXPL)
 917 #define HAVE_CEXPL 1
 918 long double complex cexpl (long double complex z);
 919
 920 long double complex
 921 cexpl (long double complex z)
 922 {
 923   long double a, b;
 924   long double complex v;
 925
 926   a = REALPART (z);
 927   b = IMAGPART (z);
 928   COMPLEX_ASSIGN (v, cosl (b), sinl (b));
 929   return expl (a) * v;
 930 }
 931 #endif
 932
 933
 934 /* log(z) = log (cabs(z)) + i*carg(z)  */
 935 #if !defined(HAVE_CLOGF)
 936 #define HAVE_CLOGF 1
 937 float complex clogf (float complex z);
 938
 939 float complex
 940 clogf (float complex z)
 941 {
 942   float complex v;
 943
 944   COMPLEX_ASSIGN (v, logf (cabsf (z)), cargf (z));
 945   return v;
 946 }
 947 #endif
 948
 949 #if !defined(HAVE_CLOG)
 950 #define HAVE_CLOG 1
 951 double complex clog (double complex z);
 952
 953 double complex
 954 clog (double complex z)
 955 {
 956   double complex v;
 957
 958   COMPLEX_ASSIGN (v, log (cabs (z)), carg (z));
 959   return v;
 960 }
 961 #endif
 962
 963 #if !defined(HAVE_CLOGL) && defined(HAVE_LOGL) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
 964 #define HAVE_CLOGL 1
 965 long double complex clogl (long double complex z);
 966
 967 long double complex
 968 clogl (long double complex z)
 969 {
 970   long double complex v;
 971
 972   COMPLEX_ASSIGN (v, logl (cabsl (z)), cargl (z));
 973   return v;
 974 }
 975 #endif
 976
 977
 978 /* log10(z) = log10 (cabs(z)) + i*carg(z)  */
 979 #if !defined(HAVE_CLOG10F)
 980 #define HAVE_CLOG10F 1
 981 float complex clog10f (float complex z);
 982
 983 float complex
 984 clog10f (float complex z)
 985 {
 986   float complex v;
 987
 988   COMPLEX_ASSIGN (v, log10f (cabsf (z)), cargf (z));
 989   return v;
 990 }
 991 #endif
 992
 993 #if !defined(HAVE_CLOG10)
 994 #define HAVE_CLOG10 1
 995 double complex clog10 (double complex z);
 996
 997 double complex
 998 clog10 (double complex z)
 999 {
1000   double complex v;
1001
1002   COMPLEX_ASSIGN (v, log10 (cabs (z)), carg (z));
1003   return v;
1004 }
1005 #endif
1006
1007 #if !defined(HAVE_CLOG10L) && defined(HAVE_LOG10L) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
1008 #define HAVE_CLOG10L 1
1009 long double complex clog10l (long double complex z);
1010
1011 long double complex
1012 clog10l (long double complex z)
1013 {
1014   long double complex v;
1015
1016   COMPLEX_ASSIGN (v, log10l (cabsl (z)), cargl (z));
1017   return v;
1018 }
1019 #endif
1020
1021
1022 /* pow(base, power) = cexp (power * clog (base))  */
1023 #if !defined(HAVE_CPOWF)
1024 #define HAVE_CPOWF 1
1025 float complex cpowf (float complex base, float complex power);
1026
1027 float complex
1028 cpowf (float complex base, float complex power)
1029 {
1030   return cexpf (power * clogf (base));
1031 }
1032 #endif
1033
1034 #if !defined(HAVE_CPOW)
1035 #define HAVE_CPOW 1
1036 double complex cpow (double complex base, double complex power);
1037
1038 double complex
1039 cpow (double complex base, double complex power)
1040 {
1041   return cexp (power * clog (base));
1042 }
1043 #endif
1044
1045 #if !defined(HAVE_CPOWL) && defined(HAVE_CEXPL) && defined(HAVE_CLOGL)
1046 #define HAVE_CPOWL 1
1047 long double complex cpowl (long double complex base, long double complex power);
1048
1049 long double complex
1050 cpowl (long double complex base, long double complex power)
1051 {
1052   return cexpl (power * clogl (base));
1053 }
1054 #endif
1055
1056
1057 /* sqrt(z).  Algorithm pulled from glibc.  */
1058 #if !defined(HAVE_CSQRTF)
1059 #define HAVE_CSQRTF 1
1060 float complex csqrtf (float complex z);
1061
1062 float complex
1063 csqrtf (float complex z)
1064 {
1065   float re, im;
1066   float complex v;
1067
1068   re = REALPART (z);
1069   im = IMAGPART (z);
1070   if (im == 0)
1071     {
1072       if (re < 0)
1073         {
1074           COMPLEX_ASSIGN (v, 0, copysignf (sqrtf (-re), im));
1075         }
1076       else
1077         {
1078           COMPLEX_ASSIGN (v, fabsf (sqrtf (re)), copysignf (0, im));
1079         }
1080     }
1081   else if (re == 0)
1082     {
1083       float r;
1084
1085       r = sqrtf (0.5 * fabsf (im));
1086
1087       COMPLEX_ASSIGN (v, r, copysignf (r, im));
1088     }
1089   else
1090     {
1091       float d, r, s;
1092
1093       d = hypotf (re, im);
1094       /* Use the identity   2  Re res  Im res = Im x
1095          to avoid cancellation error in  d +/- Re x.  */
1096       if (re > 0)
1097         {
1098           r = sqrtf (0.5 * d + 0.5 * re);
1099           s = (0.5 * im) / r;
1100         }
1101       else
1102         {
1103           s = sqrtf (0.5 * d - 0.5 * re);
1104           r = fabsf ((0.5 * im) / s);
1105         }
1106
1107       COMPLEX_ASSIGN (v, r, copysignf (s, im));
1108     }
1109   return v;
1110 }
1111 #endif
1112
1113 #if !defined(HAVE_CSQRT)
1114 #define HAVE_CSQRT 1
1115 double complex csqrt (double complex z);
1116
1117 double complex
1118 csqrt (double complex z)
1119 {
1120   double re, im;
1121   double complex v;
1122
1123   re = REALPART (z);
1124   im = IMAGPART (z);
1125   if (im == 0)
1126     {
1127       if (re < 0)
1128         {
1129           COMPLEX_ASSIGN (v, 0, copysign (sqrt (-re), im));
1130         }
1131       else
1132         {
1133           COMPLEX_ASSIGN (v, fabs (sqrt (re)), copysign (0, im));
1134         }
1135     }
1136   else if (re == 0)
1137     {
1138       double r;
1139
1140       r = sqrt (0.5 * fabs (im));
1141
1142       COMPLEX_ASSIGN (v, r, copysign (r, im));
1143     }
1144   else
1145     {
1146       double d, r, s;
1147
1148       d = hypot (re, im);
1149       /* Use the identity   2  Re res  Im res = Im x
1150          to avoid cancellation error in  d +/- Re x.  */
1151       if (re > 0)
1152         {
1153           r = sqrt (0.5 * d + 0.5 * re);
1154           s = (0.5 * im) / r;
1155         }
1156       else
1157         {
1158           s = sqrt (0.5 * d - 0.5 * re);
1159           r = fabs ((0.5 * im) / s);
1160         }
1161
1162       COMPLEX_ASSIGN (v, r, copysign (s, im));
1163     }
1164   return v;
1165 }
1166 #endif
1167
1168 #if !defined(HAVE_CSQRTL) && defined(HAVE_COPYSIGNL) && defined(HAVE_SQRTL) && defined(HAVE_FABSL) && defined(HAVE_HYPOTL)
1169 #define HAVE_CSQRTL 1
1170 long double complex csqrtl (long double complex z);
1171
1172 long double complex
1173 csqrtl (long double complex z)
1174 {
1175   long double re, im;
1176   long double complex v;
1177
1178   re = REALPART (z);
1179   im = IMAGPART (z);
1180   if (im == 0)
1181     {
1182       if (re < 0)
1183         {
1184           COMPLEX_ASSIGN (v, 0, copysignl (sqrtl (-re), im));
1185         }
1186       else
1187         {
1188           COMPLEX_ASSIGN (v, fabsl (sqrtl (re)), copysignl (0, im));
1189         }
1190     }
1191   else if (re == 0)
1192     {
1193       long double r;
1194
1195       r = sqrtl (0.5 * fabsl (im));
1196
1197       COMPLEX_ASSIGN (v, copysignl (r, im), r);
1198     }
1199   else
1200     {
1201       long double d, r, s;
1202
1203       d = hypotl (re, im);
1204       /* Use the identity   2  Re res  Im res = Im x
1205          to avoid cancellation error in  d +/- Re x.  */
1206       if (re > 0)
1207         {
1208           r = sqrtl (0.5 * d + 0.5 * re);
1209           s = (0.5 * im) / r;
1210         }
1211       else
1212         {
1213           s = sqrtl (0.5 * d - 0.5 * re);
1214           r = fabsl ((0.5 * im) / s);
1215         }
1216
1217       COMPLEX_ASSIGN (v, r, copysignl (s, im));
1218     }
1219   return v;
1220 }
1221 #endif
1222
1223
1224 /* sinh(a + i b) = sinh(a) cos(b) + i cosh(a) sin(b)  */
1225 #if !defined(HAVE_CSINHF)
1226 #define HAVE_CSINHF 1
1227 float complex csinhf (float complex a);
1228
1229 float complex
1230 csinhf (float complex a)
1231 {
1232   float r, i;
1233   float complex v;
1234
1235   r = REALPART (a);
1236   i = IMAGPART (a);
1237   COMPLEX_ASSIGN (v, sinhf (r) * cosf (i), coshf (r) * sinf (i));
1238   return v;
1239 }
1240 #endif
1241
1242 #if !defined(HAVE_CSINH)
1243 #define HAVE_CSINH 1
1244 double complex csinh (double complex a);
1245
1246 double complex
1247 csinh (double complex a)
1248 {
1249   double r, i;
1250   double complex v;
1251
1252   r = REALPART (a);
1253   i = IMAGPART (a);
1254   COMPLEX_ASSIGN (v, sinh (r) * cos (i), cosh (r) * sin (i));
1255   return v;
1256 }
1257 #endif
1258
1259 #if !defined(HAVE_CSINHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1260 #define HAVE_CSINHL 1
1261 long double complex csinhl (long double complex a);
1262
1263 long double complex
1264 csinhl (long double complex a)
1265 {
1266   long double r, i;
1267   long double complex v;
1268
1269   r = REALPART (a);
1270   i = IMAGPART (a);
1271   COMPLEX_ASSIGN (v, sinhl (r) * cosl (i), coshl (r) * sinl (i));
1272   return v;
1273 }
1274 #endif
1275
1276
1277 /* cosh(a + i b) = cosh(a) cos(b) + i sinh(a) sin(b)  */
1278 #if !defined(HAVE_CCOSHF)
1279 #define HAVE_CCOSHF 1
1280 float complex ccoshf (float complex a);
1281
1282 float complex
1283 ccoshf (float complex a)
1284 {
1285   float r, i;
1286   float complex v;
1287
1288   r = REALPART (a);
1289   i = IMAGPART (a);
1290   COMPLEX_ASSIGN (v, coshf (r) * cosf (i), sinhf (r) * sinf (i));
1291   return v;
1292 }
1293 #endif
1294
1295 #if !defined(HAVE_CCOSH)
1296 #define HAVE_CCOSH 1
1297 double complex ccosh (double complex a);
1298
1299 double complex
1300 ccosh (double complex a)
1301 {
1302   double r, i;
1303   double complex v;
1304
1305   r = REALPART (a);
1306   i = IMAGPART (a);
1307   COMPLEX_ASSIGN (v, cosh (r) * cos (i),  sinh (r) * sin (i));
1308   return v;
1309 }
1310 #endif
1311
1312 #if !defined(HAVE_CCOSHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1313 #define HAVE_CCOSHL 1
1314 long double complex ccoshl (long double complex a);
1315
1316 long double complex
1317 ccoshl (long double complex a)
1318 {
1319   long double r, i;
1320   long double complex v;
1321
1322   r = REALPART (a);
1323   i = IMAGPART (a);
1324   COMPLEX_ASSIGN (v, coshl (r) * cosl (i), sinhl (r) * sinl (i));
1325   return v;
1326 }
1327 #endif
1328
1329
1330 /* tanh(a + i b) = (tanh(a) + i tan(b)) / (1 + i tanh(a) tan(b))  */
1331 #if !defined(HAVE_CTANHF)
1332 #define HAVE_CTANHF 1
1333 float complex ctanhf (float complex a);
1334
1335 float complex
1336 ctanhf (float complex a)
1337 {
1338   float rt, it;
1339   float complex n, d;
1340
1341   rt = tanhf (REALPART (a));
1342   it = tanf (IMAGPART (a));
1343   COMPLEX_ASSIGN (n, rt, it);
1344   COMPLEX_ASSIGN (d, 1, rt * it);
1345
1346   return n / d;
1347 }
1348 #endif
1349
1350 #if !defined(HAVE_CTANH)
1351 #define HAVE_CTANH 1
1352 double complex ctanh (double complex a);
1353 double complex
1354 ctanh (double complex a)
1355 {
1356   double rt, it;
1357   double complex n, d;
1358
1359   rt = tanh (REALPART (a));
1360   it = tan (IMAGPART (a));
1361   COMPLEX_ASSIGN (n, rt, it);
1362   COMPLEX_ASSIGN (d, 1, rt * it);
1363
1364   return n / d;
1365 }
1366 #endif
1367
1368 #if !defined(HAVE_CTANHL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1369 #define HAVE_CTANHL 1
1370 long double complex ctanhl (long double complex a);
1371
1372 long double complex
1373 ctanhl (long double complex a)
1374 {
1375   long double rt, it;
1376   long double complex n, d;
1377
1378   rt = tanhl (REALPART (a));
1379   it = tanl (IMAGPART (a));
1380   COMPLEX_ASSIGN (n, rt, it);
1381   COMPLEX_ASSIGN (d, 1, rt * it);
1382
1383   return n / d;
1384 }
1385 #endif
1386
1387
1388 /* sin(a + i b) = sin(a) cosh(b) + i cos(a) sinh(b)  */
1389 #if !defined(HAVE_CSINF)
1390 #define HAVE_CSINF 1
1391 float complex csinf (float complex a);
1392
1393 float complex
1394 csinf (float complex a)
1395 {
1396   float r, i;
1397   float complex v;
1398
1399   r = REALPART (a);
1400   i = IMAGPART (a);
1401   COMPLEX_ASSIGN (v, sinf (r) * coshf (i), cosf (r) * sinhf (i));
1402   return v;
1403 }
1404 #endif
1405
1406 #if !defined(HAVE_CSIN)
1407 #define HAVE_CSIN 1
1408 double complex csin (double complex a);
1409
1410 double complex
1411 csin (double complex a)
1412 {
1413   double r, i;
1414   double complex v;
1415
1416   r = REALPART (a);
1417   i = IMAGPART (a);
1418   COMPLEX_ASSIGN (v, sin (r) * cosh (i), cos (r) * sinh (i));
1419   return v;
1420 }
1421 #endif
1422
1423 #if !defined(HAVE_CSINL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1424 #define HAVE_CSINL 1
1425 long double complex csinl (long double complex a);
1426
1427 long double complex
1428 csinl (long double complex a)
1429 {
1430   long double r, i;
1431   long double complex v;
1432
1433   r = REALPART (a);
1434   i = IMAGPART (a);
1435   COMPLEX_ASSIGN (v, sinl (r) * coshl (i), cosl (r) * sinhl (i));
1436   return v;
1437 }
1438 #endif
1439
1440
1441 /* cos(a + i b) = cos(a) cosh(b) - i sin(a) sinh(b)  */
1442 #if !defined(HAVE_CCOSF)
1443 #define HAVE_CCOSF 1
1444 float complex ccosf (float complex a);
1445
1446 float complex
1447 ccosf (float complex a)
1448 {
1449   float r, i;
1450   float complex v;
1451
1452   r = REALPART (a);
1453   i = IMAGPART (a);
1454   COMPLEX_ASSIGN (v, cosf (r) * coshf (i), - (sinf (r) * sinhf (i)));
1455   return v;
1456 }
1457 #endif
1458
1459 #if !defined(HAVE_CCOS)
1460 #define HAVE_CCOS 1
1461 double complex ccos (double complex a);
1462
1463 double complex
1464 ccos (double complex a)
1465 {
1466   double r, i;
1467   double complex v;
1468
1469   r = REALPART (a);
1470   i = IMAGPART (a);
1471   COMPLEX_ASSIGN (v, cos (r) * cosh (i), - (sin (r) * sinh (i)));
1472   return v;
1473 }
1474 #endif
1475
1476 #if !defined(HAVE_CCOSL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1477 #define HAVE_CCOSL 1
1478 long double complex ccosl (long double complex a);
1479
1480 long double complex
1481 ccosl (long double complex a)
1482 {
1483   long double r, i;
1484   long double complex v;
1485
1486   r = REALPART (a);
1487   i = IMAGPART (a);
1488   COMPLEX_ASSIGN (v, cosl (r) * coshl (i), - (sinl (r) * sinhl (i)));
1489   return v;
1490 }
1491 #endif
1492
1493
1494 /* tan(a + i b) = (tan(a) + i tanh(b)) / (1 - i tan(a) tanh(b))  */
1495 #if !defined(HAVE_CTANF)
1496 #define HAVE_CTANF 1
1497 float complex ctanf (float complex a);
1498
1499 float complex
1500 ctanf (float complex a)
1501 {
1502   float rt, it;
1503   float complex n, d;
1504
1505   rt = tanf (REALPART (a));
1506   it = tanhf (IMAGPART (a));
1507   COMPLEX_ASSIGN (n, rt, it);
1508   COMPLEX_ASSIGN (d, 1, - (rt * it));
1509
1510   return n / d;
1511 }
1512 #endif
1513
1514 #if !defined(HAVE_CTAN)
1515 #define HAVE_CTAN 1
1516 double complex ctan (double complex a);
1517
1518 double complex
1519 ctan (double complex a)
1520 {
1521   double rt, it;
1522   double complex n, d;
1523
1524   rt = tan (REALPART (a));
1525   it = tanh (IMAGPART (a));
1526   COMPLEX_ASSIGN (n, rt, it);
1527   COMPLEX_ASSIGN (d, 1, - (rt * it));
1528
1529   return n / d;
1530 }
1531 #endif
1532
1533 #if !defined(HAVE_CTANL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1534 #define HAVE_CTANL 1
1535 long double complex ctanl (long double complex a);
1536
1537 long double complex
1538 ctanl (long double complex a)
1539 {
1540   long double rt, it;
1541   long double complex n, d;
1542
1543   rt = tanl (REALPART (a));
1544   it = tanhl (IMAGPART (a));
1545   COMPLEX_ASSIGN (n, rt, it);
1546   COMPLEX_ASSIGN (d, 1, - (rt * it));
1547
1548   return n / d;
1549 }
1550 #endif
1551
1552
1553 /* Complex ASIN.  Returns wrongly NaN for infinite arguments.
1554    Algorithm taken from Abramowitz & Stegun.  */
1555
1556 #if !defined(HAVE_CASINF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1557 #define HAVE_CASINF 1
1558 complex float casinf (complex float z);
1559
1560 complex float
1561 casinf (complex float z)
1562 {
1563   return -I*clogf (I*z + csqrtf (1.0f-z*z));
1564 }
1565 #endif
1566
1567
1568 #if !defined(HAVE_CASIN) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1569 #define HAVE_CASIN 1
1570 complex double casin (complex double z);
1571
1572 complex double
1573 casin (complex double z)
1574 {
1575   return -I*clog (I*z + csqrt (1.0-z*z));
1576 }
1577 #endif
1578
1579
1580 #if !defined(HAVE_CASINL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1581 #define HAVE_CASINL 1
1582 complex long double casinl (complex long double z);
1583
1584 complex long double
1585 casinl (complex long double z)
1586 {
1587   return -I*clogl (I*z + csqrtl (1.0L-z*z));
1588 }
1589 #endif
1590
1591
1592 /* Complex ACOS.  Returns wrongly NaN for infinite arguments.
1593    Algorithm taken from Abramowitz & Stegun.  */
1594
1595 #if !defined(HAVE_CACOSF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1596 #define HAVE_CACOSF 1
1597 complex float cacosf (complex float z);
1598
1599 complex float
1600 cacosf (complex float z)
1601 {
1602   return -I*clogf (z + I*csqrtf (1.0f-z*z));
1603 }
1604 #endif
1605
1606
1607 #if !defined(HAVE_CACOS) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1608 #define HAVE_CACOS 1
1609 complex double cacos (complex double z);
1610
1611 complex double
1612 cacos (complex double z)
1613 {
1614   return -I*clog (z + I*csqrt (1.0-z*z));
1615 }
1616 #endif
1617
1618
1619 #if !defined(HAVE_CACOSL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1620 #define HAVE_CACOSL 1
1621 complex long double cacosl (complex long double z);
1622
1623 complex long double
1624 cacosl (complex long double z)
1625 {
1626   return -I*clogl (z + I*csqrtl (1.0L-z*z));
1627 }
1628 #endif
1629
1630
1631 /* Complex ATAN.  Returns wrongly NaN for infinite arguments.
1632    Algorithm taken from Abramowitz & Stegun.  */
1633
1634 #if !defined(HAVE_CATANF) && defined(HAVE_CLOGF)
1635 #define HAVE_CACOSF 1
1636 complex float catanf (complex float z);
1637
1638 complex float
1639 catanf (complex float z)
1640 {
1641   return I*clogf ((I+z)/(I-z))/2.0f;
1642 }
1643 #endif
1644
1645
1646 #if !defined(HAVE_CATAN) && defined(HAVE_CLOG)
1647 #define HAVE_CACOS 1
1648 complex double catan (complex double z);
1649
1650 complex double
1651 catan (complex double z)
1652 {
1653   return I*clog ((I+z)/(I-z))/2.0;
1654 }
1655 #endif
1656
1657
1658 #if !defined(HAVE_CATANL) && defined(HAVE_CLOGL)
1659 #define HAVE_CACOSL 1
1660 complex long double catanl (complex long double z);
1661
1662 complex long double
1663 catanl (complex long double z)
1664 {
1665   return I*clogl ((I+z)/(I-z))/2.0L;
1666 }
1667 #endif
1668
1669
1670 /* Complex ASINH.  Returns wrongly NaN for infinite arguments.
1671    Algorithm taken from Abramowitz & Stegun.  */
1672
1673 #if !defined(HAVE_CASINHF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1674 #define HAVE_CASINHF 1
1675 complex float casinhf (complex float z);
1676
1677 complex float
1678 casinhf (complex float z)
1679 {
1680   return clogf (z + csqrtf (z*z+1.0f));
1681 }
1682 #endif
1683
1684
1685 #if !defined(HAVE_CASINH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1686 #define HAVE_CASINH 1
1687 complex double casinh (complex double z);
1688
1689 complex double
1690 casinh (complex double z)
1691 {
1692   return clog (z + csqrt (z*z+1.0));
1693 }
1694 #endif
1695
1696
1697 #if !defined(HAVE_CASINHL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1698 #define HAVE_CASINHL 1
1699 complex long double casinhl (complex long double z);
1700
1701 complex long double
1702 casinhl (complex long double z)
1703 {
1704   return clogl (z + csqrtl (z*z+1.0L));
1705 }
1706 #endif
1707
1708
1709 /* Complex ACOSH.  Returns wrongly NaN for infinite arguments.
1710    Algorithm taken from Abramowitz & Stegun.  */
1711
1712 #if !defined(HAVE_CACOSHF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1713 #define HAVE_CACOSHF 1
1714 complex float cacoshf (complex float z);
1715
1716 complex float
1717 cacoshf (complex float z)
1718 {
1719   return clogf (z + csqrtf (z-1.0f) * csqrtf (z+1.0f));
1720 }
1721 #endif
1722
1723
1724 #if !defined(HAVE_CACOSH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1725 #define HAVE_CACOSH 1
1726 complex double cacosh (complex double z);
1727
1728 complex double
1729 cacosh (complex double z)
1730 {
1731   return clog (z + csqrt (z-1.0) * csqrt (z+1.0));
1732 }
1733 #endif
1734
1735
1736 #if !defined(HAVE_CACOSHL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1737 #define HAVE_CACOSHL 1
1738 complex long double cacoshl (complex long double z);
1739
1740 complex long double
1741 cacoshl (complex long double z)
1742 {
1743   return clogl (z + csqrtl (z-1.0L) * csqrtl (z+1.0L));
1744 }
1745 #endif
1746
1747
1748 /* Complex ATANH.  Returns wrongly NaN for infinite arguments.
1749    Algorithm taken from Abramowitz & Stegun.  */
1750
1751 #if !defined(HAVE_CATANHF) && defined(HAVE_CLOGF)
1752 #define HAVE_CATANHF 1
1753 complex float catanhf (complex float z);
1754
1755 complex float
1756 catanhf (complex float z)
1757 {
1758   return clogf ((1.0f+z)/(1.0f-z))/2.0f;
1759 }
1760 #endif
1761
1762
1763 #if !defined(HAVE_CATANH) && defined(HAVE_CLOG)
1764 #define HAVE_CATANH 1
1765 complex double catanh (complex double z);
1766
1767 complex double
1768 catanh (complex double z)
1769 {
1770   return clog ((1.0+z)/(1.0-z))/2.0;
1771 }
1772 #endif
1773
1774 #if !defined(HAVE_CATANHL) && defined(HAVE_CLOGL)
1775 #define HAVE_CATANHL 1
1776 complex long double catanhl (complex long double z);
1777
1778 complex long double
1779 catanhl (complex long double z)
1780 {
1781   return clogl ((1.0L+z)/(1.0L-z))/2.0L;
1782 }
1783 #endif
1784
1785
1786 #if !defined(HAVE_TGAMMA)
1787 #define HAVE_TGAMMA 1
1788 double tgamma (double);
1789
1790 /* Fallback tgamma() function. Uses the algorithm from
1791    http://www.netlib.org/specfun/gamma and references therein.  */
1792
1793 #undef SQRTPI
1794 #define SQRTPI 0.9189385332046727417803297
1795
1796 #undef PI
1797 #define PI 3.1415926535897932384626434
1798
1799 double
1800 tgamma (double x)
1801 {
1802   int i, n, parity;
1803   double fact, res, sum, xden, xnum, y, y1, ysq, z;
1804
1805   static double p[8] = {
1806     -1.71618513886549492533811e0,  2.47656508055759199108314e1,
1807     -3.79804256470945635097577e2,  6.29331155312818442661052e2,
1808      8.66966202790413211295064e2, -3.14512729688483675254357e4,
1809     -3.61444134186911729807069e4,  6.64561438202405440627855e4 };
1810
1811   static double q[8] = {
1812     -3.08402300119738975254353e1,  3.15350626979604161529144e2,
1813     -1.01515636749021914166146e3, -3.10777167157231109440444e3,
1814      2.25381184209801510330112e4,  4.75584627752788110767815e3,
1815     -1.34659959864969306392456e5, -1.15132259675553483497211e5 };
1816
1817   static double c[7] = {             -1.910444077728e-03,
1818      8.4171387781295e-04,            -5.952379913043012e-04,
1819      7.93650793500350248e-04,        -2.777777777777681622553e-03,
1820      8.333333333333333331554247e-02,  5.7083835261e-03 };
1821
1822   static const double xminin = 2.23e-308;
1823   static const double xbig = 171.624;
1824   static const double xnan = __builtin_nan ("0x0"), xinf = __builtin_inf ();
1825   static double eps = 0;
1826
1827   if (eps == 0)
1828     eps = nextafter (1., 2.) - 1.;
1829
1830   parity = 0;
1831   fact = 1;
1832   n = 0;
1833   y = x;
1834
1835   if (isnan (x))
1836     return x;
1837
1838   if (y <= 0)
1839     {
1840       y = -x;
1841       y1 = trunc (y);
1842       res = y - y1;
1843
1844       if (res != 0)
1845         {
1846           if (y1 != trunc (y1*0.5l)*2)
1847             parity = 1;
1848           fact = -PI / sin (PI*res);
1849           y = y + 1;
1850         }
1851       else
1852         return x == 0 ? copysign (xinf, x) : xnan;
1853     }
1854
1855   if (y < eps)
1856     {
1857       if (y >= xminin)
1858         res = 1 / y;
1859       else
1860         return xinf;
1861     }
1862   else if (y < 13)
1863     {
1864       y1 = y;
1865       if (y < 1)
1866         {
1867           z = y;
1868           y = y + 1;
1869         }
1870       else
1871         {
1872           n = (int)y - 1;
1873           y = y - n;
1874           z = y - 1;
1875         }
1876
1877       xnum = 0;
1878       xden = 1;
1879       for (i = 0; i < 8; i++)
1880         {
1881           xnum = (xnum + p[i]) * z;
1882           xden = xden * z + q[i];
1883         }
1884
1885       res = xnum / xden + 1;
1886
1887       if (y1 < y)
1888         res = res / y1;
1889       else if (y1 > y)
1890         for (i = 1; i <= n; i++)
1891           {
1892             res = res * y;
1893             y = y + 1;
1894           }
1895     }
1896   else
1897     {
1898       if (y < xbig)
1899         {
1900           ysq = y * y;
1901           sum = c[6];
1902           for (i = 0; i < 6; i++)
1903             sum = sum / ysq + c[i];
1904
1905           sum = sum/y - y + SQRTPI;
1906           sum = sum + (y - 0.5) * log (y);
1907           res = exp (sum);
1908         }
1909       else
1910         return x < 0 ? xnan : xinf;
1911     }
1912
1913   if (parity)
1914     res = -res;
1915   if (fact != 1)
1916     res = fact / res;
1917
1918   return res;
1919 }
1920 #endif
1921
1922
1923
1924 #if !defined(HAVE_LGAMMA)
1925 #define HAVE_LGAMMA 1
1926 double lgamma (double);
1927
1928 /* Fallback lgamma() function. Uses the algorithm from
1929    http://www.netlib.org/specfun/algama and references therein,
1930    except for negative arguments (where netlib would return +Inf)
1931    where we use the following identity:
1932        lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y)
1933  */
1934
1935 double
1936 lgamma (double y)
1937 {
1938
1939 #undef SQRTPI
1940 #define SQRTPI 0.9189385332046727417803297
1941
1942 #undef PI
1943 #define PI 3.1415926535897932384626434
1944
1945 #define PNT68  0.6796875
1946 #define D1    -0.5772156649015328605195174
1947 #define D2     0.4227843350984671393993777
1948 #define D4     1.791759469228055000094023
1949
1950   static double p1[8] = {
1951               4.945235359296727046734888e0, 2.018112620856775083915565e2,
1952               2.290838373831346393026739e3, 1.131967205903380828685045e4,
1953               2.855724635671635335736389e4, 3.848496228443793359990269e4,
1954               2.637748787624195437963534e4, 7.225813979700288197698961e3 };
1955   static double q1[8] = {
1956               6.748212550303777196073036e1,  1.113332393857199323513008e3,
1957               7.738757056935398733233834e3,  2.763987074403340708898585e4,
1958               5.499310206226157329794414e4,  6.161122180066002127833352e4,
1959               3.635127591501940507276287e4,  8.785536302431013170870835e3 };
1960   static double p2[8] = {
1961               4.974607845568932035012064e0,  5.424138599891070494101986e2,
1962               1.550693864978364947665077e4,  1.847932904445632425417223e5,
1963               1.088204769468828767498470e6,  3.338152967987029735917223e6,
1964               5.106661678927352456275255e6,  3.074109054850539556250927e6 };
1965   static double q2[8] = {
1966               1.830328399370592604055942e2,  7.765049321445005871323047e3,
1967               1.331903827966074194402448e5,  1.136705821321969608938755e6,
1968               5.267964117437946917577538e6,  1.346701454311101692290052e7,
1969               1.782736530353274213975932e7,  9.533095591844353613395747e6 };
1970   static double p4[8] = {
1971               1.474502166059939948905062e4,  2.426813369486704502836312e6,
1972               1.214755574045093227939592e8,  2.663432449630976949898078e9,
1973               2.940378956634553899906876e10, 1.702665737765398868392998e11,
1974               4.926125793377430887588120e11, 5.606251856223951465078242e11 };
1975   static double q4[8] = {
1976               2.690530175870899333379843e3,  6.393885654300092398984238e5,
1977               4.135599930241388052042842e7,  1.120872109616147941376570e9,
1978               1.488613728678813811542398e10, 1.016803586272438228077304e11,
1979               3.417476345507377132798597e11, 4.463158187419713286462081e11 };
1980   static double  c[7] = {
1981              -1.910444077728e-03,            8.4171387781295e-04,
1982              -5.952379913043012e-04,         7.93650793500350248e-04,
1983              -2.777777777777681622553e-03,   8.333333333333333331554247e-02,
1984               5.7083835261e-03 };
1985
1986   static double xbig = 2.55e305, xinf = __builtin_inf (), eps = 0,
1987                 frtbig = 2.25e76;
1988
1989   int i;
1990   double corr, res, xden, xm1, xm2, xm4, xnum, ysq;
1991
1992   if (eps == 0)
1993     eps = __builtin_nextafter (1., 2.) - 1.;
1994
1995   if ((y > 0) && (y <= xbig))
1996     {
1997       if (y <= eps)
1998         res = -log (y);
1999       else if (y <= 1.5)
2000         {
2001           if (y < PNT68)
2002             {
2003               corr = -log (y);
2004               xm1 = y;
2005             }
2006           else
2007             {
2008               corr = 0;
2009               xm1 = (y - 0.5) - 0.5;
2010             }
2011
2012           if ((y <= 0.5) || (y >= PNT68))
2013             {
2014               xden = 1;
2015               xnum = 0;
2016               for (i = 0; i < 8; i++)
2017                 {
2018                   xnum = xnum*xm1 + p1[i];
2019                   xden = xden*xm1 + q1[i];
2020                 }
2021               res = corr + (xm1 * (D1 + xm1*(xnum/xden)));
2022             }
2023           else
2024             {
2025               xm2 = (y - 0.5) - 0.5;
2026               xden = 1;
2027               xnum = 0;
2028               for (i = 0; i < 8; i++)
2029                 {
2030                   xnum = xnum*xm2 + p2[i];
2031                   xden = xden*xm2 + q2[i];
2032                 }
2033               res = corr + xm2 * (D2 + xm2*(xnum/xden));
2034             }
2035         }
2036       else if (y <= 4)
2037         {
2038           xm2 = y - 2;
2039           xden = 1;
2040           xnum = 0;
2041           for (i = 0; i < 8; i++)
2042             {
2043               xnum = xnum*xm2 + p2[i];
2044               xden = xden*xm2 + q2[i];
2045             }
2046           res = xm2 * (D2 + xm2*(xnum/xden));
2047         }
2048       else if (y <= 12)
2049         {
2050           xm4 = y - 4;
2051           xden = -1;
2052           xnum = 0;
2053           for (i = 0; i < 8; i++)
2054             {
2055               xnum = xnum*xm4 + p4[i];
2056               xden = xden*xm4 + q4[i];
2057             }
2058           res = D4 + xm4*(xnum/xden);
2059         }
2060       else
2061         {
2062           res = 0;
2063           if (y <= frtbig)
2064             {
2065               res = c[6];
2066               ysq = y * y;
2067               for (i = 0; i < 6; i++)
2068                 res = res / ysq + c[i];
2069             }
2070           res = res/y;
2071           corr = log (y);
2072           res = res + SQRTPI - 0.5*corr;
2073           res = res + y*(corr-1);
2074         }
2075     }
2076   else if (y < 0 && __builtin_floor (y) != y)
2077     {
2078       /* lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y)
2079          For abs(y) very close to zero, we use a series expansion to
2080          the first order in y to avoid overflow.  */
2081       if (y > -1.e-100)
2082         res = -2 * log (fabs (y)) - lgamma (-y);
2083       else
2084         res = log (PI / fabs (y * sin (PI * y))) - lgamma (-y);
2085     }
2086   else
2087     res = xinf;
2088
2089   return res;
2090 }
2091 #endif
2092
2093
2094 #if defined(HAVE_TGAMMA) && !defined(HAVE_TGAMMAF)
2095 #define HAVE_TGAMMAF 1
2096 float tgammaf (float);
2097
2098 float
2099 tgammaf (float x)
2100 {
2101   return (float) tgamma ((double) x);
2102 }
2103 #endif
2104
2105 #if defined(HAVE_LGAMMA) && !defined(HAVE_LGAMMAF)
2106 #define HAVE_LGAMMAF 1
2107 float lgammaf (float);
2108
2109 float
2110 lgammaf (float x)
2111 {
2112   return (float) lgamma ((double) x);
2113 }
2114 #endif