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 /* Note that if fpclassify is not defined, then NaN is not handled */
 562
 563 /* Algorithm by Steven G. Kargl.  */
 564
 565 #if !defined(HAVE_ROUNDL)
 566 #define HAVE_ROUNDL 1
 567 long double roundl (long double x);
 568
 569 #if defined(HAVE_CEILL)
 570 /* Round to nearest integral value.  If the argument is halfway between two
 571    integral values then round away from zero.  */
 572
 573 long double
 574 roundl (long double x)
 575 {
 576    long double t;
 577    if (!isfinite (x))
 578      return (x);
 579
 580    if (x >= 0.0)
 581     {
 582       t = ceill (x);
 583       if (t - x > 0.5)
 584         t -= 1.0;
 585       return (t);
 586     }
 587    else
 588     {
 589       t = ceill (-x);
 590       if (t + x > 0.5)
 591         t -= 1.0;
 592       return (-t);
 593     }
 594 }
 595 #else
 596
 597 /* Poor version of roundl for system that don't have ceill.  */
 598 long double
 599 roundl (long double x)
 600 {
 601   if (x > DBL_MAX || x < -DBL_MAX)
 602     {
 603 #ifdef HAVE_NEXTAFTERL
 604       static long double prechalf = nexafterl (0.5L, LDBL_MAX);
 605 #else
 606       static long double prechalf = 0.5L;
 607 #endif
 608       return (GFC_INTEGER_LARGEST) (x + (x > 0 ? prechalf : -prechalf));
 609     }
 610   else
 611     /* Use round().  */
 612     return round ((double) x);
 613 }
 614
 615 #endif
 616 #endif
 617
 618 #ifndef HAVE_ROUND
 619 #define HAVE_ROUND 1
 620 /* Round to nearest integral value.  If the argument is halfway between two
 621    integral values then round away from zero.  */
 622 double round (double x);
 623
 624 double
 625 round (double x)
 626 {
 627    double t;
 628    if (!isfinite (x))
 629      return (x);
 630
 631    if (x >= 0.0)
 632     {
 633       t = floor (x);
 634       if (t - x <= -0.5)
 635         t += 1.0;
 636       return (t);
 637     }
 638    else
 639     {
 640       t = floor (-x);
 641       if (t + x <= -0.5)
 642         t += 1.0;
 643       return (-t);
 644     }
 645 }
 646 #endif
 647
 648 #ifndef HAVE_ROUNDF
 649 #define HAVE_ROUNDF 1
 650 /* Round to nearest integral value.  If the argument is halfway between two
 651    integral values then round away from zero.  */
 652 float roundf (float x);
 653
 654 float
 655 roundf (float x)
 656 {
 657    float t;
 658    if (!isfinite (x))
 659      return (x);
 660
 661    if (x >= 0.0)
 662     {
 663       t = floorf (x);
 664       if (t - x <= -0.5)
 665         t += 1.0;
 666       return (t);
 667     }
 668    else
 669     {
 670       t = floorf (-x);
 671       if (t + x <= -0.5)
 672         t += 1.0;
 673       return (-t);
 674     }
 675 }
 676 #endif
 677
 678
 679 /* lround{f,,l} and llround{f,,l} functions.  */
 680
 681 #if !defined(HAVE_LROUNDF) && defined(HAVE_ROUNDF)
 682 #define HAVE_LROUNDF 1
 683 long int lroundf (float x);
 684
 685 long int
 686 lroundf (float x)
 687 {
 688   return (long int) roundf (x);
 689 }
 690 #endif
 691
 692 #if !defined(HAVE_LROUND) && defined(HAVE_ROUND)
 693 #define HAVE_LROUND 1
 694 long int lround (double x);
 695
 696 long int
 697 lround (double x)
 698 {
 699   return (long int) round (x);
 700 }
 701 #endif
 702
 703 #if !defined(HAVE_LROUNDL) && defined(HAVE_ROUNDL)
 704 #define HAVE_LROUNDL 1
 705 long int lroundl (long double x);
 706
 707 long int
 708 lroundl (long double x)
 709 {
 710   return (long long int) roundl (x);
 711 }
 712 #endif
 713
 714 #if !defined(HAVE_LLROUNDF) && defined(HAVE_ROUNDF)
 715 #define HAVE_LLROUNDF 1
 716 long long int llroundf (float x);
 717
 718 long long int
 719 llroundf (float x)
 720 {
 721   return (long long int) roundf (x);
 722 }
 723 #endif
 724
 725 #if !defined(HAVE_LLROUND) && defined(HAVE_ROUND)
 726 #define HAVE_LLROUND 1
 727 long long int llround (double x);
 728
 729 long long int
 730 llround (double x)
 731 {
 732   return (long long int) round (x);
 733 }
 734 #endif
 735
 736 #if !defined(HAVE_LLROUNDL) && defined(HAVE_ROUNDL)
 737 #define HAVE_LLROUNDL 1
 738 long long int llroundl (long double x);
 739
 740 long long int
 741 llroundl (long double x)
 742 {
 743   return (long long int) roundl (x);
 744 }
 745 #endif
 746
 747
 748 #ifndef HAVE_LOG10L
 749 #define HAVE_LOG10L 1
 750 /* log10 function for long double variables. The version provided here
 751    reduces the argument until it fits into a double, then use log10.  */
 752 long double log10l (long double x);
 753
 754 long double
 755 log10l (long double x)
 756 {
 757 #if LDBL_MAX_EXP > DBL_MAX_EXP
 758   if (x > DBL_MAX)
 759     {
 760       double val;
 761       int p2_result = 0;
 762       if (x > 0x1p16383L) { p2_result += 16383; x /= 0x1p16383L; }
 763       if (x > 0x1p8191L) { p2_result += 8191; x /= 0x1p8191L; }
 764       if (x > 0x1p4095L) { p2_result += 4095; x /= 0x1p4095L; }
 765       if (x > 0x1p2047L) { p2_result += 2047; x /= 0x1p2047L; }
 766       if (x > 0x1p1023L) { p2_result += 1023; x /= 0x1p1023L; }
 767       val = log10 ((double) x);
 768       return (val + p2_result * .30102999566398119521373889472449302L);
 769     }
 770 #endif
 771 #if LDBL_MIN_EXP < DBL_MIN_EXP
 772   if (x < DBL_MIN)
 773     {
 774       double val;
 775       int p2_result = 0;
 776       if (x < 0x1p-16380L) { p2_result += 16380; x /= 0x1p-16380L; }
 777       if (x < 0x1p-8189L) { p2_result += 8189; x /= 0x1p-8189L; }
 778       if (x < 0x1p-4093L) { p2_result += 4093; x /= 0x1p-4093L; }
 779       if (x < 0x1p-2045L) { p2_result += 2045; x /= 0x1p-2045L; }
 780       if (x < 0x1p-1021L) { p2_result += 1021; x /= 0x1p-1021L; }
 781       val = fabs (log10 ((double) x));
 782       return (- val - p2_result * .30102999566398119521373889472449302L);
 783     }
 784 #endif
 785     return log10 (x);
 786 }
 787 #endif
 788
 789
 790 #ifndef HAVE_FLOORL
 791 #define HAVE_FLOORL 1
 792 long double floorl (long double x);
 793
 794 long double
 795 floorl (long double x)
 796 {
 797   /* Zero, possibly signed.  */
 798   if (x == 0)
 799     return x;
 800
 801   /* Large magnitude.  */
 802   if (x > DBL_MAX || x < (-DBL_MAX))
 803     return x;
 804
 805   /* Small positive values.  */
 806   if (x >= 0 && x < DBL_MIN)
 807     return 0;
 808
 809   /* Small negative values.  */
 810   if (x < 0 && x > (-DBL_MIN))
 811     return -1;
 812
 813   return floor (x);
 814 }
 815 #endif
 816
 817
 818 #ifndef HAVE_FMODL
 819 #define HAVE_FMODL 1
 820 long double fmodl (long double x, long double y);
 821
 822 long double
 823 fmodl (long double x, long double y)
 824 {
 825   if (y == 0.0L)
 826     return 0.0L;
 827
 828   /* Need to check that the result has the same sign as x and magnitude
 829      less than the magnitude of y.  */
 830   return x - floorl (x / y) * y;
 831 }
 832 #endif
 833
 834
 835 #if !defined(HAVE_CABSF)
 836 #define HAVE_CABSF 1
 837 float cabsf (float complex z);
 838
 839 float
 840 cabsf (float complex z)
 841 {
 842   return hypotf (REALPART (z), IMAGPART (z));
 843 }
 844 #endif
 845
 846 #if !defined(HAVE_CABS)
 847 #define HAVE_CABS 1
 848 double cabs (double complex z);
 849
 850 double
 851 cabs (double complex z)
 852 {
 853   return hypot (REALPART (z), IMAGPART (z));
 854 }
 855 #endif
 856
 857 #if !defined(HAVE_CABSL) && defined(HAVE_HYPOTL)
 858 #define HAVE_CABSL 1
 859 long double cabsl (long double complex z);
 860
 861 long double
 862 cabsl (long double complex z)
 863 {
 864   return hypotl (REALPART (z), IMAGPART (z));
 865 }
 866 #endif
 867
 868
 869 #if !defined(HAVE_CARGF)
 870 #define HAVE_CARGF 1
 871 float cargf (float complex z);
 872
 873 float
 874 cargf (float complex z)
 875 {
 876   return atan2f (IMAGPART (z), REALPART (z));
 877 }
 878 #endif
 879
 880 #if !defined(HAVE_CARG)
 881 #define HAVE_CARG 1
 882 double carg (double complex z);
 883
 884 double
 885 carg (double complex z)
 886 {
 887   return atan2 (IMAGPART (z), REALPART (z));
 888 }
 889 #endif
 890
 891 #if !defined(HAVE_CARGL) && defined(HAVE_ATAN2L)
 892 #define HAVE_CARGL 1
 893 long double cargl (long double complex z);
 894
 895 long double
 896 cargl (long double complex z)
 897 {
 898   return atan2l (IMAGPART (z), REALPART (z));
 899 }
 900 #endif
 901
 902
 903 /* exp(z) = exp(a)*(cos(b) + i sin(b))  */
 904 #if !defined(HAVE_CEXPF)
 905 #define HAVE_CEXPF 1
 906 float complex cexpf (float complex z);
 907
 908 float complex
 909 cexpf (float complex z)
 910 {
 911   float a, b;
 912   float complex v;
 913
 914   a = REALPART (z);
 915   b = IMAGPART (z);
 916   COMPLEX_ASSIGN (v, cosf (b), sinf (b));
 917   return expf (a) * v;
 918 }
 919 #endif
 920
 921 #if !defined(HAVE_CEXP)
 922 #define HAVE_CEXP 1
 923 double complex cexp (double complex z);
 924
 925 double complex
 926 cexp (double complex z)
 927 {
 928   double a, b;
 929   double complex v;
 930
 931   a = REALPART (z);
 932   b = IMAGPART (z);
 933   COMPLEX_ASSIGN (v, cos (b), sin (b));
 934   return exp (a) * v;
 935 }
 936 #endif
 937
 938 #if !defined(HAVE_CEXPL) && defined(HAVE_COSL) && defined(HAVE_SINL) && defined(EXPL)
 939 #define HAVE_CEXPL 1
 940 long double complex cexpl (long double complex z);
 941
 942 long double complex
 943 cexpl (long double complex z)
 944 {
 945   long double a, b;
 946   long double complex v;
 947
 948   a = REALPART (z);
 949   b = IMAGPART (z);
 950   COMPLEX_ASSIGN (v, cosl (b), sinl (b));
 951   return expl (a) * v;
 952 }
 953 #endif
 954
 955
 956 /* log(z) = log (cabs(z)) + i*carg(z)  */
 957 #if !defined(HAVE_CLOGF)
 958 #define HAVE_CLOGF 1
 959 float complex clogf (float complex z);
 960
 961 float complex
 962 clogf (float complex z)
 963 {
 964   float complex v;
 965
 966   COMPLEX_ASSIGN (v, logf (cabsf (z)), cargf (z));
 967   return v;
 968 }
 969 #endif
 970
 971 #if !defined(HAVE_CLOG)
 972 #define HAVE_CLOG 1
 973 double complex clog (double complex z);
 974
 975 double complex
 976 clog (double complex z)
 977 {
 978   double complex v;
 979
 980   COMPLEX_ASSIGN (v, log (cabs (z)), carg (z));
 981   return v;
 982 }
 983 #endif
 984
 985 #if !defined(HAVE_CLOGL) && defined(HAVE_LOGL) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
 986 #define HAVE_CLOGL 1
 987 long double complex clogl (long double complex z);
 988
 989 long double complex
 990 clogl (long double complex z)
 991 {
 992   long double complex v;
 993
 994   COMPLEX_ASSIGN (v, logl (cabsl (z)), cargl (z));
 995   return v;
 996 }
 997 #endif
 998
 999
