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