libgfortran/intrinsics/c99_functions.c

   1 /* Implementation of various C99 functions
   2    Copyright (C) 2004 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 2 of the License, or (at your option) any later version.
  10
  11 In addition to the permissions in the GNU General Public License, the
  12 Free Software Foundation gives you unlimited permission to link the
  13 compiled version of this file into combinations with other programs,
  14 and to distribute those combinations without any restriction coming
  15 from the use of this file.  (The General Public License restrictions
  16 do apply in other respects; for example, they cover modification of
  17 the file, and distribution when not linked into a combine
  18 executable.)
  19
  20 Libgfortran is distributed in the hope that it will be useful,
  21 but WITHOUT ANY WARRANTY; without even the implied warranty of
  22 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  23 GNU General Public License for more details.
  24
  25 You should have received a copy of the GNU General Public
  26 License along with libgfortran; see the file COPYING.  If not,
  27 write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
  28 Boston, MA 02110-1301, USA.  */
  29
  30 #include "config.h"
  31
  32 #define C99_PROTOS_H WE_DONT_WANT_PROTOS_NOW
  33 #include "libgfortran.h"
  34
  35 /* IRIX's <math.h> declares a non-C99 compliant implementation of cabs,
  36    which takes two floating point arguments instead of a single complex.
  37    If <complex.h> is missing this prevents building of c99_functions.c.
  38    To work around this we redirect cabs{,f,l} calls to __gfc_cabs{,f,l}.  */
  39
  40 #if defined(__sgi__) && !defined(HAVE_COMPLEX_H)
  41 #undef HAVE_CABS
  42 #undef HAVE_CABSF
  43 #undef HAVE_CABSL
  44 #define cabs __gfc_cabs
  45 #define cabsf __gfc_cabsf
  46 #define cabsl __gfc_cabsl
  47 #endif
  48
  49 /* Tru64's <math.h> declares a non-C99 compliant implementation of cabs,
  50    which takes two floating point arguments instead of a single complex.
  51    To work around this we redirect cabs{,f,l} calls to __gfc_cabs{,f,l}.  */
  52
  53 #ifdef __osf__
  54 #undef HAVE_CABS
  55 #undef HAVE_CABSF
  56 #undef HAVE_CABSL
  57 #define cabs __gfc_cabs
  58 #define cabsf __gfc_cabsf
  59 #define cabsl __gfc_cabsl
  60 #endif
  61
  62 /* Prototypes to silence -Wstrict-prototypes -Wmissing-prototypes.  */
  63
  64 float cabsf(float complex);
  65 double cabs(double complex);
  66 long double cabsl(long double complex);
  67
  68 float cargf(float complex);
  69 double carg(double complex);
  70 long double cargl(long double complex);
  71
  72 float complex clog10f(float complex);
  73 double complex clog10(double complex);
  74 long double complex clog10l(long double complex);
  75
  76
  77 /* Wrappers for systems without the various C99 single precision Bessel
  78    functions.  */
  79
  80 #if defined(HAVE_J0) && ! defined(HAVE_J0F)
  81 #define HAVE_J0F 1
  82 extern float j0f (float);
  83
  84 float
  85 j0f (float x)
  86 {
  87   return (float) j0 ((double) x);
  88 }
  89 #endif
  90
  91 #if defined(HAVE_J1) && !defined(HAVE_J1F)
  92 #define HAVE_J1F 1
  93 extern float j1f (float);
  94
  95 float j1f (float x)
  96 {
  97   return (float) j1 ((double) x);
  98 }
  99 #endif
 100
 101 #if defined(HAVE_JN) && !defined(HAVE_JNF)
 102 #define HAVE_JNF 1
 103 extern float jnf (int, float);
 104
 105 float
 106 jnf (int n, float x)
 107 {
 108   return (float) jn (n, (double) x);
 109 }
 110 #endif
 111
 112 #if defined(HAVE_Y0) && !defined(HAVE_Y0F)
 113 #define HAVE_Y0F 1
 114 extern float y0f (float);
 115
 116 float
 117 y0f (float x)
 118 {
 119   return (float) y0 ((double) x);
 120 }
 121 #endif
 122
 123 #if defined(HAVE_Y1) && !defined(HAVE_Y1F)
 124 #define HAVE_Y1F 1
 125 extern float y1f (float);
 126
 127 float
 128 y1f (float x)
 129 {
 130   return (float) y1 ((double) x);
 131 }
 132 #endif
 133
 134 #if defined(HAVE_YN) && !defined(HAVE_YNF)
 135 #define HAVE_YNF 1
 136 extern float ynf (int, float);
 137
 138 float
 139 ynf (int n, float x)
 140 {
 141   return (float) yn (n, (double) x);
 142 }
 143 #endif
 144
 145
 146 /* Wrappers for systems without the C99 erff() and erfcf() functions.  */
 147
 148 #if defined(HAVE_ERF) && !defined(HAVE_ERFF)
 149 #define HAVE_ERFF 1
 150 extern float erff (float);
 151
 152 float
 153 erff (float x)
 154 {
 155   return (float) erf ((double) x);
 156 }
 157 #endif
 158
 159 #if defined(HAVE_ERFC) && !defined(HAVE_ERFCF)
 160 #define HAVE_ERFCF 1
 161 extern float erfcf (float);
 162
 163 float
 164 erfcf (float x)
 165 {
 166   return (float) erfc ((double) x);
 167 }
 168 #endif
 169
 170
 171 #ifndef HAVE_ACOSF
 172 #define HAVE_ACOSF 1
 173 float
 174 acosf(float x)
 175 {
 176   return (float) acos(x);
 177 }
 178 #endif
 179
 180 #if HAVE_ACOSH && !HAVE_ACOSHF
 181 float
 182 acoshf (float x)
 183 {
 184   return (float) acosh ((double) x);
 185 }
 186 #endif
 187
 188 #ifndef HAVE_ASINF
 189 #define HAVE_ASINF 1
 190 float
 191 asinf(float x)
 192 {
 193   return (float) asin(x);
 194 }
 195 #endif
 196
 197 #if HAVE_ASINH && !HAVE_ASINHF
 198 float
 199 asinhf (float x)
 200 {
 201   return (float) asinh ((double) x);
 202 }
 203 #endif
 204
 205 #ifndef HAVE_ATAN2F
 206 #define HAVE_ATAN2F 1
 207 float
 208 atan2f(float y, float x)
 209 {
 210   return (float) atan2(y, x);
 211 }
 212 #endif
 213
 214 #ifndef HAVE_ATANF
 215 #define HAVE_ATANF 1
 216 float
 217 atanf(float x)
 218 {
 219   return (float) atan(x);
 220 }
 221 #endif
 222
 223 #if HAVE_ATANH && !HAVE_ATANHF
 224 float
 225 atanhf (float x)
 226 {
 227   return (float) atanh ((double) x);
 228 }
 229 #endif
 230
 231 #ifndef HAVE_CEILF
 232 #define HAVE_CEILF 1
 233 float
 234 ceilf(float x)
 235 {
 236   return (float) ceil(x);
 237 }
 238 #endif
 239
 240 #ifndef HAVE_COPYSIGNF
 241 #define HAVE_COPYSIGNF 1
 242 float
 243 copysignf(float x, float y)
 244 {
 245   return (float) copysign(x, y);
 246 }
 247 #endif
 248
 249 #ifndef HAVE_COSF
 250 #define HAVE_COSF 1
 251 float
 252 cosf(float x)
 253 {
 254   return (float) cos(x);
 255 }
 256 #endif
 257
 258 #ifndef HAVE_COSHF
 259 #define HAVE_COSHF 1
 260 float
 261 coshf(float x)
 262 {
 263   return (float) cosh(x);
 264 }
 265 #endif
 266
 267 #ifndef HAVE_EXPF
 268 #define HAVE_EXPF 1
 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
 279 fabsf(float x)
 280 {
 281   return (float) fabs(x);
 282 }
 283 #endif
 284
 285 #ifndef HAVE_FLOORF
 286 #define HAVE_FLOORF 1
 287 float
 288 floorf(float x)
 289 {
 290   return (float) floor(x);
 291 }
 292 #endif
 293
 294 #ifndef HAVE_FMODF
 295 #define HAVE_FMODF 1
 296 float
 297 fmodf (float x, float y)
 298 {
 299   return (float) fmod (x, y);
 300 }
 301 #endif
 302
 303 #ifndef HAVE_FREXPF
 304 #define HAVE_FREXPF 1
 305 float
 306 frexpf(float x, int *exp)
 307 {
 308   return (float) frexp(x, exp);
 309 }
 310 #endif
 311
 312 #ifndef HAVE_HYPOTF
 313 #define HAVE_HYPOTF 1
 314 float
 315 hypotf(float x, float y)
 316 {
 317   return (float) hypot(x, y);
 318 }
 319 #endif
 320
 321 #ifndef HAVE_LOGF
 322 #define HAVE_LOGF 1
 323 float
 324 logf(float x)
 325 {
 326   return (float) log(x);
 327 }
 328 #endif
 329
 330 #ifndef HAVE_LOG10F
 331 #define HAVE_LOG10F 1
 332 float
 333 log10f(float x)
 334 {
 335   return (float) log10(x);
 336 }
 337 #endif
 338
 339 #ifndef HAVE_SCALBN
 340 #define HAVE_SCALBN 1
 341 double
 342 scalbn(double x, int y)
 343 {
 344   return x * pow(FLT_RADIX, y);
 345 }
 346 #endif
 347
 348 #ifndef HAVE_SCALBNF
 349 #define HAVE_SCALBNF 1
 350 float
 351 scalbnf(float x, int y)
 352 {
 353   return (float) scalbn(x, y);
 354 }
 355 #endif
 356
 357 #ifndef HAVE_SINF
 358 #define HAVE_SINF 1
 359 float
 360 sinf(float x)
 361 {
 362   return (float) sin(x);
 363 }
 364 #endif
 365
 366 #ifndef HAVE_SINHF
 367 #define HAVE_SINHF 1
 368 float
 369 sinhf(float x)
 370 {
 371   return (float) sinh(x);
 372 }
 373 #endif
 374
 375 #ifndef HAVE_SQRTF
 376 #define HAVE_SQRTF 1
 377 float
 378 sqrtf(float x)
 379 {
 380   return (float) sqrt(x);
 381 }
 382 #endif
 383
 384 #ifndef HAVE_TANF
 385 #define HAVE_TANF 1
 386 float
 387 tanf(float x)
 388 {
 389   return (float) tan(x);
 390 }
 391 #endif
 392
 393 #ifndef HAVE_TANHF
 394 #define HAVE_TANHF 1
 395 float
 396 tanhf(float x)
 397 {
 398   return (float) tanh(x);
 399 }
 400 #endif
 401
 402 #ifndef HAVE_TRUNC
 403 #define HAVE_TRUNC 1
 404 double
 405 trunc(double x)
 406 {
 407   if (!isfinite (x))
 408     return x;
 409
 410   if (x < 0.0)
 411     return - floor (-x);
 412   else
 413     return floor (x);
 414 }
 415 #endif
 416
 417 #ifndef HAVE_TRUNCF
 418 #define HAVE_TRUNCF 1
 419 float
 420 truncf(float x)
 421 {
 422   return (float) trunc (x);
 423 }
 424 #endif
 425
 426 #ifndef HAVE_NEXTAFTERF
 427 #define HAVE_NEXTAFTERF 1
 428 /* This is a portable implementation of nextafterf that is intended to be
 429    independent of the floating point format or its in memory representation.
 430    This implementation works correctly with denormalized values.  */
 431 float
 432 nextafterf(float x, float y)
 433 {
 434   /* This variable is marked volatile to avoid excess precision problems
 435      on some platforms, including IA-32.  */
 436   volatile float delta;
 437   float absx, denorm_min;
 438
 439   if (isnan(x) || isnan(y))
 440     return x + y;
 441   if (x == y)
 442     return x;
 443   if (!isfinite (x))
 444     return x > 0 ? __FLT_MAX__ : - __FLT_MAX__;
 445
 446   /* absx = fabsf (x);  */
 447   absx = (x < 0.0) ? -x : x;
 448
 449   /* __FLT_DENORM_MIN__ is non-zero iff the target supports denormals.  */
 450   if (__FLT_DENORM_MIN__ == 0.0f)
 451     denorm_min = __FLT_MIN__;
 452   else
 453     denorm_min = __FLT_DENORM_MIN__;
 454
 455   if (absx < __FLT_MIN__)
 456     delta = denorm_min;
 457   else
 458     {
 459       float frac;
 460       int exp;
 461
 462       /* Discard the fraction from x.  */
 463       frac = frexpf (absx, &exp);
 464       delta = scalbnf (0.5f, exp);
 465
 466       /* Scale x by the epsilon of the representation.  By rights we should
 467          have been able to combine this with scalbnf, but some targets don't
 468          get that correct with denormals.  */
 469       delta *= __FLT_EPSILON__;
 470
 471       /* If we're going to be reducing the absolute value of X, and doing so
 472          would reduce the exponent of X, then the delta to be applied is
 473          one exponent smaller.  */
 474       if (frac == 0.5f && (y < x) == (x > 0))
 475         delta *= 0.5f;
 476
 477       /* If that underflows to zero, then we're back to the minimum.  */
 478       if (delta == 0.0f)
 479         delta = denorm_min;
 480     }
 481
 482   if (y < x)
 483     delta = -delta;
 484
 485   return x + delta;
 486 }
 487 #endif
 488
 489
 490 #ifndef HAVE_POWF
 491 #define HAVE_POWF 1
 492 float
 493 powf(float x, float y)
 494 {
 495   return (float) pow(x, y);
 496 }
 497 #endif
 498
 499 /* Note that if fpclassify is not defined, then NaN is not handled */
 500
 501 /* Algorithm by Steven G. Kargl.  */
 502
 503 #if !defined(HAVE_ROUNDL)
 504 #define HAVE_ROUNDL 1
 505 #if defined(HAVE_CEILL)
 506 /* Round to nearest integral value.  If the argument is halfway between two
 507    integral values then round away from zero.  */
 508
 509 long double
 510 roundl(long double x)
 511 {
 512    long double t;
 513    if (!isfinite (x))
 514      return (x);
 515
 516    if (x >= 0.0)
 517     {
 518       t = ceill(x);
 519       if (t - x > 0.5)
 520         t -= 1.0;
 521       return (t);
 522     }
 523    else
 524     {
 525       t = ceill(-x);
 526       if (t + x > 0.5)
 527         t -= 1.0;
 528       return (-t);
 529     }
 530 }
 531 #else
 532
 533 /* Poor version of roundl for system that don't have ceill.  */
 534 long double
 535 roundl(long double x)
 536 {
 537   if (x > DBL_MAX || x < -DBL_MAX)
 538     {
 539 #ifdef HAVE_NEXTAFTERL
 540       static long double prechalf = nexafterl (0.5L, LDBL_MAX);
 541 #else
 542       static long double prechalf = 0.5L;
 543 #endif
 544       return (GFC_INTEGER_LARGEST) (x + (x > 0 ? prechalf : -prechalf));
 545     }
 546   else
 547     /* Use round().  */
 548     return round((double) x);
 549 }
 550
 551 #endif
 552 #endif
 553
 554 #ifndef HAVE_ROUND
 555 #define HAVE_ROUND 1
 556 /* Round to nearest integral value.  If the argument is halfway between two
 557    integral values then round away from zero.  */
 558
 559 double
 560 round(double x)
 561 {
 562    double t;
 563    if (!isfinite (x))
 564      return (x);
 565
 566    if (x >= 0.0)
 567     {
 568       t = ceil(x);
 569       if (t - x > 0.5)
 570         t -= 1.0;
 571       return (t);
 572     }
 573    else
 574     {
 575       t = ceil(-x);
 576       if (t + x > 0.5)
 577         t -= 1.0;
 578       return (-t);
 579     }
 580 }
 581 #endif
 582
 583 #ifndef HAVE_ROUNDF
 584 #define HAVE_ROUNDF 1
 585 /* Round to nearest integral value.  If the argument is halfway between two
 586    integral values then round away from zero.  */
 587
 588 float
 589 roundf(float x)
 590 {
 591    float t;
 592    if (!isfinite (x))
 593      return (x);
 594
 595    if (x >= 0.0)
 596     {
 597       t = ceilf(x);
 598       if (t - x > 0.5)
 599         t -= 1.0;
 600       return (t);
 601     }
 602    else
 603     {
 604       t = ceilf(-x);
 605       if (t + x > 0.5)
 606         t -= 1.0;
 607       return (-t);
 608     }
 609 }
 610 #endif
 611
 612
 613 /* lround{f,,l} and llround{f,,l} functions.  */
 614
 615 #if !defined(HAVE_LROUNDF) && defined(HAVE_ROUNDF)
 616 #define HAVE_LROUNDF 1
 617 long int
 618 lroundf (float x)
 619 {
 620   return (long int) roundf (x);
 621 }
 622 #endif
 623
 624 #if !defined(HAVE_LROUND) && defined(HAVE_ROUND)
 625 #define HAVE_LROUND 1
 626 long int
 627 lround (double x)
 628 {
 629   return (long int) round (x);
 630 }
 631 #endif
 632
 633 #if !defined(HAVE_LROUNDL) && defined(HAVE_ROUNDL)
 634 #define HAVE_LROUNDL 1
 635 long int
 636 lroundl (long double x)
 637 {
 638   return (long long int) roundl (x);
 639 }
 640 #endif
 641
 642 #if !defined(HAVE_LLROUNDF) && defined(HAVE_ROUNDF)
 643 #define HAVE_LLROUNDF 1
 644 long long int
 645 llroundf (float x)
 646 {
 647   return (long long int) roundf (x);
 648 }
 649 #endif
 650
 651 #if !defined(HAVE_LLROUND) && defined(HAVE_ROUND)
 652 #define HAVE_LLROUND 1
 653 long long int
 654 llround (double x)
 655 {
 656   return (long long int) round (x);
 657 }
 658 #endif
 659
 660 #if !defined(HAVE_LLROUNDL) && defined(HAVE_ROUNDL)
 661 #define HAVE_LLROUNDL 1
 662 long long int
 663 llroundl (long double x)
 664 {
 665   return (long long int) roundl (x);
 666 }
 667 #endif
 668
 669
 670 #ifndef HAVE_LOG10L
 671 #define HAVE_LOG10L 1
 672 /* log10 function for long double variables. The version provided here
 673    reduces the argument until it fits into a double, then use log10.  */
 674 long double
 675 log10l(long double x)
 676 {
 677 #if LDBL_MAX_EXP > DBL_MAX_EXP
 678   if (x > DBL_MAX)
 679     {
 680       double val;
 681       int p2_result = 0;
 682       if (x > 0x1p16383L) { p2_result += 16383; x /= 0x1p16383L; }
 683       if (x > 0x1p8191L) { p2_result += 8191; x /= 0x1p8191L; }
 684       if (x > 0x1p4095L) { p2_result += 4095; x /= 0x1p4095L; }
 685       if (x > 0x1p2047L) { p2_result += 2047; x /= 0x1p2047L; }
 686       if (x > 0x1p1023L) { p2_result += 1023; x /= 0x1p1023L; }
 687       val = log10 ((double) x);
 688       return (val + p2_result * .30102999566398119521373889472449302L);
 689     }
 690 #endif
 691 #if LDBL_MIN_EXP < DBL_MIN_EXP
 692   if (x < DBL_MIN)
 693     {
 694       double val;
 695       int p2_result = 0;
 696       if (x < 0x1p-16380L) { p2_result += 16380; x /= 0x1p-16380L; }
 697       if (x < 0x1p-8189L) { p2_result += 8189; x /= 0x1p-8189L; }
 698       if (x < 0x1p-4093L) { p2_result += 4093; x /= 0x1p-4093L; }
 699       if (x < 0x1p-2045L) { p2_result += 2045; x /= 0x1p-2045L; }
 700       if (x < 0x1p-1021L) { p2_result += 1021; x /= 0x1p-1021L; }
 701       val = fabs(log10 ((double) x));
 702       return (- val - p2_result * .30102999566398119521373889472449302L);
 703     }
 704 #endif
 705     return log10 (x);
 706 }
 707 #endif
 708
 709
 710 #ifndef HAVE_FLOORL
 711 #define HAVE_FLOORL 1
 712 long double
 713 floorl (long double x)
 714 {
 715   /* Zero, possibly signed.  */
 716   if (x == 0)
 717     return x;
 718
 719   /* Large magnitude.  */
 720   if (x > DBL_MAX || x < (-DBL_MAX))
 721     return x;
 722
 723   /* Small positive values.  */
 724   if (x >= 0 && x < DBL_MIN)
 725     return 0;
 726
 727   /* Small negative values.  */
 728   if (x < 0 && x > (-DBL_MIN))
 729     return -1;
 730
 731   return floor (x);
 732 }
 733 #endif
 734
 735
 736 #ifndef HAVE_FMODL
 737 #define HAVE_FMODL 1
 738 long double
 739 fmodl (long double x, long double y)
 740 {
 741   if (y == 0.0L)
 742     return 0.0L;
 743
 744   /* Need to check that the result has the same sign as x and magnitude
 745      less than the magnitude of y.  */
 746   return x - floorl (x / y) * y;
 747 }
 748 #endif
 749
 750
 751 #if !defined(HAVE_CABSF)
 752 #define HAVE_CABSF 1
 753 float
 754 cabsf (float complex z)
 755 {
 756   return hypotf (REALPART (z), IMAGPART (z));
 757 }
 758 #endif
 759
 760 #if !defined(HAVE_CABS)
 761 #define HAVE_CABS 1
 762 double
 763 cabs (double complex z)
 764 {
 765   return hypot (REALPART (z), IMAGPART (z));
 766 }
 767 #endif
 768
 769 #if !defined(HAVE_CABSL) && defined(HAVE_HYPOTL)
 770 #define HAVE_CABSL 1
 771 long double
 772 cabsl (long double complex z)
 773 {
 774   return hypotl (REALPART (z), IMAGPART (z));
 775 }
 776 #endif
 777
 778
 779 #if !defined(HAVE_CARGF)
 780 #define HAVE_CARGF 1
 781 float
 782 cargf (float complex z)
 783 {
 784   return atan2f (IMAGPART (z), REALPART (z));
 785 }
 786 #endif
 787
 788 #if !defined(HAVE_CARG)
 789 #define HAVE_CARG 1
 790 double
 791 carg (double complex z)
 792 {
 793   return atan2 (IMAGPART (z), REALPART (z));
 794 }
 795 #endif
 796
 797 #if !defined(HAVE_CARGL) && defined(HAVE_ATAN2L)
 798 #define HAVE_CARGL 1
 799 long double
 800 cargl (long double complex z)
 801 {
 802   return atan2l (IMAGPART (z), REALPART (z));
 803 }
 804 #endif
 805
 806
 807 /* exp(z) = exp(a)*(cos(b) + i sin(b))  */
 808 #if !defined(HAVE_CEXPF)
 809 #define HAVE_CEXPF 1
 810 float complex
 811 cexpf (float complex z)
 812 {
 813   float a, b;
 814   float complex v;
 815
 816   a = REALPART (z);
 817   b = IMAGPART (z);
 818   COMPLEX_ASSIGN (v, cosf (b), sinf (b));
 819   return expf (a) * v;
 820 }
 821 #endif
 822
 823 #if !defined(HAVE_CEXP)
 824 #define HAVE_CEXP 1
 825 double complex
 826 cexp (double complex z)
 827 {
 828   double a, b;
 829   double complex v;
 830
 831   a = REALPART (z);
 832   b = IMAGPART (z);
 833   COMPLEX_ASSIGN (v, cos (b), sin (b));
 834   return exp (a) * v;
 835 }
 836 #endif
 837
 838 #if !defined(HAVE_CEXPL) && defined(HAVE_COSL) && defined(HAVE_SINL) && defined(EXPL)
 839 #define HAVE_CEXPL 1
 840 long double complex
 841 cexpl (long double complex z)
 842 {
 843   long double a, b;
 844   long double complex v;
 845
 846   a = REALPART (z);
 847   b = IMAGPART (z);
 848   COMPLEX_ASSIGN (v, cosl (b), sinl (b));
 849   return expl (a) * v;
 850 }
 851 #endif
 852
 853
 854 /* log(z) = log (cabs(z)) + i*carg(z)  */
 855 #if !defined(HAVE_CLOGF)
 856 #define HAVE_CLOGF 1
 857 float complex
 858 clogf (float complex z)
 859 {
 860   float complex v;
 861
 862   COMPLEX_ASSIGN (v, logf (cabsf (z)), cargf (z));
 863   return v;
 864 }
 865 #endif
 866
 867 #if !defined(HAVE_CLOG)
 868 #define HAVE_CLOG 1
 869 double complex
 870 clog (double complex z)
 871 {
 872   double complex v;
 873
 874   COMPLEX_ASSIGN (v, log (cabs (z)), carg (z));
 875   return v;
 876 }
 877 #endif
 878
 879 #if !defined(HAVE_CLOGL) && defined(HAVE_LOGL) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
 880 #define HAVE_CLOGL 1
 881 long double complex
 882 clogl (long double complex z)
 883 {
 884   long double complex v;
 885
 886   COMPLEX_ASSIGN (v, logl (cabsl (z)), cargl (z));
 887   return v;
 888 }
 889 #endif
 890
 891
 892 /* log10(z) = log10 (cabs(z)) + i*carg(z)  */
 893 #if !defined(HAVE_CLOG10F)
 894 #define HAVE_CLOG10F 1
 895 float complex
 896 clog10f (float complex z)
 897 {
 898   float complex v;
 899
 900   COMPLEX_ASSIGN (v, log10f (cabsf (z)), cargf (z));
 901   return v;
 902 }
 903 #endif
 904
 905 #if !defined(HAVE_CLOG10)
 906 #define HAVE_CLOG10 1
 907 double complex
 908 clog10 (double complex z)
 909 {
 910   double complex v;
 911
 912   COMPLEX_ASSIGN (v, log10 (cabs (z)), carg (z));
 913   return v;
 914 }
 915 #endif
 916
 917 #if !defined(HAVE_CLOG10L) && defined(HAVE_LOG10L) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
 918 #define HAVE_CLOG10L 1
 919 long double complex
 920 clog10l (long double complex z)
 921 {
 922   long double complex v;
 923
 924   COMPLEX_ASSIGN (v, log10l (cabsl (z)), cargl (z));
 925   return v;
 926 }
 927 #endif
 928
 929
 930 /* pow(base, power) = cexp (power * clog (base))  */
 931 #if !defined(HAVE_CPOWF)
 932 #define HAVE_CPOWF 1
 933 float complex
 934 cpowf (float complex base, float complex power)
 935 {
 936   return cexpf (power * clogf (base));
 937 }
 938 #endif
 939
 940 #if !defined(HAVE_CPOW)
 941 #define HAVE_CPOW 1
 942 double complex
 943 cpow (double complex base, double complex power)
 944 {
 945   return cexp (power * clog (base));
 946 }
 947 #endif
 948
 949 #if !defined(HAVE_CPOWL) && defined(HAVE_CEXPL) && defined(HAVE_CLOGL)
 950 #define HAVE_CPOWL 1
 951 long double complex
 952 cpowl (long double complex base, long double complex power)
 953 {
 954   return cexpl (power * clogl (base));
 955 }
 956 #endif
 957
 958
 959 /* sqrt(z).  Algorithm pulled from glibc.  */
 960 #if !defined(HAVE_CSQRTF)
 961 #define HAVE_CSQRTF 1
 962 float complex
 963 csqrtf (float complex z)
 964 {
 965   float re, im;
 966   float complex v;
 967
 968   re = REALPART (z);
 969   im = IMAGPART (z);
 970   if (im == 0)
 971     {
 972       if (re < 0)
 973         {
 974           COMPLEX_ASSIGN (v, 0, copysignf (sqrtf (-re), im));
 975         }
 976       else
 977         {
 978           COMPLEX_ASSIGN (v, fabsf (sqrtf (re)), copysignf (0, im));
 979         }
 980     }
 981   else if (re == 0)
 982     {
 983       float r;
 984
 985       r = sqrtf (0.5 * fabsf (im));
 986
 987       COMPLEX_ASSIGN (v, r, copysignf (r, im));
 988     }
 989   else
 990     {
 991       float d, r, s;
 992
 993       d = hypotf (re, im);
 994       /* Use the identity   2  Re res  Im res = Im x
 995          to avoid cancellation error in  d +/- Re x.  */
 996       if (re > 0)
 997         {
 998           r = sqrtf (0.5 * d + 0.5 * re);
 999           s = (0.5 * im) / r;
1000         }
1001       else
1002         {
1003           s = sqrtf (0.5 * d - 0.5 * re);
1004           r = fabsf ((0.5 * im) / s);
1005         }
1006
1007       COMPLEX_ASSIGN (v, r, copysignf (s, im));
1008     }
1009   return v;
1010 }
1011 #endif
1012
1013 #if !defined(HAVE_CSQRT)
1014 #define HAVE_CSQRT 1
1015 double complex
1016 csqrt (double complex z)
1017 {
1018   double re, im;
1019   double complex v;
1020
1021   re = REALPART (z);
1022   im = IMAGPART (z);
1023   if (im == 0)
1024     {
1025       if (re < 0)
1026         {
1027           COMPLEX_ASSIGN (v, 0, copysign (sqrt (-re), im));
1028         }
1029       else
1030         {
1031           COMPLEX_ASSIGN (v, fabs (sqrt (re)), copysign (0, im));
1032         }
1033     }
1034   else if (re == 0)
1035     {
1036       double r;
1037
1038       r = sqrt (0.5 * fabs (im));
1039
1040       COMPLEX_ASSIGN (v, r, copysign (r, im));
1041     }
1042   else
1043     {
1044       double d, r, s;
1045
1046       d = hypot (re, im);
1047       /* Use the identity   2  Re res  Im res = Im x
1048          to avoid cancellation error in  d +/- Re x.  */
1049       if (re > 0)
1050         {
1051           r = sqrt (0.5 * d + 0.5 * re);
1052           s = (0.5 * im) / r;
1053         }
1054       else
1055         {
1056           s = sqrt (0.5 * d - 0.5 * re);
1057           r = fabs ((0.5 * im) / s);
1058         }
1059
1060       COMPLEX_ASSIGN (v, r, copysign (s, im));
1061     }
1062   return v;
1063 }
1064 #endif
1065
1066 #if !defined(HAVE_CSQRTL) && defined(HAVE_COPYSIGNL) && defined(HAVE_SQRTL) && defined(HAVE_FABSL) && defined(HAVE_HYPOTL)
1067 #define HAVE_CSQRTL 1
1068 long double complex
1069 csqrtl (long double complex z)
1070 {
1071   long double re, im;
1072   long double complex v;
1073
1074   re = REALPART (z);
1075   im = IMAGPART (z);
1076   if (im == 0)
1077     {
1078       if (re < 0)
1079         {
1080           COMPLEX_ASSIGN (v, 0, copysignl (sqrtl (-re), im));
1081         }
1082       else
1083         {
1084           COMPLEX_ASSIGN (v, fabsl (sqrtl (re)), copysignl (0, im));
1085         }
1086     }
1087   else if (re == 0)
1088     {
1089       long double r;
1090
1091       r = sqrtl (0.5 * fabsl (im));
1092
1093       COMPLEX_ASSIGN (v, copysignl (r, im), r);
1094     }
1095   else
1096     {
1097       long double d, r, s;
1098
1099       d = hypotl (re, im);
1100       /* Use the identity   2  Re res  Im res = Im x
1101          to avoid cancellation error in  d +/- Re x.  */
1102       if (re > 0)
1103         {
1104           r = sqrtl (0.5 * d + 0.5 * re);
1105           s = (0.5 * im) / r;
1106         }
1107       else
1108         {
1109           s = sqrtl (0.5 * d - 0.5 * re);
1110           r = fabsl ((0.5 * im) / s);
1111         }
1112
1113       COMPLEX_ASSIGN (v, r, copysignl (s, im));
1114     }
1115   return v;
1116 }
1117 #endif
1118
1119
1120 /* sinh(a + i b) = sinh(a) cos(b) + i cosh(a) sin(b)  */
1121 #if !defined(HAVE_CSINHF)
1122 #define HAVE_CSINHF 1
1123 float complex
1124 csinhf (float complex a)
1125 {
1126   float r, i;
1127   float complex v;
1128
1129   r = REALPART (a);
1130   i = IMAGPART (a);
1131   COMPLEX_ASSIGN (v, sinhf (r) * cosf (i), coshf (r) * sinf (i));
1132   return v;
1133 }
1134 #endif
1135
1136 #if !defined(HAVE_CSINH)
1137 #define HAVE_CSINH 1
1138 double complex
1139 csinh (double complex a)
1140 {
1141   double r, i;
1142   double complex v;
1143
1144   r = REALPART (a);
1145   i = IMAGPART (a);
1146   COMPLEX_ASSIGN (v, sinh (r) * cos (i), cosh (r) * sin (i));
1147   return v;
1148 }
1149 #endif
1150
1151 #if !defined(HAVE_CSINHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1152 #define HAVE_CSINHL 1
1153 long double complex
1154 csinhl (long double complex a)
1155 {
1156   long double r, i;
1157   long double complex v;
1158
1159   r = REALPART (a);
1160   i = IMAGPART (a);
1161   COMPLEX_ASSIGN (v, sinhl (r) * cosl (i), coshl (r) * sinl (i));
1162   return v;
1163 }
1164 #endif
1165
1166
1167 /* cosh(a + i b) = cosh(a) cos(b) - i sinh(a) sin(b)  */
1168 #if !defined(HAVE_CCOSHF)
1169 #define HAVE_CCOSHF 1
1170 float complex
1171 ccoshf (float complex a)
1172 {
1173   float r, i;
1174   float complex v;
1175
1176   r = REALPART (a);
1177   i = IMAGPART (a);
1178   COMPLEX_ASSIGN (v, coshf (r) * cosf (i), - (sinhf (r) * sinf (i)));
1179   return v;
1180 }
1181 #endif
1182
1183 #if !defined(HAVE_CCOSH)
1184 #define HAVE_CCOSH 1
1185 double complex
1186 ccosh (double complex a)
1187 {
1188   double r, i;
1189   double complex v;
1190
1191   r = REALPART (a);
1192   i = IMAGPART (a);
1193   COMPLEX_ASSIGN (v, cosh (r) * cos (i), - (sinh (r) * sin (i)));
1194   return v;
1195 }
1196 #endif
1197
1198 #if !defined(HAVE_CCOSHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1199 #define HAVE_CCOSHL 1
1200 long double complex
1201 ccoshl (long double complex a)
1202 {
1203   long double r, i;
1204   long double complex v;
1205
1206   r = REALPART (a);
1207   i = IMAGPART (a);
1208   COMPLEX_ASSIGN (v, coshl (r) * cosl (i), - (sinhl (r) * sinl (i)));
1209   return v;
1210 }
1211 #endif
1212
1213
1214 /* tanh(a + i b) = (tanh(a) + i tan(b)) / (1 - i tanh(a) tan(b))  */
1215 #if !defined(HAVE_CTANHF)
1216 #define HAVE_CTANHF 1
1217 float complex
1218 ctanhf (float complex a)
1219 {
1220   float rt, it;
1221   float complex n, d;
1222
1223   rt = tanhf (REALPART (a));
1224   it = tanf (IMAGPART (a));
1225   COMPLEX_ASSIGN (n, rt, it);
1226   COMPLEX_ASSIGN (d, 1, - (rt * it));
1227
1228   return n / d;
1229 }
1230 #endif
1231
1232 #if !defined(HAVE_CTANH)
1233 #define HAVE_CTANH 1
1234 double complex
1235 ctanh (double complex a)
1236 {
1237   double rt, it;
1238   double complex n, d;
1239
1240   rt = tanh (REALPART (a));
1241   it = tan (IMAGPART (a));
1242   COMPLEX_ASSIGN (n, rt, it);
1243   COMPLEX_ASSIGN (d, 1, - (rt * it));
1244
1245   return n / d;
1246 }
1247 #endif
1248
1249 #if !defined(HAVE_CTANHL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1250 #define HAVE_CTANHL 1
1251 long double complex
1252 ctanhl (long double complex a)
1253 {
1254   long double rt, it;
1255   long double complex n, d;
1256
1257   rt = tanhl (REALPART (a));
1258   it = tanl (IMAGPART (a));
1259   COMPLEX_ASSIGN (n, rt, it);
1260   COMPLEX_ASSIGN (d, 1, - (rt * it));
1261
1262   return n / d;
1263 }
1264 #endif
1265
1266
1267 /* sin(a + i b) = sin(a) cosh(b) + i cos(a) sinh(b)  */
1268 #if !defined(HAVE_CSINF)
1269 #define HAVE_CSINF 1
1270 float complex
1271 csinf (float complex a)
1272 {
1273   float r, i;
1274   float complex v;
1275
1276   r = REALPART (a);
1277   i = IMAGPART (a);
1278   COMPLEX_ASSIGN (v, sinf (r) * coshf (i), cosf (r) * sinhf (i));
1279   return v;
1280 }
1281 #endif
1282
1283 #if !defined(HAVE_CSIN)
1284 #define HAVE_CSIN 1
1285 double complex
1286 csin (double complex a)
1287 {
1288   double r, i;
1289   double complex v;
1290
1291   r = REALPART (a);
1292   i = IMAGPART (a);
1293   COMPLEX_ASSIGN (v, sin (r) * cosh (i), cos (r) * sinh (i));
1294   return v;
1295 }
1296 #endif
1297
1298 #if !defined(HAVE_CSINL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1299 #define HAVE_CSINL 1
1300 long double complex
1301 csinl (long double complex a)
1302 {
1303   long double r, i;
1304   long double complex v;
1305
1306   r = REALPART (a);
1307   i = IMAGPART (a);
1308   COMPLEX_ASSIGN (v, sinl (r) * coshl (i), cosl (r) * sinhl (i));
1309   return v;
1310 }
1311 #endif
1312
1313
1314 /* cos(a + i b) = cos(a) cosh(b) - i sin(a) sinh(b)  */
1315 #if !defined(HAVE_CCOSF)
1316 #define HAVE_CCOSF 1
1317 float complex
1318 ccosf (float complex a)
1319 {
1320   float r, i;
1321   float complex v;
1322
1323   r = REALPART (a);
1324   i = IMAGPART (a);
1325   COMPLEX_ASSIGN (v, cosf (r) * coshf (i), - (sinf (r) * sinhf (i)));
1326   return v;
1327 }
1328 #endif
1329
1330 #if !defined(HAVE_CCOS)
1331 #define HAVE_CCOS 1
1332 double complex
1333 ccos (double complex a)
1334 {
1335   double r, i;
1336   double complex v;
1337
1338   r = REALPART (a);
1339   i = IMAGPART (a);
1340   COMPLEX_ASSIGN (v, cos (r) * cosh (i), - (sin (r) * sinh (i)));
1341   return v;
1342 }
1343 #endif
1344
1345 #if !defined(HAVE_CCOSL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1346 #define HAVE_CCOSL 1
1347 long double complex
1348 ccosl (long double complex a)
1349 {
1350   long double r, i;
1351   long double complex v;
1352
1353   r = REALPART (a);
1354   i = IMAGPART (a);
1355   COMPLEX_ASSIGN (v, cosl (r) * coshl (i), - (sinl (r) * sinhl (i)));
1356   return v;
1357 }
1358 #endif
1359
1360
1361 /* tan(a + i b) = (tan(a) + i tanh(b)) / (1 - i tan(a) tanh(b))  */
1362 #if !defined(HAVE_CTANF)
1363 #define HAVE_CTANF 1
1364 float complex
1365 ctanf (float complex a)
1366 {
1367   float rt, it;
1368   float complex n, d;
1369
1370   rt = tanf (REALPART (a));
1371   it = tanhf (IMAGPART (a));
1372   COMPLEX_ASSIGN (n, rt, it);
1373   COMPLEX_ASSIGN (d, 1, - (rt * it));
1374
1375   return n / d;
1376 }
1377 #endif
1378
1379 #if !defined(HAVE_CTAN)
1380 #define HAVE_CTAN 1
1381 double complex
1382 ctan (double complex a)
1383 {
1384   double rt, it;
1385   double complex n, d;
1386
1387   rt = tan (REALPART (a));
1388   it = tanh (IMAGPART (a));
1389   COMPLEX_ASSIGN (n, rt, it);
1390   COMPLEX_ASSIGN (d, 1, - (rt * it));
1391
1392   return n / d;
1393 }
1394 #endif
1395
1396 #if !defined(HAVE_CTANL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1397 #define HAVE_CTANL 1
1398 long double complex
1399 ctanl (long double complex a)
1400 {
1401   long double rt, it;
1402   long double complex n, d;
1403
1404   rt = tanl (REALPART (a));
1405   it = tanhl (IMAGPART (a));
1406   COMPLEX_ASSIGN (n, rt, it);
1407   COMPLEX_ASSIGN (d, 1, - (rt * it));
1408
1409   return n / d;
1410 }
1411 #endif
1412