1000 /* log10(z) = log10 (cabs(z)) + i*carg(z)  */
1001 #if !defined(HAVE_CLOG10F)
1002 #define HAVE_CLOG10F 1
1003 float complex clog10f (float complex z);
1004
1005 float complex
1006 clog10f (float complex z)
1007 {
1008   float complex v;
1009
1010   COMPLEX_ASSIGN (v, log10f (cabsf (z)), cargf (z));
1011   return v;
1012 }
1013 #endif
1014
1015 #if !defined(HAVE_CLOG10)
1016 #define HAVE_CLOG10 1
1017 double complex clog10 (double complex z);
1018
1019 double complex
1020 clog10 (double complex z)
1021 {
1022   double complex v;
1023
1024   COMPLEX_ASSIGN (v, log10 (cabs (z)), carg (z));
1025   return v;
1026 }
1027 #endif
1028
1029 #if !defined(HAVE_CLOG10L) && defined(HAVE_LOG10L) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
1030 #define HAVE_CLOG10L 1
1031 long double complex clog10l (long double complex z);
1032
1033 long double complex
1034 clog10l (long double complex z)
1035 {
1036   long double complex v;
1037
1038   COMPLEX_ASSIGN (v, log10l (cabsl (z)), cargl (z));
1039   return v;
1040 }
1041 #endif
1042
1043
1044 /* pow(base, power) = cexp (power * clog (base))  */
1045 #if !defined(HAVE_CPOWF)
1046 #define HAVE_CPOWF 1
1047 float complex cpowf (float complex base, float complex power);
1048
1049 float complex
1050 cpowf (float complex base, float complex power)
1051 {
1052   return cexpf (power * clogf (base));
1053 }
1054 #endif
1055
1056 #if !defined(HAVE_CPOW)
1057 #define HAVE_CPOW 1
1058 double complex cpow (double complex base, double complex power);
1059
1060 double complex
1061 cpow (double complex base, double complex power)
1062 {
1063   return cexp (power * clog (base));
1064 }
1065 #endif
1066
1067 #if !defined(HAVE_CPOWL) && defined(HAVE_CEXPL) && defined(HAVE_CLOGL)
1068 #define HAVE_CPOWL 1
1069 long double complex cpowl (long double complex base, long double complex power);
1070
1071 long double complex
1072 cpowl (long double complex base, long double complex power)
1073 {
1074   return cexpl (power * clogl (base));
1075 }
1076 #endif
1077
1078
1079 /* sqrt(z).  Algorithm pulled from glibc.  */
1080 #if !defined(HAVE_CSQRTF)
1081 #define HAVE_CSQRTF 1
1082 float complex csqrtf (float complex z);
1083
1084 float complex
1085 csqrtf (float complex z)
1086 {
1087   float re, im;
1088   float complex v;
1089
1090   re = REALPART (z);
1091   im = IMAGPART (z);
1092   if (im == 0)
1093     {
1094       if (re < 0)
1095         {
1096           COMPLEX_ASSIGN (v, 0, copysignf (sqrtf (-re), im));
1097         }
1098       else
1099         {
1100           COMPLEX_ASSIGN (v, fabsf (sqrtf (re)), copysignf (0, im));
1101         }
1102     }
1103   else if (re == 0)
1104     {
1105       float r;
1106
1107       r = sqrtf (0.5 * fabsf (im));
1108
1109       COMPLEX_ASSIGN (v, r, copysignf (r, im));
1110     }
1111   else
1112     {
1113       float d, r, s;
1114
1115       d = hypotf (re, im);
1116       /* Use the identity   2  Re res  Im res = Im x
1117          to avoid cancellation error in  d +/- Re x.  */
1118       if (re > 0)
1119         {
1120           r = sqrtf (0.5 * d + 0.5 * re);
1121           s = (0.5 * im) / r;
1122         }
1123       else
1124         {
1125           s = sqrtf (0.5 * d - 0.5 * re);
1126           r = fabsf ((0.5 * im) / s);
1127         }
1128
1129       COMPLEX_ASSIGN (v, r, copysignf (s, im));
1130     }
1131   return v;
1132 }
1133 #endif
1134
1135 #if !defined(HAVE_CSQRT)
1136 #define HAVE_CSQRT 1
1137 double complex csqrt (double complex z);
1138
1139 double complex
1140 csqrt (double complex z)
1141 {
1142   double re, im;
1143   double complex v;
1144
1145   re = REALPART (z);
1146   im = IMAGPART (z);
1147   if (im == 0)
1148     {
1149       if (re < 0)
1150         {
1151           COMPLEX_ASSIGN (v, 0, copysign (sqrt (-re), im));
1152         }
1153       else
1154         {
1155           COMPLEX_ASSIGN (v, fabs (sqrt (re)), copysign (0, im));
1156         }
1157     }
1158   else if (re == 0)
1159     {
1160       double r;
1161
1162       r = sqrt (0.5 * fabs (im));
1163
1164       COMPLEX_ASSIGN (v, r, copysign (r, im));
1165     }
1166   else
1167     {
1168       double d, r, s;
1169
1170       d = hypot (re, im);
1171       /* Use the identity   2  Re res  Im res = Im x
1172          to avoid cancellation error in  d +/- Re x.  */
1173       if (re > 0)
1174         {
1175           r = sqrt (0.5 * d + 0.5 * re);
1176           s = (0.5 * im) / r;
1177         }
1178       else
1179         {
1180           s = sqrt (0.5 * d - 0.5 * re);
1181           r = fabs ((0.5 * im) / s);
1182         }
1183
1184       COMPLEX_ASSIGN (v, r, copysign (s, im));
1185     }
1186   return v;
1187 }
1188 #endif
1189
1190 #if !defined(HAVE_CSQRTL) && defined(HAVE_COPYSIGNL) && defined(HAVE_SQRTL) && defined(HAVE_FABSL) && defined(HAVE_HYPOTL)
1191 #define HAVE_CSQRTL 1
1192 long double complex csqrtl (long double complex z);
1193
1194 long double complex
1195 csqrtl (long double complex z)
1196 {
1197   long double re, im;
1198   long double complex v;
1199
1200   re = REALPART (z);
1201   im = IMAGPART (z);
1202   if (im == 0)
1203     {
1204       if (re < 0)
1205         {
1206           COMPLEX_ASSIGN (v, 0, copysignl (sqrtl (-re), im));
1207         }
1208       else
1209         {
1210           COMPLEX_ASSIGN (v, fabsl (sqrtl (re)), copysignl (0, im));
1211         }
1212     }
1213   else if (re == 0)
1214     {
1215       long double r;
1216
1217       r = sqrtl (0.5 * fabsl (im));
1218
1219       COMPLEX_ASSIGN (v, copysignl (r, im), r);
1220     }
1221   else
1222     {
1223       long double d, r, s;
1224
1225       d = hypotl (re, im);
1226       /* Use the identity   2  Re res  Im res = Im x
1227          to avoid cancellation error in  d +/- Re x.  */
1228       if (re > 0)
1229         {
1230           r = sqrtl (0.5 * d + 0.5 * re);
1231           s = (0.5 * im) / r;
1232         }
1233       else
1234         {
1235           s = sqrtl (0.5 * d - 0.5 * re);
1236           r = fabsl ((0.5 * im) / s);
1237         }
1238
1239       COMPLEX_ASSIGN (v, r, copysignl (s, im));
1240     }
1241   return v;
1242 }
1243 #endif
1244
1245
1246 /* sinh(a + i b) = sinh(a) cos(b) + i cosh(a) sin(b)  */
1247 #if !defined(HAVE_CSINHF)
1248 #define HAVE_CSINHF 1
1249 float complex csinhf (float complex a);
1250
1251 float complex
1252 csinhf (float complex a)
1253 {
1254   float r, i;
1255   float complex v;
1256
1257   r = REALPART (a);
1258   i = IMAGPART (a);
1259   COMPLEX_ASSIGN (v, sinhf (r) * cosf (i), coshf (r) * sinf (i));
1260   return v;
1261 }
1262 #endif
1263
1264 #if !defined(HAVE_CSINH)
1265 #define HAVE_CSINH 1
1266 double complex csinh (double complex a);
1267
1268 double complex
1269 csinh (double complex a)
1270 {
1271   double r, i;
1272   double complex v;
1273
1274   r = REALPART (a);
1275   i = IMAGPART (a);
1276   COMPLEX_ASSIGN (v, sinh (r) * cos (i), cosh (r) * sin (i));
1277   return v;
1278 }
1279 #endif
1280
1281 #if !defined(HAVE_CSINHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1282 #define HAVE_CSINHL 1
1283 long double complex csinhl (long double complex a);
1284
1285 long double complex
1286 csinhl (long double complex a)
1287 {
1288   long double r, i;
1289   long double complex v;
1290
1291   r = REALPART (a);
1292   i = IMAGPART (a);
1293   COMPLEX_ASSIGN (v, sinhl (r) * cosl (i), coshl (r) * sinl (i));
1294   return v;
1295 }
1296 #endif
1297
1298
1299 /* cosh(a + i b) = cosh(a) cos(b) + i sinh(a) sin(b)  */
1300 #if !defined(HAVE_CCOSHF)
1301 #define HAVE_CCOSHF 1
1302 float complex ccoshf (float complex a);
1303
1304 float complex
1305 ccoshf (float complex a)
1306 {
1307   float r, i;
1308   float complex v;
1309
1310   r = REALPART (a);
1311   i = IMAGPART (a);
1312   COMPLEX_ASSIGN (v, coshf (r) * cosf (i), sinhf (r) * sinf (i));
1313   return v;
1314 }
1315 #endif
1316
1317 #if !defined(HAVE_CCOSH)
1318 #define HAVE_CCOSH 1
1319 double complex ccosh (double complex a);
1320
1321 double complex
1322 ccosh (double complex a)
1323 {
1324   double r, i;
1325   double complex v;
1326
1327   r = REALPART (a);
1328   i = IMAGPART (a);
1329   COMPLEX_ASSIGN (v, cosh (r) * cos (i),  sinh (r) * sin (i));
1330   return v;
1331 }
1332 #endif
1333
1334 #if !defined(HAVE_CCOSHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1335 #define HAVE_CCOSHL 1
1336 long double complex ccoshl (long double complex a);
1337
1338 long double complex
1339 ccoshl (long double complex a)
1340 {
1341   long double r, i;
1342   long double complex v;
1343
1344   r = REALPART (a);
1345   i = IMAGPART (a);
1346   COMPLEX_ASSIGN (v, coshl (r) * cosl (i), sinhl (r) * sinl (i));
1347   return v;
1348 }
1349 #endif
1350
1351
1352 /* tanh(a + i b) = (tanh(a) + i tan(b)) / (1 + i tanh(a) tan(b))  */
1353 #if !defined(HAVE_CTANHF)
1354 #define HAVE_CTANHF 1
1355 float complex ctanhf (float complex a);
1356
1357 float complex
1358 ctanhf (float complex a)
1359 {
1360   float rt, it;
1361   float complex n, d;
1362
1363   rt = tanhf (REALPART (a));
1364   it = tanf (IMAGPART (a));
1365   COMPLEX_ASSIGN (n, rt, it);
1366   COMPLEX_ASSIGN (d, 1, rt * it);
1367
1368   return n / d;
1369 }
1370 #endif
1371
1372 #if !defined(HAVE_CTANH)
1373 #define HAVE_CTANH 1
1374 double complex ctanh (double complex a);
1375 double complex
1376 ctanh (double complex a)
1377 {
1378   double rt, it;
1379   double complex n, d;
1380
1381   rt = tanh (REALPART (a));
1382   it = tan (IMAGPART (a));
1383   COMPLEX_ASSIGN (n, rt, it);
1384   COMPLEX_ASSIGN (d, 1, rt * it);
1385
1386   return n / d;
1387 }
1388 #endif
1389
1390 #if !defined(HAVE_CTANHL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1391 #define HAVE_CTANHL 1
1392 long double complex ctanhl (long double complex a);
1393
1394 long double complex
1395 ctanhl (long double complex a)
1396 {
1397   long double rt, it;
1398   long double complex n, d;
1399
1400   rt = tanhl (REALPART (a));
1401   it = tanl (IMAGPART (a));
1402   COMPLEX_ASSIGN (n, rt, it);
1403   COMPLEX_ASSIGN (d, 1, rt * it);
1404
1405   return n / d;
1406 }
1407 #endif
1408
1409
1410 /* sin(a + i b) = sin(a) cosh(b) + i cos(a) sinh(b)  */
1411 #if !defined(HAVE_CSINF)
1412 #define HAVE_CSINF 1
1413 float complex csinf (float complex a);
1414
1415 float complex
1416 csinf (float complex a)
1417 {
1418   float r, i;
1419   float complex v;
1420
1421   r = REALPART (a);
1422   i = IMAGPART (a);
1423   COMPLEX_ASSIGN (v, sinf (r) * coshf (i), cosf (r) * sinhf (i));
1424   return v;
1425 }
1426 #endif
1427
1428 #if !defined(HAVE_CSIN)
1429 #define HAVE_CSIN 1
1430 double complex csin (double complex a);
1431
1432 double complex
1433 csin (double complex a)
1434 {
1435   double r, i;
1436   double complex v;
1437
1438   r = REALPART (a);
1439   i = IMAGPART (a);
1440   COMPLEX_ASSIGN (v, sin (r) * cosh (i), cos (r) * sinh (i));
1441   return v;
1442 }
1443 #endif
1444
1445 #if !defined(HAVE_CSINL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1446 #define HAVE_CSINL 1
1447 long double complex csinl (long double complex a);
1448
1449 long double complex
1450 csinl (long double complex a)
1451 {
1452   long double r, i;
1453   long double complex v;
1454
1455   r = REALPART (a);
1456   i = IMAGPART (a);
1457   COMPLEX_ASSIGN (v, sinl (r) * coshl (i), cosl (r) * sinhl (i));
1458   return v;
1459 }
1460 #endif
1461
1462
1463 /* cos(a + i b) = cos(a) cosh(b) - i sin(a) sinh(b)  */
1464 #if !defined(HAVE_CCOSF)
1465 #define HAVE_CCOSF 1
1466 float complex ccosf (float complex a);
1467
1468 float complex
1469 ccosf (float complex a)
1470 {
1471   float r, i;
1472   float complex v;
1473
1474   r = REALPART (a);
1475   i = IMAGPART (a);
1476   COMPLEX_ASSIGN (v, cosf (r) * coshf (i), - (sinf (r) * sinhf (i)));
1477   return v;
1478 }
1479 #endif
1480
1481 #if !defined(HAVE_CCOS)
1482 #define HAVE_CCOS 1
1483 double complex ccos (double complex a);
1484
1485 double complex
1486 ccos (double complex a)
1487 {
1488   double r, i;
1489   double complex v;
1490
1491   r = REALPART (a);
1492   i = IMAGPART (a);
1493   COMPLEX_ASSIGN (v, cos (r) * cosh (i), - (sin (r) * sinh (i)));
1494   return v;
1495 }
1496 #endif
1497
1498 #if !defined(HAVE_CCOSL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1499 #define HAVE_CCOSL 1
1500 long double complex ccosl (long double complex a);
1501
1502 long double complex
1503 ccosl (long double complex a)
1504 {
1505   long double r, i;
1506   long double complex v;
1507
1508   r = REALPART (a);
1509   i = IMAGPART (a);
1510   COMPLEX_ASSIGN (v, cosl (r) * coshl (i), - (sinl (r) * sinhl (i)));
1511   return v;
1512 }
1513 #endif
1514
1515
1516 /* tan(a + i b) = (tan(a) + i tanh(b)) / (1 - i tan(a) tanh(b))  */
1517 #if !defined(HAVE_CTANF)
1518 #define HAVE_CTANF 1
1519 float complex ctanf (float complex a);
1520
1521 float complex
1522 ctanf (float complex a)
1523 {
1524   float rt, it;
1525   float complex n, d;
1526
1527   rt = tanf (REALPART (a));
1528   it = tanhf (IMAGPART (a));
1529   COMPLEX_ASSIGN (n, rt, it);
1530   COMPLEX_ASSIGN (d, 1, - (rt * it));
1531
1532   return n / d;
1533 }
1534 #endif
1535
1536 #if !defined(HAVE_CTAN)
1537 #define HAVE_CTAN 1
1538 double complex ctan (double complex a);
1539
1540 double complex
1541 ctan (double complex a)
1542 {
1543   double rt, it;
1544   double complex n, d;
1545
1546   rt = tan (REALPART (a));
1547   it = tanh (IMAGPART (a));
1548   COMPLEX_ASSIGN (n, rt, it);
1549   COMPLEX_ASSIGN (d, 1, - (rt * it));
1550
1551   return n / d;
1552 }
1553 #endif
1554
1555 #if !defined(HAVE_CTANL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1556 #define HAVE_CTANL 1
1557 long double complex ctanl (long double complex a);
1558
1559 long double complex
1560 ctanl (long double complex a)
1561 {
1562   long double rt, it;
1563   long double complex n, d;
1564
1565   rt = tanl (REALPART (a));
1566   it = tanhl (IMAGPART (a));
1567   COMPLEX_ASSIGN (n, rt, it);
1568   COMPLEX_ASSIGN (d, 1, - (rt * it));
1569
1570   return n / d;
1571 }
1572 #endif
1573
1574
1575 /* Complex ASIN.  Returns wrongly NaN for infinite arguments.
1576    Algorithm taken from Abramowitz & Stegun.  */
1577
1578 #if !defined(HAVE_CASINF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1579 #define HAVE_CASINF 1
1580 complex float casinf (complex float z);
1581
1582 complex float
1583 casinf (complex float z)
1584 {
1585   return -I*clogf (I*z + csqrtf (1.0f-z*z));
1586 }
1587 #endif
1588
1589
1590 #if !defined(HAVE_CASIN) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1591 #define HAVE_CASIN 1
1592 complex double casin (complex double z);
1593
1594 complex double
1595 casin (complex double z)
1596 {
1597   return -I*clog (I*z + csqrt (1.0-z*z));
1598 }
1599 #endif
1600
1601
1602 #if !defined(HAVE_CASINL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1603 #define HAVE_CASINL 1
1604 complex long double casinl (complex long double z);
1605
1606 complex long double
1607 casinl (complex long double z)
1608 {
1609   return -I*clogl (I*z + csqrtl (1.0L-z*z));
1610 }
1611 #endif
1612
1613
1614 /* Complex ACOS.  Returns wrongly NaN for infinite arguments.
1615    Algorithm taken from Abramowitz & Stegun.  */
1616
1617 #if !defined(HAVE_CACOSF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1618 #define HAVE_CACOSF 1
1619 complex float cacosf (complex float z);
1620
1621 complex float
1622 cacosf (complex float z)
1623 {
1624   return -I*clogf (z + I*csqrtf (1.0f-z*z));
1625 }
1626 #endif
1627
1628
1629 #if !defined(HAVE_CACOS) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1630 #define HAVE_CACOS 1
1631 complex double cacos (complex double z);
1632
1633 complex double
1634 cacos (complex double z)
1635 {
1636   return -I*clog (z + I*csqrt (1.0-z*z));
1637 }
1638 #endif
1639
1640
1641 #if !defined(HAVE_CACOSL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1642 #define HAVE_CACOSL 1
1643 complex long double cacosl (complex long double z);
1644
1645 complex long double
1646 cacosl (complex long double z)
1647 {
1648   return -I*clogl (z + I*csqrtl (1.0L-z*z));
1649 }
1650 #endif
1651
1652
1653 /* Complex ATAN.  Returns wrongly NaN for infinite arguments.
1654    Algorithm taken from Abramowitz & Stegun.  */
1655
1656 #if !defined(HAVE_CATANF) && defined(HAVE_CLOGF)
1657 #define HAVE_CACOSF 1
1658 complex float catanf (complex float z);
1659
1660 complex float
1661 catanf (complex float z)
1662 {
1663   return I*clogf ((I+z)/(I-z))/2.0f;
1664 }
1665 #endif
1666
1667
1668 #if !defined(HAVE_CATAN) && defined(HAVE_CLOG)
1669 #define HAVE_CACOS 1
1670 complex double catan (complex double z);
1671
1672 complex double
1673 catan (complex double z)
1674 {
1675   return I*clog ((I+z)/(I-z))/2.0;
1676 }
1677 #endif
1678
1679
1680 #if !defined(HAVE_CATANL) && defined(HAVE_CLOGL)
1681 #define HAVE_CACOSL 1
1682 complex long double catanl (complex long double z);
1683
1684 complex long double
1685 catanl (complex long double z)
1686 {
1687   return I*clogl ((I+z)/(I-z))/2.0L;
1688 }
1689 #endif
1690
1691
1692 /* Complex ASINH.  Returns wrongly NaN for infinite arguments.
1693    Algorithm taken from Abramowitz & Stegun.  */
1694
1695 #if !defined(HAVE_CASINHF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1696 #define HAVE_CASINHF 1
1697 complex float casinhf (complex float z);
1698
1699 complex float
1700 casinhf (complex float z)
1701 {
1702   return clogf (z + csqrtf (z*z+1.0f));
1703 }
1704 #endif
1705
1706
1707 #if !defined(HAVE_CASINH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1708 #define HAVE_CASINH 1
1709 complex double casinh (complex double z);
1710
1711 complex double
1712 casinh (complex double z)
1713 {
1714   return clog (z + csqrt (z*z+1.0));
1715 }
1716 #endif
1717
1718
1719 #if !defined(HAVE_CASINHL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1720 #define HAVE_CASINHL 1
1721 complex long double casinhl (complex long double z);
1722
1723 complex long double
1724 casinhl (complex long double z)
1725 {
1726   return clogl (z + csqrtl (z*z+1.0L));
1727 }
1728 #endif
1729
1730
1731 /* Complex ACOSH.  Returns wrongly NaN for infinite arguments.
1732    Algorithm taken from Abramowitz & Stegun.  */
1733
1734 #if !defined(HAVE_CACOSHF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1735 #define HAVE_CACOSHF 1
1736 complex float cacoshf (complex float z);
1737
1738 complex float
1739 cacoshf (complex float z)
1740 {
1741   return clogf (z + csqrtf (z-1.0f) * csqrtf (z+1.0f));
1742 }
1743 #endif
1744
1745
1746 #if !defined(HAVE_CACOSH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1747 #define HAVE_CACOSH 1
1748 complex double cacosh (complex double z);
1749
1750 complex double
1751 cacosh (complex double z)
1752 {
1753   return clog (z + csqrt (z-1.0) * csqrt (z+1.0));
1754 }
1755 #endif
1756
1757
1758 #if !defined(HAVE_CACOSHL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1759 #define HAVE_CACOSHL 1
1760 complex long double cacoshl (complex long double z);
1761
1762 complex long double
1763 cacoshl (complex long double z)
1764 {
1765   return clogl (z + csqrtl (z-1.0L) * csqrtl (z+1.0L));
1766 }
1767 #endif
1768
1769
1770 /* Complex ATANH.  Returns wrongly NaN for infinite arguments.
1771    Algorithm taken from Abramowitz & Stegun.  */
1772
1773 #if !defined(HAVE_CATANHF) && defined(HAVE_CLOGF)
1774 #define HAVE_CATANHF 1
1775 complex float catanhf (complex float z);
1776
1777 complex float
1778 catanhf (complex float z)
1779 {
1780   return clogf ((1.0f+z)/(1.0f-z))/2.0f;
1781 }
1782 #endif
1783
1784
1785 #if !defined(HAVE_CATANH) && defined(HAVE_CLOG)
1786 #define HAVE_CATANH 1
1787 complex double catanh (complex double z);
1788
1789 complex double
1790 catanh (complex double z)
1791 {
1792   return clog ((1.0+z)/(1.0-z))/2.0;
1793 }
1794 #endif
1795
1796 #if !defined(HAVE_CATANHL) && defined(HAVE_CLOGL)
1797 #define HAVE_CATANHL 1
1798 complex long double catanhl (complex long double z);
1799
1800 complex long double
1801 catanhl (complex long double z)
1802 {
1803   return clogl ((1.0L+z)/(1.0L-z))/2.0L;
1804 }
1805 #endif
1806
1807
1808 #if !defined(HAVE_TGAMMA)
1809 #define HAVE_TGAMMA 1
1810 double tgamma (double);
1811
1812 /* Fallback tgamma() function. Uses the algorithm from
1813    http://www.netlib.org/specfun/gamma and references therein.  */
1814
1815 #undef SQRTPI
1816 #define SQRTPI 0.9189385332046727417803297
1817
1818 #undef PI
1819 #define PI 3.1415926535897932384626434
1820
1821 double
1822 tgamma (double x)
1823 {
1824   int i, n, parity;
1825   double fact, res, sum, xden, xnum, y, y1, ysq, z;
1826
1827   static double p[8] = {
1828     -1.71618513886549492533811e0,  2.47656508055759199108314e1,
1829     -3.79804256470945635097577e2,  6.29331155312818442661052e2,
1830      8.66966202790413211295064e2, -3.14512729688483675254357e4,
1831     -3.61444134186911729807069e4,  6.64561438202405440627855e4 };
1832
1833   static double q[8] = {
1834     -3.08402300119738975254353e1,  3.15350626979604161529144e2,
1835     -1.01515636749021914166146e3, -3.10777167157231109440444e3,
1836      2.25381184209801510330112e4,  4.75584627752788110767815e3,
1837     -1.34659959864969306392456e5, -1.15132259675553483497211e5 };
1838
1839   static double c[7] = {             -1.910444077728e-03,
1840      8.4171387781295e-04,            -5.952379913043012e-04,
1841      7.93650793500350248e-04,        -2.777777777777681622553e-03,
1842      8.333333333333333331554247e-02,  5.7083835261e-03 };
1843
1844   static const double xminin = 2.23e-308;
1845   static const double xbig = 171.624;
1846   static const double xnan = __builtin_nan ("0x0"), xinf = __builtin_inf ();
1847   static double eps = 0;
1848
1849   if (eps == 0)
1850     eps = nextafter (1., 2.) - 1.;
1851
1852   parity = 0;
1853   fact = 1;
1854   n = 0;
1855   y = x;
1856
1857   if (__builtin_isnan (x))
1858     return x;
1859
1860   if (y <= 0)
1861     {
1862       y = -x;
1863       y1 = trunc (y);
1864       res = y - y1;
1865
1866       if (res != 0)
1867         {
1868           if (y1 != trunc (y1*0.5l)*2)
1869             parity = 1;
1870           fact = -PI / sin (PI*res);
1871           y = y + 1;
1872         }
1873       else
1874         return x == 0 ? copysign (xinf, x) : xnan;
1875     }
1876
1877   if (y < eps)
1878     {
1879       if (y >= xminin)
1880         res = 1 / y;
1881       else
1882         return xinf;
1883     }
1884   else if (y < 13)
1885     {
1886       y1 = y;
1887       if (y < 1)
1888         {
1889           z = y;
1890           y = y + 1;
1891         }
1892       else
1893         {
1894           n = (int)y - 1;
1895           y = y - n;
1896           z = y - 1;
1897         }
1898
1899       xnum = 0;
1900       xden = 1;
1901       for (i = 0; i < 8; i++)
1902         {
1903           xnum = (xnum + p[i]) * z;
1904           xden = xden * z + q[i];
1905         }
1906
1907       res = xnum / xden + 1;
1908
1909       if (y1 < y)
1910         res = res / y1;
1911       else if (y1 > y)
1912         for (i = 1; i <= n; i++)
1913           {
1914             res = res * y;
1915             y = y + 1;
1916           }
1917     }
1918   else
1919     {
1920       if (y < xbig)
1921         {
1922           ysq = y * y;
1923           sum = c[6];
1924           for (i = 0; i < 6; i++)
1925             sum = sum / ysq + c[i];
1926
1927           sum = sum/y - y + SQRTPI;
1928           sum = sum + (y - 0.5) * log (y);
1929           res = exp (sum);
1930         }
1931       else
1932         return x < 0 ? xnan : xinf;
1933     }
1934
1935   if (parity)
1936     res = -res;
1937   if (fact != 1)
1938     res = fact / res;
1939
1940   return res;
1941 }
1942 #endif
1943
1944
1945
1946 #if !defined(HAVE_LGAMMA)
1947 #define HAVE_LGAMMA 1
1948 double lgamma (double);
1949
1950 /* Fallback lgamma() function. Uses the algorithm from
1951    http://www.netlib.org/specfun/algama and references therein,
1952    except for negative arguments (where netlib would return +Inf)
1953    where we use the following identity:
1954        lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y)
1955  */
1956
1957 double
1958 lgamma (double y)
1959 {
1960
1961 #undef SQRTPI
1962 #define SQRTPI 0.9189385332046727417803297
1963
1964 #undef PI
1965 #define PI 3.1415926535897932384626434
1966
1967 #define PNT68  0.6796875
1968 #define D1    -0.5772156649015328605195174
1969 #define D2     0.4227843350984671393993777
1970 #define D4     1.791759469228055000094023
1971
1972   static double p1[8] = {
1973               4.945235359296727046734888e0, 2.018112620856775083915565e2,
1974               2.290838373831346393026739e3, 1.131967205903380828685045e4,
1975               2.855724635671635335736389e4, 3.848496228443793359990269e4,
1976               2.637748787624195437963534e4, 7.225813979700288197698961e3 };
1977   static double q1[8] = {
1978               6.748212550303777196073036e1,  1.113332393857199323513008e3,
1979               7.738757056935398733233834e3,  2.763987074403340708898585e4,
1980               5.499310206226157329794414e4,  6.161122180066002127833352e4,
1981               3.635127591501940507276287e4,  8.785536302431013170870835e3 };
1982   static double p2[8] = {
1983               4.974607845568932035012064e0,  5.424138599891070494101986e2,
1984               1.550693864978364947665077e4,  1.847932904445632425417223e5,
1985               1.088204769468828767498470e6,  3.338152967987029735917223e6,
1986               5.106661678927352456275255e6,  3.074109054850539556250927e6 };
1987   static double q2[8] = {
1988               1.830328399370592604055942e2,  7.765049321445005871323047e3,
1989               1.331903827966074194402448e5,  1.136705821321969608938755e6,
1990               5.267964117437946917577538e6,  1.346701454311101692290052e7,
1991               1.782736530353274213975932e7,  9.533095591844353613395747e6 };
1992   static double p4[8] = {
1993               1.474502166059939948905062e4,  2.426813369486704502836312e6,
1994               1.214755574045093227939592e8,  2.663432449630976949898078e9,
1995               2.940378956634553899906876e10, 1.702665737765398868392998e11,
1996               4.926125793377430887588120e11, 5.606251856223951465078242e11 };
1997   static double q4[8] = {
1998               2.690530175870899333379843e3,  6.393885654300092398984238e5,
1999               4.135599930241388052042842e7,  1.120872109616147941376570e9,
2000               1.488613728678813811542398e10, 1.016803586272438228077304e11,
2001               3.417476345507377132798597e11, 4.463158187419713286462081e11 };
2002   static double  c[7] = {
2003              -1.910444077728e-03,            8.4171387781295e-04,
2004              -5.952379913043012e-04,         7.93650793500350248e-04,
2005              -2.777777777777681622553e-03,   8.333333333333333331554247e-02,
2006               5.7083835261e-03 };
2007
2008   static double xbig = 2.55e305, xinf = __builtin_inf (), eps = 0,
2009                 frtbig = 2.25e76;
2010
2011   int i;
2012   double corr, res, xden, xm1, xm2, xm4, xnum, ysq;
2013
2014   if (eps == 0)
2015     eps = __builtin_nextafter (1., 2.) - 1.;
2016
2017   if ((y > 0) && (y <= xbig))
2018     {
2019       if (y <= eps)
2020         res = -log (y);
2021       else if (y <= 1.5)
2022         {
2023           if (y < PNT68)
2024             {
2025               corr = -log (y);
2026               xm1 = y;
2027             }
2028           else
2029             {
2030               corr = 0;
2031               xm1 = (y - 0.5) - 0.5;
2032             }
2033
2034           if ((y <= 0.5) || (y >= PNT68))
2035             {
2036               xden = 1;
2037               xnum = 0;
2038               for (i = 0; i < 8; i++)
2039                 {
2040                   xnum = xnum*xm1 + p1[i];
2041                   xden = xden*xm1 + q1[i];
2042                 }
2043               res = corr + (xm1 * (D1 + xm1*(xnum/xden)));
2044             }
2045           else
2046             {
2047               xm2 = (y - 0.5) - 0.5;
2048               xden = 1;
2049               xnum = 0;
2050               for (i = 0; i < 8; i++)
2051                 {
2052                   xnum = xnum*xm2 + p2[i];
2053                   xden = xden*xm2 + q2[i];
2054                 }
2055               res = corr + xm2 * (D2 + xm2*(xnum/xden));
2056             }
2057         }
2058       else if (y <= 4)
2059         {
2060           xm2 = y - 2;
2061           xden = 1;
2062           xnum = 0;
2063           for (i = 0; i < 8; i++)
2064             {
2065               xnum = xnum*xm2 + p2[i];
2066               xden = xden*xm2 + q2[i];
2067             }
2068           res = xm2 * (D2 + xm2*(xnum/xden));
2069         }
2070       else if (y <= 12)
2071         {
2072           xm4 = y - 4;
2073           xden = -1;
2074           xnum = 0;
2075           for (i = 0; i < 8; i++)
2076             {
2077               xnum = xnum*xm4 + p4[i];
2078               xden = xden*xm4 + q4[i];
2079             }
2080           res = D4 + xm4*(xnum/xden);
2081         }
2082       else
2083         {
2084           res = 0;
2085           if (y <= frtbig)
2086             {
2087               res = c[6];
2088               ysq = y * y;
2089               for (i = 0; i < 6; i++)
2090                 res = res / ysq + c[i];
2091             }
2092           res = res/y;
2093           corr = log (y);
2094           res = res + SQRTPI - 0.5*corr;
2095           res = res + y*(corr-1);
2096         }
2097     }
2098   else if (y < 0 && __builtin_floor (y) != y)
2099     {
2100       /* lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y)
2101          For abs(y) very close to zero, we use a series expansion to
2102          the first order in y to avoid overflow.  */
2103       if (y > -1.e-100)
2104         res = -2 * log (fabs (y)) - lgamma (-y);
2105       else
2106         res = log (PI / fabs (y * sin (PI * y))) - lgamma (-y);
2107     }
2108   else
2109     res = xinf;
2110
2111   return res;
2112 }
2113 #endif
2114
2115
2116 #if defined(HAVE_TGAMMA) && !defined(HAVE_TGAMMAF)
2117 #define HAVE_TGAMMAF 1
2118 float tgammaf (float);
2119
2120 float
2121 tgammaf (float x)
2122 {
2123   return (float) tgamma ((double) x);
2124 }
2125 #endif
2126
2127 #if defined(HAVE_LGAMMA) && !defined(HAVE_LGAMMAF)
2128 #define HAVE_LGAMMAF 1
2129 float lgammaf (float);
2130
2131 float
2132 lgammaf (float x)
2133 {
2134   return (float) lgamma ((double) x);
2135 }
2136 #endif