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