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 #if !defined(HAVE_TGAMMA)
1416 #define HAVE_TGAMMA 1
1417
1418 extern double tgamma (double);
1419
1420 /* Fallback tgamma() function. Uses the algorithm from
1421    http://www.netlib.org/specfun/gamma and references therein.  */
1422
1423 #undef SQRTPI
1424 #define SQRTPI 0.9189385332046727417803297
1425
1426 #undef PI
1427 #define PI 3.1415926535897932384626434
1428
1429 double
1430 tgamma (double x)
1431 {
1432   int i, n, parity;
1433   double fact, res, sum, xden, xnum, y, y1, ysq, z;
1434
1435   static double p[8] = {
1436     -1.71618513886549492533811e0,  2.47656508055759199108314e1,
1437     -3.79804256470945635097577e2,  6.29331155312818442661052e2,
1438      8.66966202790413211295064e2, -3.14512729688483675254357e4,
1439     -3.61444134186911729807069e4,  6.64561438202405440627855e4 };
1440
1441   static double q[8] = {
1442     -3.08402300119738975254353e1,  3.15350626979604161529144e2,
1443     -1.01515636749021914166146e3, -3.10777167157231109440444e3,
1444      2.25381184209801510330112e4,  4.75584627752788110767815e3,
1445     -1.34659959864969306392456e5, -1.15132259675553483497211e5 };
1446
1447   static double c[7] = {             -1.910444077728e-03,
1448      8.4171387781295e-04,            -5.952379913043012e-04,
1449      7.93650793500350248e-04,        -2.777777777777681622553e-03,
1450      8.333333333333333331554247e-02,  5.7083835261e-03 };
1451
1452   static const double xminin = 2.23e-308;
1453   static const double xbig = 171.624;
1454   static const double xnan = __builtin_nan ("0x0"), xinf = __builtin_inf ();
1455   static double eps = 0;
1456
1457   if (eps == 0)
1458     eps = nextafter(1., 2.) - 1.;
1459
1460   parity = 0;
1461   fact = 1;
1462   n = 0;
1463   y = x;
1464
1465   if (__builtin_isnan (x))
1466     return x;
1467
1468   if (y <= 0)
1469     {
1470       y = -x;
1471       y1 = trunc(y);
1472       res = y - y1;
1473
1474       if (res != 0)
1475         {
1476           if (y1 != trunc(y1*0.5l)*2)
1477             parity = 1;
1478           fact = -PI / sin(PI*res);
1479           y = y + 1;
1480         }
1481       else
1482         return x == 0 ? copysign (xinf, x) : xnan;
1483     }
1484
1485   if (y < eps)
1486     {
1487       if (y >= xminin)
1488         res = 1 / y;
1489       else
1490         return xinf;
1491     }
1492   else if (y < 13)
1493     {
1494       y1 = y;
1495       if (y < 1)
1496         {
1497           z = y;
1498           y = y + 1;
1499         }
1500       else
1501         {
1502           n = (int)y - 1;
1503           y = y - n;
1504           z = y - 1;
1505         }
1506
1507       xnum = 0;
1508       xden = 1;
1509       for (i = 0; i < 8; i++)
1510         {
1511           xnum = (xnum + p[i]) * z;
1512           xden = xden * z + q[i];
1513         }
1514
1515       res = xnum / xden + 1;
1516
1517       if (y1 < y)
1518         res = res / y1;
1519       else if (y1 > y)
1520         for (i = 1; i <= n; i++)
1521           {
1522             res = res * y;
1523             y = y + 1;
1524           }
1525     }
1526   else
1527     {
1528       if (y < xbig)
1529         {
1530           ysq = y * y;
1531           sum = c[6];
1532           for (i = 0; i < 6; i++)
1533             sum = sum / ysq + c[i];
1534
1535           sum = sum/y - y + SQRTPI;
1536           sum = sum + (y - 0.5) * log(y);
1537           res = exp(sum);
1538         }
1539       else
1540         return x < 0 ? xnan : xinf;
1541     }
1542
1543   if (parity)
1544     res = -res;
1545   if (fact != 1)
1546     res = fact / res;
1547
1548   return res;
1549 }
1550 #endif
1551
1552
1553
1554 #if !defined(HAVE_LGAMMA)
1555 #define HAVE_LGAMMA 1
1556
1557 extern double lgamma (double);
1558
1559 /* Fallback lgamma() function. Uses the algorithm from
1560    http://www.netlib.org/specfun/algama and references therein,
1561    except for negative arguments (where netlib would return +Inf)
1562    where we use the following identity:
1563        lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y)
1564  */
1565
1566 double
1567 lgamma (double y)
1568 {
1569
1570 #undef SQRTPI
1571 #define SQRTPI 0.9189385332046727417803297
1572
1573 #undef PI
1574 #define PI 3.1415926535897932384626434
1575
1576 #define PNT68  0.6796875
1577 #define D1    -0.5772156649015328605195174
1578 #define D2     0.4227843350984671393993777
1579 #define D4     1.791759469228055000094023
1580
1581   static double p1[8] = {
1582               4.945235359296727046734888e0, 2.018112620856775083915565e2,
1583               2.290838373831346393026739e3, 1.131967205903380828685045e4,
1584               2.855724635671635335736389e4, 3.848496228443793359990269e4,
1585               2.637748787624195437963534e4, 7.225813979700288197698961e3 };
1586   static double q1[8] = {
1587               6.748212550303777196073036e1,  1.113332393857199323513008e3,
1588               7.738757056935398733233834e3,  2.763987074403340708898585e4,
1589               5.499310206226157329794414e4,  6.161122180066002127833352e4,
1590               3.635127591501940507276287e4,  8.785536302431013170870835e3 };
1591   static double p2[8] = {
1592               4.974607845568932035012064e0,  5.424138599891070494101986e2,
1593               1.550693864978364947665077e4,  1.847932904445632425417223e5,
1594               1.088204769468828767498470e6,  3.338152967987029735917223e6,
1595               5.106661678927352456275255e6,  3.074109054850539556250927e6 };
1596   static double q2[8] = {
1597               1.830328399370592604055942e2,  7.765049321445005871323047e3,
1598               1.331903827966074194402448e5,  1.136705821321969608938755e6,
1599               5.267964117437946917577538e6,  1.346701454311101692290052e7,
1600               1.782736530353274213975932e7,  9.533095591844353613395747e6 };
1601   static double p4[8] = {
1602               1.474502166059939948905062e4,  2.426813369486704502836312e6,
1603               1.214755574045093227939592e8,  2.663432449630976949898078e9,
1604               2.940378956634553899906876e10, 1.702665737765398868392998e11,
1605               4.926125793377430887588120e11, 5.606251856223951465078242e11 };
1606   static double q4[8] = {
1607               2.690530175870899333379843e3,  6.393885654300092398984238e5,
1608               4.135599930241388052042842e7,  1.120872109616147941376570e9,
1609               1.488613728678813811542398e10, 1.016803586272438228077304e11,
1610               3.417476345507377132798597e11, 4.463158187419713286462081e11 };
1611   static double  c[7] = {
1612              -1.910444077728e-03,            8.4171387781295e-04,
1613              -5.952379913043012e-04,         7.93650793500350248e-04,
1614              -2.777777777777681622553e-03,   8.333333333333333331554247e-02,
1615               5.7083835261e-03 };
1616
1617   static double xbig = 2.55e305, xinf = __builtin_inf (), eps = 0,
1618                 frtbig = 2.25e76;
1619
1620   int i;
1621   double corr, res, xden, xm1, xm2, xm4, xnum, ysq;
1622
1623   if (eps == 0)
1624     eps = __builtin_nextafter(1., 2.) - 1.;
1625
1626   if ((y > 0) && (y <= xbig))
1627     {
1628       if (y <= eps)
1629         res = -log(y);
1630       else if (y <= 1.5)
1631         {
1632           if (y < PNT68)
1633             {
1634               corr = -log(y);
1635               xm1 = y;
1636             }
1637           else
1638             {
1639               corr = 0;
1640               xm1 = (y - 0.5) - 0.5;
1641             }
1642
1643           if ((y <= 0.5) || (y >= PNT68))
1644             {
1645               xden = 1;
1646               xnum = 0;
1647               for (i = 0; i < 8; i++)
1648                 {
1649                   xnum = xnum*xm1 + p1[i];
1650                   xden = xden*xm1 + q1[i];
1651                 }
1652               res = corr + (xm1 * (D1 + xm1*(xnum/xden)));
1653             }
1654           else
1655             {
1656               xm2 = (y - 0.5) - 0.5;
1657               xden = 1;
1658               xnum = 0;
1659               for (i = 0; i < 8; i++)
1660                 {
1661                   xnum = xnum*xm2 + p2[i];
1662                   xden = xden*xm2 + q2[i];
1663                 }
1664               res = corr + xm2 * (D2 + xm2*(xnum/xden));
1665             }
1666         }
1667       else if (y <= 4)
1668         {
1669           xm2 = y - 2;
1670           xden = 1;
1671           xnum = 0;
1672           for (i = 0; i < 8; i++)
1673             {
1674               xnum = xnum*xm2 + p2[i];
1675               xden = xden*xm2 + q2[i];
1676             }
1677           res = xm2 * (D2 + xm2*(xnum/xden));
1678         }
1679       else if (y <= 12)
1680         {
1681           xm4 = y - 4;
1682           xden = -1;
1683           xnum = 0;
1684           for (i = 0; i < 8; i++)
1685             {
1686               xnum = xnum*xm4 + p4[i];
1687               xden = xden*xm4 + q4[i];
1688             }
1689           res = D4 + xm4*(xnum/xden);
1690         }
1691       else
1692         {
1693           res = 0;
1694           if (y <= frtbig)
1695             {
1696               res = c[6];
1697               ysq = y * y;
1698               for (i = 0; i < 6; i++)
1699                 res = res / ysq + c[i];
1700             }
1701           res = res/y;
1702           corr = log(y);
1703           res = res + SQRTPI - 0.5*corr;
1704           res = res + y*(corr-1);
1705         }
1706     }
1707   else if (y < 0 && __builtin_floor (y) != y)
1708     {
1709       /* lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y)
1710          For abs(y) very close to zero, we use a series expansion to
1711          the first order in y to avoid overflow.  */
1712       if (y > -1.e-100)
1713         res = -2 * log (fabs (y)) - lgamma (-y);
1714       else
1715         res = log (PI / fabs (y * sin (PI * y))) - lgamma (-y);
1716     }
1717   else
1718     res = xinf;
1719
1720   return res;
1721 }
1722 #endif
1723
1724
1725 #if defined(HAVE_TGAMMA) && !defined(HAVE_TGAMMAF)
1726 #define HAVE_TGAMMAF 1
1727 extern float tgammaf (float);
1728
1729 float
1730 tgammaf (float x)
1731 {
1732   return (float) tgamma ((double) x);
1733 }
1734 #endif
1735
1736 #if defined(HAVE_LGAMMA) && !defined(HAVE_LGAMMAF)
1737 #define HAVE_LGAMMAF 1
1738 extern float lgammaf (float);
1739
1740 float
1741 lgammaf (float x)
1742 {
1743   return (float) lgamma ((double) x);
1744 }
1745 #endif