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