libgfortran/intrinsics/c99_functions.c

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