1/* Implementation of various C99 functions 2 Copyright (C) 2004-2020 Free Software Foundation, Inc. 3 4This file is part of the GNU Fortran 95 runtime library (libgfortran). 5 6Libgfortran is free software; you can redistribute it and/or 7modify it under the terms of the GNU General Public 8License as published by the Free Software Foundation; either 9version 3 of the License, or (at your option) any later version. 10 11Libgfortran is distributed in the hope that it will be useful, 12but WITHOUT ANY WARRANTY; without even the implied warranty of 13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14GNU General Public License for more details. 15 16Under Section 7 of GPL version 3, you are granted additional 17permissions described in the GCC Runtime Library Exception, version 183.1, as published by the Free Software Foundation. 19 20You should have received a copy of the GNU General Public License and 21a copy of the GCC Runtime Library Exception along with this program; 22see 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 58float j0f (float); 59 60float 61j0f (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 69float j1f (float); 70 71float 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 79float jnf (int, float); 80 81float 82jnf (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 90float y0f (float); 91 92float 93y0f (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 101float y1f (float); 102 103float 104y1f (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 112float ynf (int, float); 113 114float 115ynf (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 126float erff (float); 127 128float 129erff (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 137float erfcf (float); 138 139float 140erfcf (float x) 141{ 142 return (float) erfc ((double) x); 143} 144#endif 145 146 147#ifndef HAVE_ACOSF 148#define HAVE_ACOSF 1 149float acosf (float x); 150 151float 152acosf (float x) 153{ 154 return (float) acos (x); 155} 156#endif 157 158#if HAVE_ACOSH && !HAVE_ACOSHF 159float acoshf (float x); 160 161float 162acoshf (float x) 163{ 164 return (float) acosh ((double) x); 165} 166#endif 167 168#ifndef HAVE_ASINF 169#define HAVE_ASINF 1 170float asinf (float x); 171 172float 173asinf (float x) 174{ 175 return (float) asin (x); 176} 177#endif 178 179#if HAVE_ASINH && !HAVE_ASINHF 180float asinhf (float x); 181 182float 183asinhf (float x) 184{ 185 return (float) asinh ((double) x); 186} 187#endif 188 189#ifndef HAVE_ATAN2F 190#define HAVE_ATAN2F 1 191float atan2f (float y, float x); 192 193float 194atan2f (float y, float x) 195{ 196 return (float) atan2 (y, x); 197} 198#endif 199 200#ifndef HAVE_ATANF 201#define HAVE_ATANF 1 202float atanf (float x); 203 204float 205atanf (float x) 206{ 207 return (float) atan (x); 208} 209#endif 210 211#if HAVE_ATANH && !HAVE_ATANHF 212float atanhf (float x); 213 214float 215atanhf (float x) 216{ 217 return (float) atanh ((double) x); 218} 219#endif 220 221#ifndef HAVE_CEILF 222#define HAVE_CEILF 1 223float ceilf (float x); 224 225float 226ceilf (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 234double copysign (double x, double y); 235 236double 237copysign (double x, double y) 238{ 239 return __builtin_copysign (x, y); 240} 241#endif 242 243#ifndef HAVE_COPYSIGNF 244#define HAVE_COPYSIGNF 1 245float copysignf (float x, float y); 246 247float 248copysignf (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 256long double copysignl (long double x, long double y); 257 258long double 259copysignl (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 267float cosf (float x); 268 269float 270cosf (float x) 271{ 272 return (float) cos (x); 273} 274#endif 275 276#ifndef HAVE_COSHF 277#define HAVE_COSHF 1 278float coshf (float x); 279 280float 281coshf (float x) 282{ 283 return (float) cosh (x); 284} 285#endif 286 287#ifndef HAVE_EXPF 288#define HAVE_EXPF 1 289float expf (float x); 290 291float 292expf (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 300double fabs (double x); 301 302double 303fabs (double x) 304{ 305 return __builtin_fabs (x); 306} 307#endif 308 309#ifndef HAVE_FABSF 310#define HAVE_FABSF 1 311float fabsf (float x); 312 313float 314fabsf (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 322long double fabsl (long double x); 323 324long double 325fabsl (long double x) 326{ 327 return __builtin_fabsl (x); 328} 329#endif 330 331#ifndef HAVE_FLOORF 332#define HAVE_FLOORF 1 333float floorf (float x); 334 335float 336floorf (float x) 337{ 338 return (float) floor (x); 339} 340#endif 341 342#ifndef HAVE_FMODF 343#define HAVE_FMODF 1 344float fmodf (float x, float y); 345 346float 347fmodf (float x, float y) 348{ 349 return (float) fmod (x, y); 350} 351#endif 352 353#ifndef HAVE_FREXPF 354#define HAVE_FREXPF 1 355float frexpf (float x, int *exp); 356 357float 358frexpf (float x, int *exp) 359{ 360 return (float) frexp (x, exp); 361} 362#endif 363 364#ifndef HAVE_HYPOTF 365#define HAVE_HYPOTF 1 366float hypotf (float x, float y); 367 368float 369hypotf (float x, float y) 370{ 371 return (float) hypot (x, y); 372} 373#endif 374 375#ifndef HAVE_LOGF 376#define HAVE_LOGF 1 377float logf (float x); 378 379float 380logf (float x) 381{ 382 return (float) log (x); 383} 384#endif 385 386#ifndef HAVE_LOG10F 387#define HAVE_LOG10F 1 388float log10f (float x); 389 390float 391log10f (float x) 392{ 393 return (float) log10 (x); 394} 395#endif 396 397#ifndef HAVE_SCALBN 398#define HAVE_SCALBN 1 399double scalbn (double x, int y); 400 401double 402scalbn (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 414float scalbnf (float x, int y); 415 416float 417scalbnf (float x, int y) 418{ 419 return (float) scalbn (x, y); 420} 421#endif 422 423#ifndef HAVE_SINF 424#define HAVE_SINF 1 425float sinf (float x); 426 427float 428sinf (float x) 429{ 430 return (float) sin (x); 431} 432#endif 433 434#ifndef HAVE_SINHF 435#define HAVE_SINHF 1 436float sinhf (float x); 437 438float 439sinhf (float x) 440{ 441 return (float) sinh (x); 442} 443#endif 444 445#ifndef HAVE_SQRTF 446#define HAVE_SQRTF 1 447float sqrtf (float x); 448 449float 450sqrtf (float x) 451{ 452 return (float) sqrt (x); 453} 454#endif 455 456#ifndef HAVE_TANF 457#define HAVE_TANF 1 458float tanf (float x); 459 460float 461tanf (float x) 462{ 463 return (float) tan (x); 464} 465#endif 466 467#ifndef HAVE_TANHF 468#define HAVE_TANHF 1 469float tanhf (float x); 470 471float 472tanhf (float x) 473{ 474 return (float) tanh (x); 475} 476#endif 477 478#ifndef HAVE_TRUNC 479#define HAVE_TRUNC 1 480double trunc (double x); 481 482double 483trunc (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 497float truncf (float x); 498 499float 500truncf (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. */ 511float nextafterf (float x, float y); 512 513float 514nextafterf (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 574float powf (float x, float y); 575 576float 577powf (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. */ 588double round (double x); 589 590double 591round (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 619long 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 625long double 626roundl (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. */ 650long double 651roundl (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. */ 674float roundf (float x); 675 676float 677roundf (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 705long int lroundf (float x); 706 707long int 708lroundf (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 716long int lround (double x); 717 718long int 719lround (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 727long int lroundl (long double x); 728 729long int 730lroundl (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 738long long int llroundf (float x); 739 740long long int 741llroundf (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 749long long int llround (double x); 750 751long long int 752llround (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 760long long int llroundl (long double x); 761 762long long int 763llroundl (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. */ 774long double log10l (long double x); 775 776long double 777log10l (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 814long double floorl (long double x); 815 816long double 817floorl (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 842long double fmodl (long double x, long double y); 843 844long double 845fmodl (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 859float cabsf (float complex z); 860 861float 862cabsf (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 870double cabs (double complex z); 871 872double 873cabs (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 881long double cabsl (long double complex z); 882 883long double 884cabsl (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 893float cargf (float complex z); 894 895float 896cargf (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 904double carg (double complex z); 905 906double 907carg (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 915long double cargl (long double complex z); 916 917long double 918cargl (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 928float complex cexpf (float complex z); 929 930float complex 931cexpf (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 945double complex cexp (double complex z); 946 947double complex 948cexp (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 962long double complex cexpl (long double complex z); 963 964long double complex 965cexpl (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 981float complex clogf (float complex z); 982 983float complex 984clogf (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 995double complex clog (double complex z); 996 997double complex 998clog (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 1009long double complex clogl (long double complex z); 1010 1011long double complex 1012clogl (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 1025float complex clog10f (float complex z); 1026 1027float complex 1028clog10f (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 1039double complex clog10 (double complex z); 1040 1041double complex 1042clog10 (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 1053long double complex clog10l (long double complex z); 1054 1055long double complex 1056clog10l (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 1069float complex cpowf (float complex base, float complex power); 1070 1071float complex 1072cpowf (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 1080double complex cpow (double complex base, double complex power); 1081 1082double complex 1083cpow (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 1091long double complex cpowl (long double complex base, long double complex power); 1092 1093long double complex 1094cpowl (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 1104float complex csqrtf (float complex z); 1105 1106float complex 1107csqrtf (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 1159double complex csqrt (double complex z); 1160 1161double complex 1162csqrt (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 1214long double complex csqrtl (long double complex z); 1215 1216long double complex 1217csqrtl (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 1271float complex csinhf (float complex a); 1272 1273float complex 1274csinhf (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 1288double complex csinh (double complex a); 1289 1290double complex 1291csinh (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 1305long double complex csinhl (long double complex a); 1306 1307long double complex 1308csinhl (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 1324float complex ccoshf (float complex a); 1325 1326float complex 1327ccoshf (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 1341double complex ccosh (double complex a); 1342 1343double complex 1344ccosh (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 1358long double complex ccoshl (long double complex a); 1359 1360long double complex 1361ccoshl (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 1377float complex ctanhf (float complex a); 1378 1379float complex 1380ctanhf (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 1396double complex ctanh (double complex a); 1397double complex 1398ctanh (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 1414long double complex ctanhl (long double complex a); 1415 1416long double complex 1417ctanhl (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 1435float complex csinf (float complex a); 1436 1437float complex 1438csinf (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 1452double complex csin (double complex a); 1453 1454double complex 1455csin (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 1469long double complex csinl (long double complex a); 1470 1471long double complex 1472csinl (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 1488float complex ccosf (float complex a); 1489 1490float complex 1491ccosf (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 1505double complex ccos (double complex a); 1506 1507double complex 1508ccos (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 1522long double complex ccosl (long double complex a); 1523 1524long double complex 1525ccosl (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 1541float complex ctanf (float complex a); 1542 1543float complex 1544ctanf (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 1560double complex ctan (double complex a); 1561 1562double complex 1563ctan (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 1579long double complex ctanl (long double complex a); 1580 1581long double complex 1582ctanl (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 1602complex float casinf (complex float z); 1603 1604complex float 1605casinf (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 1614complex double casin (complex double z); 1615 1616complex double 1617casin (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 1626complex long double casinl (complex long double z); 1627 1628complex long double 1629casinl (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 1641complex float cacosf (complex float z); 1642 1643complex float 1644cacosf (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 1653complex double cacos (complex double z); 1654 1655complex double 1656cacos (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 1665complex long double cacosl (complex long double z); 1666 1667complex long double 1668cacosl (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 1680complex float catanf (complex float z); 1681 1682complex float 1683catanf (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 1692complex double catan (complex double z); 1693 1694complex double 1695catan (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 1704complex long double catanl (complex long double z); 1705 1706complex long double 1707catanl (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 1719complex float casinhf (complex float z); 1720 1721complex float 1722casinhf (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 1731complex double casinh (complex double z); 1732 1733complex double 1734casinh (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 1743complex long double casinhl (complex long double z); 1744 1745complex long double 1746casinhl (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 1758complex float cacoshf (complex float z); 1759 1760complex float 1761cacoshf (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 1770complex double cacosh (complex double z); 1771 1772complex double 1773cacosh (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 1782complex long double cacoshl (complex long double z); 1783 1784complex long double 1785cacoshl (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 1797complex float catanhf (complex float z); 1798 1799complex float 1800catanhf (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 1809complex double catanh (complex double z); 1810 1811complex double 1812catanh (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 1820complex long double catanhl (complex long double z); 1821 1822complex long double 1823catanhl (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 1832double 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 1843double 1844tgamma (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 1970double 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 1979double 1980lgamma (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 2140float tgammaf (float); 2141 2142float 2143tgammaf (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 2151float lgammaf (float); 2152 2153float 2154lgammaf (float x) 2155{ 2156 return (float) lgamma ((double) x); 2157} 2158#endif 2159 2160#ifndef HAVE_FMA 2161#define HAVE_FMA 1 2162double fma (double, double, double); 2163 2164double 2165fma (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 2173float fmaf (float, float, float); 2174 2175float 2176fmaf (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 2184long double fmal (long double, long double, long double); 2185 2186long double 2187fmal (long double x, long double y, long double z) 2188{ 2189 return x * y + z; 2190} 2191#endif 2192