1 /* Implementation of various C99 functions 2 Copyright (C) 2004-2020 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 2192