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 #ifndef HAVE_ROUND
 504 #define HAVE_ROUND 1
 505 /* Round to nearest integral value.  If the argument is halfway between two
 506    integral values then round away from zero.  */
 507
 508 double
 509 round(double x)
 510 {
 511    double t;
 512    if (!isfinite (x))
 513      return (x);
 514
 515    if (x >= 0.0)
 516     {
 517       t = ceil(x);
 518       if (t - x > 0.5)
 519         t -= 1.0;
 520       return (t);
 521     }
 522    else
 523     {
 524       t = ceil(-x);
 525       if (t + x > 0.5)
 526         t -= 1.0;
 527       return (-t);
 528     }
 529 }
 530 #endif
 531
 532 #ifndef HAVE_ROUNDF
 533 #define HAVE_ROUNDF 1
 534 /* Round to nearest integral value.  If the argument is halfway between two
 535    integral values then round away from zero.  */
 536
 537 float
 538 roundf(float x)
 539 {
 540    float t;
 541    if (!isfinite (x))
 542      return (x);
 543
 544    if (x >= 0.0)
 545     {
 546       t = ceilf(x);
 547       if (t - x > 0.5)
 548         t -= 1.0;
 549       return (t);
 550     }
 551    else
 552     {
 553       t = ceilf(-x);
 554       if (t + x > 0.5)
 555         t -= 1.0;
 556       return (-t);
 557     }
 558 }
 559 #endif
 560
 561 #ifndef HAVE_LOG10L
 562 #define HAVE_LOG10L 1
 563 /* log10 function for long double variables. The version provided here
 564    reduces the argument until it fits into a double, then use log10.  */
 565 long double
 566 log10l(long double x)
 567 {
 568 #if LDBL_MAX_EXP > DBL_MAX_EXP
 569   if (x > DBL_MAX)
 570     {
 571       double val;
 572       int p2_result = 0;
 573       if (x > 0x1p16383L) { p2_result += 16383; x /= 0x1p16383L; }
 574       if (x > 0x1p8191L) { p2_result += 8191; x /= 0x1p8191L; }
 575       if (x > 0x1p4095L) { p2_result += 4095; x /= 0x1p4095L; }
 576       if (x > 0x1p2047L) { p2_result += 2047; x /= 0x1p2047L; }
 577       if (x > 0x1p1023L) { p2_result += 1023; x /= 0x1p1023L; }
 578       val = log10 ((double) x);
 579       return (val + p2_result * .30102999566398119521373889472449302L);
 580     }
 581 #endif
 582 #if LDBL_MIN_EXP < DBL_MIN_EXP
 583   if (x < DBL_MIN)
 584     {
 585       double val;
 586       int p2_result = 0;
 587       if (x < 0x1p-16380L) { p2_result += 16380; x /= 0x1p-16380L; }
 588       if (x < 0x1p-8189L) { p2_result += 8189; x /= 0x1p-8189L; }
 589       if (x < 0x1p-4093L) { p2_result += 4093; x /= 0x1p-4093L; }
 590       if (x < 0x1p-2045L) { p2_result += 2045; x /= 0x1p-2045L; }
 591       if (x < 0x1p-1021L) { p2_result += 1021; x /= 0x1p-1021L; }
 592       val = fabs(log10 ((double) x));
 593       return (- val - p2_result * .30102999566398119521373889472449302L);
 594     }
 595 #endif
 596     return log10 (x);
 597 }
 598 #endif
 599
 600
 601 #ifndef HAVE_FLOORL
 602 #define HAVE_FLOORL 1
 603 long double
 604 floorl (long double x)
 605 {
 606   /* Zero, possibly signed.  */
 607   if (x == 0)
 608     return x;
 609
 610   /* Large magnitude.  */
 611   if (x > DBL_MAX || x < (-DBL_MAX))
 612     return x;
 613
 614   /* Small positive values.  */
 615   if (x >= 0 && x < DBL_MIN)
 616     return 0;
 617
 618   /* Small negative values.  */
 619   if (x < 0 && x > (-DBL_MIN))
 620     return -1;
 621
 622   return floor (x);
 623 }
 624 #endif
 625
 626
 627 #ifndef HAVE_FMODL
 628 #define HAVE_FMODL 1
 629 long double
 630 fmodl (long double x, long double y)
 631 {
 632   if (y == 0.0L)
 633     return 0.0L;
 634
 635   /* Need to check that the result has the same sign as x and magnitude
 636      less than the magnitude of y.  */
 637   return x - floorl (x / y) * y;
 638 }
 639 #endif
 640
 641
 642 #if !defined(HAVE_CABSF)
 643 #define HAVE_CABSF 1
 644 float
 645 cabsf (float complex z)
 646 {
 647   return hypotf (REALPART (z), IMAGPART (z));
 648 }
 649 #endif
 650
 651 #if !defined(HAVE_CABS)
 652 #define HAVE_CABS 1
 653 double
 654 cabs (double complex z)
 655 {
 656   return hypot (REALPART (z), IMAGPART (z));
 657 }
 658 #endif
 659
 660 #if !defined(HAVE_CABSL) && defined(HAVE_HYPOTL)
 661 #define HAVE_CABSL 1
 662 long double
 663 cabsl (long double complex z)
 664 {
 665   return hypotl (REALPART (z), IMAGPART (z));
 666 }
 667 #endif
 668
 669
 670 #if !defined(HAVE_CARGF)
 671 #define HAVE_CARGF 1
 672 float
 673 cargf (float complex z)
 674 {
 675   return atan2f (IMAGPART (z), REALPART (z));
 676 }
 677 #endif
 678
 679 #if !defined(HAVE_CARG)
 680 #define HAVE_CARG 1
 681 double
 682 carg (double complex z)
 683 {
 684   return atan2 (IMAGPART (z), REALPART (z));
 685 }
 686 #endif
 687
 688 #if !defined(HAVE_CARGL) && defined(HAVE_ATAN2L)
 689 #define HAVE_CARGL 1
 690 long double
 691 cargl (long double complex z)
 692 {
 693   return atan2l (IMAGPART (z), REALPART (z));
 694 }
 695 #endif
 696
 697
 698 /* exp(z) = exp(a)*(cos(b) + i sin(b))  */
 699 #if !defined(HAVE_CEXPF)
 700 #define HAVE_CEXPF 1
 701 float complex
 702 cexpf (float complex z)
 703 {
 704   float a, b;
 705   float complex v;
 706
 707   a = REALPART (z);
 708   b = IMAGPART (z);
 709   COMPLEX_ASSIGN (v, cosf (b), sinf (b));
 710   return expf (a) * v;
 711 }
 712 #endif
 713
 714 #if !defined(HAVE_CEXP)
 715 #define HAVE_CEXP 1
 716 double complex
 717 cexp (double complex z)
 718 {
 719   double a, b;
 720   double complex v;
 721
 722   a = REALPART (z);
 723   b = IMAGPART (z);
 724   COMPLEX_ASSIGN (v, cos (b), sin (b));
 725   return exp (a) * v;
 726 }
 727 #endif
 728
 729 #if !defined(HAVE_CEXPL) && defined(HAVE_COSL) && defined(HAVE_SINL) && defined(EXPL)
 730 #define HAVE_CEXPL 1
 731 long double complex
 732 cexpl (long double complex z)
 733 {
 734   long double a, b;
 735   long double complex v;
 736
 737   a = REALPART (z);
 738   b = IMAGPART (z);
 739   COMPLEX_ASSIGN (v, cosl (b), sinl (b));
 740   return expl (a) * v;
 741 }
 742 #endif
 743
 744
 745 /* log(z) = log (cabs(z)) + i*carg(z)  */
 746 #if !defined(HAVE_CLOGF)
 747 #define HAVE_CLOGF 1
 748 float complex
 749 clogf (float complex z)
 750 {
 751   float complex v;
 752
 753   COMPLEX_ASSIGN (v, logf (cabsf (z)), cargf (z));
 754   return v;
 755 }
 756 #endif
 757
 758 #if !defined(HAVE_CLOG)
 759 #define HAVE_CLOG 1
 760 double complex
 761 clog (double complex z)
 762 {
 763   double complex v;
 764
 765   COMPLEX_ASSIGN (v, log (cabs (z)), carg (z));
 766   return v;
 767 }
 768 #endif
 769
 770 #if !defined(HAVE_CLOGL) && defined(HAVE_LOGL) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
 771 #define HAVE_CLOGL 1
 772 long double complex
 773 clogl (long double complex z)
 774 {
 775   long double complex v;
 776
 777   COMPLEX_ASSIGN (v, logl (cabsl (z)), cargl (z));
 778   return v;
 779 }
 780 #endif
 781
 782
 783 /* log10(z) = log10 (cabs(z)) + i*carg(z)  */
 784 #if !defined(HAVE_CLOG10F)
 785 #define HAVE_CLOG10F 1
 786 float complex
 787 clog10f (float complex z)
 788 {
 789   float complex v;
 790
 791   COMPLEX_ASSIGN (v, log10f (cabsf (z)), cargf (z));
 792   return v;
 793 }
 794 #endif
 795
 796 #if !defined(HAVE_CLOG10)
 797 #define HAVE_CLOG10 1
 798 double complex
 799 clog10 (double complex z)
 800 {
 801   double complex v;
 802
 803   COMPLEX_ASSIGN (v, log10 (cabs (z)), carg (z));
 804   return v;
 805 }
 806 #endif
 807
 808 #if !defined(HAVE_CLOG10L) && defined(HAVE_LOG10L) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
 809 #define HAVE_CLOG10L 1
 810 long double complex
 811 clog10l (long double complex z)
 812 {
 813   long double complex v;
 814
 815   COMPLEX_ASSIGN (v, log10l (cabsl (z)), cargl (z));
 816   return v;
 817 }
 818 #endif
 819
 820
 821 /* pow(base, power) = cexp (power * clog (base))  */
 822 #if !defined(HAVE_CPOWF)
 823 #define HAVE_CPOWF 1
 824 float complex
 825 cpowf (float complex base, float complex power)
 826 {
 827   return cexpf (power * clogf (base));
 828 }
 829 #endif
 830
 831 #if !defined(HAVE_CPOW)
 832 #define HAVE_CPOW 1
 833 double complex
 834 cpow (double complex base, double complex power)
 835 {
 836   return cexp (power * clog (base));
 837 }
 838 #endif
 839
 840 #if !defined(HAVE_CPOWL) && defined(HAVE_CEXPL) && defined(HAVE_CLOGL)
 841 #define HAVE_CPOWL 1
 842 long double complex
 843 cpowl (long double complex base, long double complex power)
 844 {
 845   return cexpl (power * clogl (base));
 846 }
 847 #endif
 848
 849
 850 /* sqrt(z).  Algorithm pulled from glibc.  */
 851 #if !defined(HAVE_CSQRTF)
 852 #define HAVE_CSQRTF 1
 853 float complex
 854 csqrtf (float complex z)
 855 {
 856   float re, im;
 857   float complex v;
 858
 859   re = REALPART (z);
 860   im = IMAGPART (z);
 861   if (im == 0)
 862     {
 863       if (re < 0)
 864         {
 865           COMPLEX_ASSIGN (v, 0, copysignf (sqrtf (-re), im));
 866         }
 867       else
 868         {
 869           COMPLEX_ASSIGN (v, fabsf (sqrtf (re)), copysignf (0, im));
 870         }
 871     }
 872   else if (re == 0)
 873     {
 874       float r;
 875
 876       r = sqrtf (0.5 * fabsf (im));
 877
 878       COMPLEX_ASSIGN (v, r, copysignf (r, im));
 879     }
 880   else
 881     {
 882       float d, r, s;
 883
 884       d = hypotf (re, im);
 885       /* Use the identity   2  Re res  Im res = Im x
 886          to avoid cancellation error in  d +/- Re x.  */
 887       if (re > 0)
 888         {
 889           r = sqrtf (0.5 * d + 0.5 * re);
 890           s = (0.5 * im) / r;
 891         }
 892       else
 893         {
 894           s = sqrtf (0.5 * d - 0.5 * re);
 895           r = fabsf ((0.5 * im) / s);
 896         }
 897
 898       COMPLEX_ASSIGN (v, r, copysignf (s, im));
 899     }
 900   return v;
 901 }
 902 #endif
 903
 904 #if !defined(HAVE_CSQRT)
 905 #define HAVE_CSQRT 1
 906 double complex
 907 csqrt (double complex z)
 908 {
 909   double re, im;
 910   double complex v;
 911
 912   re = REALPART (z);
 913   im = IMAGPART (z);
 914   if (im == 0)
 915     {
 916       if (re < 0)
 917         {
 918           COMPLEX_ASSIGN (v, 0, copysign (sqrt (-re), im));
 919         }
 920       else
 921         {
 922           COMPLEX_ASSIGN (v, fabs (sqrt (re)), copysign (0, im));
 923         }
 924     }
 925   else if (re == 0)
 926     {
 927       double r;
 928
 929       r = sqrt (0.5 * fabs (im));
 930
 931       COMPLEX_ASSIGN (v, r, copysign (r, im));
 932     }
 933   else
 934     {
 935       double d, r, s;
 936
 937       d = hypot (re, im);
 938       /* Use the identity   2  Re res  Im res = Im x
 939          to avoid cancellation error in  d +/- Re x.  */
 940       if (re > 0)
 941         {
 942           r = sqrt (0.5 * d + 0.5 * re);
 943           s = (0.5 * im) / r;
 944         }
 945       else
 946         {
 947           s = sqrt (0.5 * d - 0.5 * re);
 948           r = fabs ((0.5 * im) / s);
 949         }
 950
 951       COMPLEX_ASSIGN (v, r, copysign (s, im));
 952     }
 953   return v;
 954 }
 955 #endif
 956
 957 #if !defined(HAVE_CSQRTL) && defined(HAVE_COPYSIGNL) && defined(HAVE_SQRTL) && defined(HAVE_FABSL) && defined(HAVE_HYPOTL)
 958 #define HAVE_CSQRTL 1
 959 long double complex
 960 csqrtl (long double complex z)
 961 {
 962   long double re, im;
 963   long double complex v;
 964
 965   re = REALPART (z);
 966   im = IMAGPART (z);
 967   if (im == 0)
 968     {
 969       if (re < 0)
 970         {
 971           COMPLEX_ASSIGN (v, 0, copysignl (sqrtl (-re), im));
 972         }
 973       else
 974         {
 975           COMPLEX_ASSIGN (v, fabsl (sqrtl (re)), copysignl (0, im));
 976         }
 977     }
 978   else if (re == 0)
 979     {
 980       long double r;
 981
 982       r = sqrtl (0.5 * fabsl (im));
 983
 984       COMPLEX_ASSIGN (v, copysignl (r, im), r);
 985     }
 986   else
 987     {
 988       long double d, r, s;
 989
 990       d = hypotl (re, im);
 991       /* Use the identity   2  Re res  Im res = Im x
 992          to avoid cancellation error in  d +/- Re x.  */
 993       if (re > 0)
 994         {
 995           r = sqrtl (0.5 * d + 0.5 * re);
 996           s = (0.5 * im) / r;
 997         }
 998       else
 999         {
1000           s = sqrtl (0.5 * d - 0.5 * re);
1001           r = fabsl ((0.5 * im) / s);
1002         }
1003
1004       COMPLEX_ASSIGN (v, r, copysignl (s, im));
1005     }
1006   return v;
1007 }
1008 #endif
1009
1010
1011 /* sinh(a + i b) = sinh(a) cos(b) + i cosh(a) sin(b)  */
1012 #if !defined(HAVE_CSINHF)
1013 #define HAVE_CSINHF 1
1014 float complex
1015 csinhf (float complex a)
1016 {
1017   float r, i;
1018   float complex v;
1019
1020   r = REALPART (a);
1021   i = IMAGPART (a);
1022   COMPLEX_ASSIGN (v, sinhf (r) * cosf (i), coshf (r) * sinf (i));
1023   return v;
1024 }
1025 #endif
1026
1027 #if !defined(HAVE_CSINH)
1028 #define HAVE_CSINH 1
1029 double complex
1030 csinh (double complex a)
1031 {
1032   double r, i;
1033   double complex v;
1034
1035   r = REALPART (a);
1036   i = IMAGPART (a);
1037   COMPLEX_ASSIGN (v, sinh (r) * cos (i), cosh (r) * sin (i));
1038   return v;
1039 }
1040 #endif
1041
1042 #if !defined(HAVE_CSINHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1043 #define HAVE_CSINHL 1
1044 long double complex
1045 csinhl (long double complex a)
1046 {
1047   long double r, i;
1048   long double complex v;
1049
1050   r = REALPART (a);
1051   i = IMAGPART (a);
1052   COMPLEX_ASSIGN (v, sinhl (r) * cosl (i), coshl (r) * sinl (i));
1053   return v;
1054 }
1055 #endif
1056
1057
1058 /* cosh(a + i b) = cosh(a) cos(b) - i sinh(a) sin(b)  */
1059 #if !defined(HAVE_CCOSHF)
1060 #define HAVE_CCOSHF 1
1061 float complex
1062 ccoshf (float complex a)
1063 {
1064   float r, i;
1065   float complex v;
1066
1067   r = REALPART (a);
1068   i = IMAGPART (a);
1069   COMPLEX_ASSIGN (v, coshf (r) * cosf (i), - (sinhf (r) * sinf (i)));
1070   return v;
1071 }
1072 #endif
1073
1074 #if !defined(HAVE_CCOSH)
1075 #define HAVE_CCOSH 1
1076 double complex
1077 ccosh (double complex a)
1078 {
1079   double r, i;
1080   double complex v;
1081
1082   r = REALPART (a);
1083   i = IMAGPART (a);
1084   COMPLEX_ASSIGN (v, cosh (r) * cos (i), - (sinh (r) * sin (i)));
1085   return v;
1086 }
1087 #endif
1088
1089 #if !defined(HAVE_CCOSHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1090 #define HAVE_CCOSHL 1
1091 long double complex
1092 ccoshl (long double complex a)
1093 {
1094   long double r, i;
1095   long double complex v;
1096
1097   r = REALPART (a);
1098   i = IMAGPART (a);
1099   COMPLEX_ASSIGN (v, coshl (r) * cosl (i), - (sinhl (r) * sinl (i)));
1100   return v;
1101 }
1102 #endif
1103
1104
1105 /* tanh(a + i b) = (tanh(a) + i tan(b)) / (1 - i tanh(a) tan(b))  */
1106 #if !defined(HAVE_CTANHF)
1107 #define HAVE_CTANHF 1
1108 float complex
1109 ctanhf (float complex a)
1110 {
1111   float rt, it;
1112   float complex n, d;
1113
1114   rt = tanhf (REALPART (a));
1115   it = tanf (IMAGPART (a));
1116   COMPLEX_ASSIGN (n, rt, it);
1117   COMPLEX_ASSIGN (d, 1, - (rt * it));
1118
1119   return n / d;
1120 }
1121 #endif
1122
1123 #if !defined(HAVE_CTANH)
1124 #define HAVE_CTANH 1
1125 double complex
1126 ctanh (double complex a)
1127 {
1128   double rt, it;
1129   double complex n, d;
1130
1131   rt = tanh (REALPART (a));
1132   it = tan (IMAGPART (a));
1133   COMPLEX_ASSIGN (n, rt, it);
1134   COMPLEX_ASSIGN (d, 1, - (rt * it));
1135
1136   return n / d;
1137 }
1138 #endif
1139
1140 #if !defined(HAVE_CTANHL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1141 #define HAVE_CTANHL 1
1142 long double complex
1143 ctanhl (long double complex a)
1144 {
1145   long double rt, it;
1146   long double complex n, d;
1147
1148   rt = tanhl (REALPART (a));
1149   it = tanl (IMAGPART (a));
1150   COMPLEX_ASSIGN (n, rt, it);
1151   COMPLEX_ASSIGN (d, 1, - (rt * it));
1152
1153   return n / d;
1154 }
1155 #endif
1156
1157
1158 /* sin(a + i b) = sin(a) cosh(b) + i cos(a) sinh(b)  */
1159 #if !defined(HAVE_CSINF)
1160 #define HAVE_CSINF 1
1161 float complex
1162 csinf (float complex a)
1163 {
1164   float r, i;
1165   float complex v;
1166
1167   r = REALPART (a);
1168   i = IMAGPART (a);
1169   COMPLEX_ASSIGN (v, sinf (r) * coshf (i), cosf (r) * sinhf (i));
1170   return v;
1171 }
1172 #endif
1173
1174 #if !defined(HAVE_CSIN)
1175 #define HAVE_CSIN 1
1176 double complex
1177 csin (double complex a)
1178 {
1179   double r, i;
1180   double complex v;
1181
1182   r = REALPART (a);
1183   i = IMAGPART (a);
1184   COMPLEX_ASSIGN (v, sin (r) * cosh (i), cos (r) * sinh (i));
1185   return v;
1186 }
1187 #endif
1188
1189 #if !defined(HAVE_CSINL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1190 #define HAVE_CSINL 1
1191 long double complex
1192 csinl (long double complex a)
1193 {
1194   long double r, i;
1195   long double complex v;
1196
1197   r = REALPART (a);
1198   i = IMAGPART (a);
1199   COMPLEX_ASSIGN (v, sinl (r) * coshl (i), cosl (r) * sinhl (i));
1200   return v;
1201 }
1202 #endif
1203
1204
1205 /* cos(a + i b) = cos(a) cosh(b) - i sin(a) sinh(b)  */
1206 #if !defined(HAVE_CCOSF)
1207 #define HAVE_CCOSF 1
1208 float complex
1209 ccosf (float complex a)
1210 {
1211   float r, i;
1212   float complex v;
1213
1214   r = REALPART (a);
1215   i = IMAGPART (a);
1216   COMPLEX_ASSIGN (v, cosf (r) * coshf (i), - (sinf (r) * sinhf (i)));
1217   return v;
1218 }
1219 #endif
1220
1221 #if !defined(HAVE_CCOS)
1222 #define HAVE_CCOS 1
1223 double complex
1224 ccos (double complex a)
1225 {
1226   double r, i;
1227   double complex v;
1228
1229   r = REALPART (a);
1230   i = IMAGPART (a);
1231   COMPLEX_ASSIGN (v, cos (r) * cosh (i), - (sin (r) * sinh (i)));
1232   return v;
1233 }
1234 #endif
1235
1236 #if !defined(HAVE_CCOSL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1237 #define HAVE_CCOSL 1
1238 long double complex
1239 ccosl (long double complex a)
1240 {
1241   long double r, i;
1242   long double complex v;
1243
1244   r = REALPART (a);
1245   i = IMAGPART (a);
1246   COMPLEX_ASSIGN (v, cosl (r) * coshl (i), - (sinl (r) * sinhl (i)));
1247   return v;
1248 }
1249 #endif
1250
1251
1252 /* tan(a + i b) = (tan(a) + i tanh(b)) / (1 - i tan(a) tanh(b))  */
1253 #if !defined(HAVE_CTANF)
1254 #define HAVE_CTANF 1
1255 float complex
1256 ctanf (float complex a)
1257 {
1258   float rt, it;
1259   float complex n, d;
1260
1261   rt = tanf (REALPART (a));
1262   it = tanhf (IMAGPART (a));
1263   COMPLEX_ASSIGN (n, rt, it);
1264   COMPLEX_ASSIGN (d, 1, - (rt * it));
1265
1266   return n / d;
1267 }
1268 #endif
1269
1270 #if !defined(HAVE_CTAN)
1271 #define HAVE_CTAN 1
1272 double complex
1273 ctan (double complex a)
1274 {
1275   double rt, it;
1276   double complex n, d;
1277
1278   rt = tan (REALPART (a));
1279   it = tanh (IMAGPART (a));
1280   COMPLEX_ASSIGN (n, rt, it);
1281   COMPLEX_ASSIGN (d, 1, - (rt * it));
1282
1283   return n / d;
1284 }
1285 #endif
1286
1287 #if !defined(HAVE_CTANL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1288 #define HAVE_CTANL 1
1289 long double complex
1290 ctanl (long double complex a)
1291 {
1292   long double rt, it;
1293   long double complex n, d;
1294
1295   rt = tanl (REALPART (a));
1296   it = tanhl (IMAGPART (a));
1297   COMPLEX_ASSIGN (n, rt, it);
1298   COMPLEX_ASSIGN (d, 1, - (rt * it));
1299
1300   return n / d;
1301 }
1302 #endif
1303