libgfortran/intrinsics/c99_functions.c

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