libgfortran/intrinsics/c99_functions.c

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