libgfortran/intrinsics/c99_functions.c

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