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