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_FREXPF
 298 #define HAVE_FREXPF 1
 299 float
 300 frexpf(float x, int *exp)
 301 {
 302   return (float) frexp(x, exp);
 303 }
 304 #endif
 305
 306 #ifndef HAVE_HYPOTF
 307 #define HAVE_HYPOTF 1
 308 float
 309 hypotf(float x, float y)
 310 {
 311   return (float) hypot(x, y);
 312 }
 313 #endif
 314
 315 #ifndef HAVE_LOGF
 316 #define HAVE_LOGF 1
 317 float
 318 logf(float x)
 319 {
 320   return (float) log(x);
 321 }
 322 #endif
 323
 324 #ifndef HAVE_LOG10F
 325 #define HAVE_LOG10F 1
 326 float
 327 log10f(float x)
 328 {
 329   return (float) log10(x);
 330 }
 331 #endif
 332
 333 #ifndef HAVE_SCALBN
 334 #define HAVE_SCALBN 1
 335 double
 336 scalbn(double x, int y)
 337 {
 338   return x * pow(FLT_RADIX, y);
 339 }
 340 #endif
 341
 342 #ifndef HAVE_SCALBNF
 343 #define HAVE_SCALBNF 1
 344 float
 345 scalbnf(float x, int y)
 346 {
 347   return (float) scalbn(x, y);
 348 }
 349 #endif
 350
 351 #ifndef HAVE_SINF
 352 #define HAVE_SINF 1
 353 float
 354 sinf(float x)
 355 {
 356   return (float) sin(x);
 357 }
 358 #endif
 359
 360 #ifndef HAVE_SINHF
 361 #define HAVE_SINHF 1
 362 float
 363 sinhf(float x)
 364 {
 365   return (float) sinh(x);
 366 }
 367 #endif
 368
 369 #ifndef HAVE_SQRTF
 370 #define HAVE_SQRTF 1
 371 float
 372 sqrtf(float x)
 373 {
 374   return (float) sqrt(x);
 375 }
 376 #endif
 377
 378 #ifndef HAVE_TANF
 379 #define HAVE_TANF 1
 380 float
 381 tanf(float x)
 382 {
 383   return (float) tan(x);
 384 }
 385 #endif
 386
 387 #ifndef HAVE_TANHF
 388 #define HAVE_TANHF 1
 389 float
 390 tanhf(float x)
 391 {
 392   return (float) tanh(x);
 393 }
 394 #endif
 395
 396 #ifndef HAVE_TRUNC
 397 #define HAVE_TRUNC 1
 398 double
 399 trunc(double x)
 400 {
 401   if (!isfinite (x))
 402     return x;
 403
 404   if (x < 0.0)
 405     return - floor (-x);
 406   else
 407     return floor (x);
 408 }
 409 #endif
 410
 411 #ifndef HAVE_TRUNCF
 412 #define HAVE_TRUNCF 1
 413 float
 414 truncf(float x)
 415 {
 416   return (float) trunc (x);
 417 }
 418 #endif
 419
 420 #ifndef HAVE_NEXTAFTERF
 421 #define HAVE_NEXTAFTERF 1
 422 /* This is a portable implementation of nextafterf that is intended to be
 423    independent of the floating point format or its in memory representation.
 424    This implementation works correctly with denormalized values.  */
 425 float
 426 nextafterf(float x, float y)
 427 {
 428   /* This variable is marked volatile to avoid excess precision problems
 429      on some platforms, including IA-32.  */
 430   volatile float delta;
 431   float absx, denorm_min;
 432
 433   if (isnan(x) || isnan(y))
 434     return x + y;
 435   if (x == y)
 436     return x;
 437   if (!isfinite (x))
 438     return x > 0 ? __FLT_MAX__ : - __FLT_MAX__;
 439
 440   /* absx = fabsf (x);  */
 441   absx = (x < 0.0) ? -x : x;
 442
 443   /* __FLT_DENORM_MIN__ is non-zero iff the target supports denormals.  */
 444   if (__FLT_DENORM_MIN__ == 0.0f)
 445     denorm_min = __FLT_MIN__;
 446   else
 447     denorm_min = __FLT_DENORM_MIN__;
 448
 449   if (absx < __FLT_MIN__)
 450     delta = denorm_min;
 451   else
 452     {
 453       float frac;
 454       int exp;
 455
 456       /* Discard the fraction from x.  */
 457       frac = frexpf (absx, &exp);
 458       delta = scalbnf (0.5f, exp);
 459
 460       /* Scale x by the epsilon of the representation.  By rights we should
 461          have been able to combine this with scalbnf, but some targets don't
 462          get that correct with denormals.  */
 463       delta *= __FLT_EPSILON__;
 464
 465       /* If we're going to be reducing the absolute value of X, and doing so
 466          would reduce the exponent of X, then the delta to be applied is
 467          one exponent smaller.  */
 468       if (frac == 0.5f && (y < x) == (x > 0))
 469         delta *= 0.5f;
 470
 471       /* If that underflows to zero, then we're back to the minimum.  */
 472       if (delta == 0.0f)
 473         delta = denorm_min;
 474     }
 475
 476   if (y < x)
 477     delta = -delta;
 478
 479   return x + delta;
 480 }
 481 #endif
 482
 483
 484 #ifndef HAVE_POWF
 485 #define HAVE_POWF 1
 486 float
 487 powf(float x, float y)
 488 {
 489   return (float) pow(x, y);
 490 }
 491 #endif
 492
 493 /* Note that if fpclassify is not defined, then NaN is not handled */
 494
 495 /* Algorithm by Steven G. Kargl.  */
 496
 497 #ifndef HAVE_ROUND
 498 #define HAVE_ROUND 1
 499 /* Round to nearest integral value.  If the argument is halfway between two
 500    integral values then round away from zero.  */
 501
 502 double
 503 round(double x)
 504 {
 505    double t;
 506    if (!isfinite (x))
 507      return (x);
 508
 509    if (x >= 0.0)
 510     {
 511       t = ceil(x);
 512       if (t - x > 0.5)
 513         t -= 1.0;
 514       return (t);
 515     }
 516    else
 517     {
 518       t = ceil(-x);
 519       if (t + x > 0.5)
 520         t -= 1.0;
 521       return (-t);
 522     }
 523 }
 524 #endif
 525
 526 #ifndef HAVE_ROUNDF
 527 #define HAVE_ROUNDF 1
 528 /* Round to nearest integral value.  If the argument is halfway between two
 529    integral values then round away from zero.  */
 530
 531 float
 532 roundf(float x)
 533 {
 534    float t;
 535    if (!isfinite (x))
 536      return (x);
 537
 538    if (x >= 0.0)
 539     {
 540       t = ceilf(x);
 541       if (t - x > 0.5)
 542         t -= 1.0;
 543       return (t);
 544     }
 545    else
 546     {
 547       t = ceilf(-x);
 548       if (t + x > 0.5)
 549         t -= 1.0;
 550       return (-t);
 551     }
 552 }
 553 #endif
 554
 555 #ifndef HAVE_LOG10L
 556 #define HAVE_LOG10L 1
 557 /* log10 function for long double variables. The version provided here
 558    reduces the argument until it fits into a double, then use log10.  */
 559 long double
 560 log10l(long double x)
 561 {
 562 #if LDBL_MAX_EXP > DBL_MAX_EXP
 563   if (x > DBL_MAX)
 564     {
 565       double val;
 566       int p2_result = 0;
 567       if (x > 0x1p16383L) { p2_result += 16383; x /= 0x1p16383L; }
 568       if (x > 0x1p8191L) { p2_result += 8191; x /= 0x1p8191L; }
 569       if (x > 0x1p4095L) { p2_result += 4095; x /= 0x1p4095L; }
 570       if (x > 0x1p2047L) { p2_result += 2047; x /= 0x1p2047L; }
 571       if (x > 0x1p1023L) { p2_result += 1023; x /= 0x1p1023L; }
 572       val = log10 ((double) x);
 573       return (val + p2_result * .30102999566398119521373889472449302L);
 574     }
 575 #endif
 576 #if LDBL_MIN_EXP < DBL_MIN_EXP
 577   if (x < DBL_MIN)
 578     {
 579       double val;
 580       int p2_result = 0;
 581       if (x < 0x1p-16380L) { p2_result += 16380; x /= 0x1p-16380L; }
 582       if (x < 0x1p-8189L) { p2_result += 8189; x /= 0x1p-8189L; }
 583       if (x < 0x1p-4093L) { p2_result += 4093; x /= 0x1p-4093L; }
 584       if (x < 0x1p-2045L) { p2_result += 2045; x /= 0x1p-2045L; }
 585       if (x < 0x1p-1021L) { p2_result += 1021; x /= 0x1p-1021L; }
 586       val = fabs(log10 ((double) x));
 587       return (- val - p2_result * .30102999566398119521373889472449302L);
 588     }
 589 #endif
 590     return log10 (x);
 591 }
 592 #endif
 593
 594
 595 #if !defined(HAVE_CABSF)
 596 #define HAVE_CABSF 1
 597 float
 598 cabsf (float complex z)
 599 {
 600   return hypotf (REALPART (z), IMAGPART (z));
 601 }
 602 #endif
 603
 604 #if !defined(HAVE_CABS)
 605 #define HAVE_CABS 1
 606 double
 607 cabs (double complex z)
 608 {
 609   return hypot (REALPART (z), IMAGPART (z));
 610 }
 611 #endif
 612
 613 #if !defined(HAVE_CABSL) && defined(HAVE_HYPOTL)
 614 #define HAVE_CABSL 1
 615 long double
 616 cabsl (long double complex z)
 617 {
 618   return hypotl (REALPART (z), IMAGPART (z));
 619 }
 620 #endif
 621
 622
 623 #if !defined(HAVE_CARGF)
 624 #define HAVE_CARGF 1
 625 float
 626 cargf (float complex z)
 627 {
 628   return atan2f (IMAGPART (z), REALPART (z));
 629 }
 630 #endif
 631
 632 #if !defined(HAVE_CARG)
 633 #define HAVE_CARG 1
 634 double
 635 carg (double complex z)
 636 {
 637   return atan2 (IMAGPART (z), REALPART (z));
 638 }
 639 #endif
 640
 641 #if !defined(HAVE_CARGL) && defined(HAVE_ATAN2L)
 642 #define HAVE_CARGL 1
 643 long double
 644 cargl (long double complex z)
 645 {
 646   return atan2l (IMAGPART (z), REALPART (z));
 647 }
 648 #endif
 649
 650
 651 /* exp(z) = exp(a)*(cos(b) + i sin(b))  */
 652 #if !defined(HAVE_CEXPF)
 653 #define HAVE_CEXPF 1
 654 float complex
 655 cexpf (float complex z)
 656 {
 657   float a, b;
 658   float complex v;
 659
 660   a = REALPART (z);
 661   b = IMAGPART (z);
 662   COMPLEX_ASSIGN (v, cosf (b), sinf (b));
 663   return expf (a) * v;
 664 }
 665 #endif
 666
 667 #if !defined(HAVE_CEXP)
 668 #define HAVE_CEXP 1
 669 double complex
 670 cexp (double complex z)
 671 {
 672   double a, b;
 673   double complex v;
 674
 675   a = REALPART (z);
 676   b = IMAGPART (z);
 677   COMPLEX_ASSIGN (v, cos (b), sin (b));
 678   return exp (a) * v;
 679 }
 680 #endif
 681
 682 #if !defined(HAVE_CEXPL) && defined(HAVE_COSL) && defined(HAVE_SINL) && defined(EXPL)
 683 #define HAVE_CEXPL 1
 684 long double complex
 685 cexpl (long double complex z)
 686 {
 687   long double a, b;
 688   long double complex v;
 689
 690   a = REALPART (z);
 691   b = IMAGPART (z);
 692   COMPLEX_ASSIGN (v, cosl (b), sinl (b));
 693   return expl (a) * v;
 694 }
 695 #endif
 696
 697
 698 /* log(z) = log (cabs(z)) + i*carg(z)  */
 699 #if !defined(HAVE_CLOGF)
 700 #define HAVE_CLOGF 1
 701 float complex
 702 clogf (float complex z)
 703 {
 704   float complex v;
 705
 706   COMPLEX_ASSIGN (v, logf (cabsf (z)), cargf (z));
 707   return v;
 708 }
 709 #endif
 710
 711 #if !defined(HAVE_CLOG)
 712 #define HAVE_CLOG 1
 713 double complex
 714 clog (double complex z)
 715 {
 716   double complex v;
 717
 718   COMPLEX_ASSIGN (v, log (cabs (z)), carg (z));
 719   return v;
 720 }
 721 #endif
 722
 723 #if !defined(HAVE_CLOGL) && defined(HAVE_LOGL) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
 724 #define HAVE_CLOGL 1
 725 long double complex
 726 clogl (long double complex z)
 727 {
 728   long double complex v;
 729
 730   COMPLEX_ASSIGN (v, logl (cabsl (z)), cargl (z));
 731   return v;
 732 }
 733 #endif
 734
 735
 736 /* log10(z) = log10 (cabs(z)) + i*carg(z)  */
 737 #if !defined(HAVE_CLOG10F)
 738 #define HAVE_CLOG10F 1
 739 float complex
 740 clog10f (float complex z)
 741 {
 742   float complex v;
 743
 744   COMPLEX_ASSIGN (v, log10f (cabsf (z)), cargf (z));
 745   return v;
 746 }
 747 #endif
 748
 749 #if !defined(HAVE_CLOG10)
 750 #define HAVE_CLOG10 1
 751 double complex
 752 clog10 (double complex z)
 753 {
 754   double complex v;
 755
 756   COMPLEX_ASSIGN (v, log10 (cabs (z)), carg (z));
 757   return v;
 758 }
 759 #endif
 760
 761 #if !defined(HAVE_CLOG10L) && defined(HAVE_LOG10L) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
 762 #define HAVE_CLOG10L 1
 763 long double complex
 764 clog10l (long double complex z)
 765 {
 766   long double complex v;
 767
 768   COMPLEX_ASSIGN (v, log10l (cabsl (z)), cargl (z));
 769   return v;
 770 }
 771 #endif
 772
 773
 774 /* pow(base, power) = cexp (power * clog (base))  */
 775 #if !defined(HAVE_CPOWF)
 776 #define HAVE_CPOWF 1
 777 float complex
 778 cpowf (float complex base, float complex power)
 779 {
 780   return cexpf (power * clogf (base));
 781 }
 782 #endif
 783
 784 #if !defined(HAVE_CPOW)
 785 #define HAVE_CPOW 1
 786 double complex
 787 cpow (double complex base, double complex power)
 788 {
 789   return cexp (power * clog (base));
 790 }
 791 #endif
 792
 793 #if !defined(HAVE_CPOWL) && defined(HAVE_CEXPL) && defined(HAVE_CLOGL)
 794 #define HAVE_CPOWL 1
 795 long double complex
 796 cpowl (long double complex base, long double complex power)
 797 {
 798   return cexpl (power * clogl (base));
 799 }
 800 #endif
 801
 802
 803 /* sqrt(z).  Algorithm pulled from glibc.  */
 804 #if !defined(HAVE_CSQRTF)
 805 #define HAVE_CSQRTF 1
 806 float complex
 807 csqrtf (float complex z)
 808 {
 809   float re, im;
 810   float complex v;
 811
 812   re = REALPART (z);
 813   im = IMAGPART (z);
 814   if (im == 0)
 815     {
 816       if (re < 0)
 817         {
 818           COMPLEX_ASSIGN (v, 0, copysignf (sqrtf (-re), im));
 819         }
 820       else
 821         {
 822           COMPLEX_ASSIGN (v, fabsf (sqrtf (re)), copysignf (0, im));
 823         }
 824     }
 825   else if (re == 0)
 826     {
 827       float r;
 828
 829       r = sqrtf (0.5 * fabsf (im));
 830
 831       COMPLEX_ASSIGN (v, r, copysignf (r, im));
 832     }
 833   else
 834     {
 835       float d, r, s;
 836
 837       d = hypotf (re, im);
 838       /* Use the identity   2  Re res  Im res = Im x
 839          to avoid cancellation error in  d +/- Re x.  */
 840       if (re > 0)
 841         {
 842           r = sqrtf (0.5 * d + 0.5 * re);
 843           s = (0.5 * im) / r;
 844         }
 845       else
 846         {
 847           s = sqrtf (0.5 * d - 0.5 * re);
 848           r = fabsf ((0.5 * im) / s);
 849         }
 850
 851       COMPLEX_ASSIGN (v, r, copysignf (s, im));
 852     }
 853   return v;
 854 }
 855 #endif
 856
 857 #if !defined(HAVE_CSQRT)
 858 #define HAVE_CSQRT 1
 859 double complex
 860 csqrt (double complex z)
 861 {
 862   double re, im;
 863   double complex v;
 864
 865   re = REALPART (z);
 866   im = IMAGPART (z);
 867   if (im == 0)
 868     {
 869       if (re < 0)
 870         {
 871           COMPLEX_ASSIGN (v, 0, copysign (sqrt (-re), im));
 872         }
 873       else
 874         {
 875           COMPLEX_ASSIGN (v, fabs (sqrt (re)), copysign (0, im));
 876         }
 877     }
 878   else if (re == 0)
 879     {
 880       double r;
 881
 882       r = sqrt (0.5 * fabs (im));
 883
 884       COMPLEX_ASSIGN (v, r, copysign (r, im));
 885     }
 886   else
 887     {
 888       double d, r, s;
 889
 890       d = hypot (re, im);
 891       /* Use the identity   2  Re res  Im res = Im x
 892          to avoid cancellation error in  d +/- Re x.  */
 893       if (re > 0)
 894         {
 895           r = sqrt (0.5 * d + 0.5 * re);
 896           s = (0.5 * im) / r;
 897         }
 898       else
 899         {
 900           s = sqrt (0.5 * d - 0.5 * re);
 901           r = fabs ((0.5 * im) / s);
 902         }
 903
 904       COMPLEX_ASSIGN (v, r, copysign (s, im));
 905     }
 906   return v;
 907 }
 908 #endif
 909
 910 #if !defined(HAVE_CSQRTL) && defined(HAVE_COPYSIGNL) && defined(HAVE_SQRTL) && defined(HAVE_FABSL) && defined(HAVE_HYPOTL)
 911 #define HAVE_CSQRTL 1
 912 long double complex
 913 csqrtl (long double complex z)
 914 {
 915   long double re, im;
 916   long double complex v;
 917
 918   re = REALPART (z);
 919   im = IMAGPART (z);
 920   if (im == 0)
 921     {
 922       if (re < 0)
 923         {
 924           COMPLEX_ASSIGN (v, 0, copysignl (sqrtl (-re), im));
 925         }
 926       else
 927         {
 928           COMPLEX_ASSIGN (v, fabsl (sqrtl (re)), copysignl (0, im));
 929         }
 930     }
 931   else if (re == 0)
 932     {
 933       long double r;
 934
 935       r = sqrtl (0.5 * fabsl (im));
 936
 937       COMPLEX_ASSIGN (v, copysignl (r, im), r);
 938     }
 939   else
 940     {
 941       long double d, r, s;
 942
 943       d = hypotl (re, im);
 944       /* Use the identity   2  Re res  Im res = Im x
 945          to avoid cancellation error in  d +/- Re x.  */
 946       if (re > 0)
 947         {
 948           r = sqrtl (0.5 * d + 0.5 * re);
 949           s = (0.5 * im) / r;
 950         }
 951       else
 952         {
 953           s = sqrtl (0.5 * d - 0.5 * re);
 954           r = fabsl ((0.5 * im) / s);
 955         }
 956
 957       COMPLEX_ASSIGN (v, r, copysignl (s, im));
 958     }
 959   return v;
 960 }
 961 #endif
 962
 963
 964 /* sinh(a + i b) = sinh(a) cos(b) + i cosh(a) sin(b)  */
 965 #if !defined(HAVE_CSINHF)
 966 #define HAVE_CSINHF 1
 967 float complex
 968 csinhf (float complex a)
 969 {
 970   float r, i;
 971   float complex v;
 972
 973   r = REALPART (a);
 974   i = IMAGPART (a);
 975   COMPLEX_ASSIGN (v, sinhf (r) * cosf (i), coshf (r) * sinf (i));
 976   return v;
 977 }
 978 #endif
 979
 980 #if !defined(HAVE_CSINH)
 981 #define HAVE_CSINH 1
 982 double complex
 983 csinh (double complex a)
 984 {
 985   double r, i;
 986   double complex v;
 987
 988   r = REALPART (a);
 989   i = IMAGPART (a);
 990   COMPLEX_ASSIGN (v, sinh (r) * cos (i), cosh (r) * sin (i));
 991   return v;
 992 }
 993 #endif
 994
 995 #if !defined(HAVE_CSINHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
 996 #define HAVE_CSINHL 1
 997 long double complex
 998 csinhl (long double complex a)
 999 {
1000   long double r, i;
1001   long double complex v;
1002
1003   r = REALPART (a);
1004   i = IMAGPART (a);
1005   COMPLEX_ASSIGN (v, sinhl (r) * cosl (i), coshl (r) * sinl (i));
1006   return v;
1007 }
1008 #endif
1009
1010
1011 /* cosh(a + i b) = cosh(a) cos(b) - i sinh(a) sin(b)  */
1012 #if !defined(HAVE_CCOSHF)
1013 #define HAVE_CCOSHF 1
1014 float complex
1015 ccoshf (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, coshf (r) * cosf (i), - (sinhf (r) * sinf (i)));
1023   return v;
1024 }
1025 #endif
1026
1027 #if !defined(HAVE_CCOSH)
1028 #define HAVE_CCOSH 1
1029 double complex
1030 ccosh (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, cosh (r) * cos (i), - (sinh (r) * sin (i)));
1038   return v;
1039 }
1040 #endif
1041
1042 #if !defined(HAVE_CCOSHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1043 #define HAVE_CCOSHL 1
1044 long double complex
1045 ccoshl (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, coshl (r) * cosl (i), - (sinhl (r) * sinl (i)));
1053   return v;
1054 }
1055 #endif
1056
1057
1058 /* tanh(a + i b) = (tanh(a) + i tan(b)) / (1 - i tanh(a) tan(b))  */
1059 #if !defined(HAVE_CTANHF)
1060 #define HAVE_CTANHF 1
1061 float complex
1062 ctanhf (float complex a)
1063 {
1064   float rt, it;
1065   float complex n, d;
1066
1067   rt = tanhf (REALPART (a));
1068   it = tanf (IMAGPART (a));
1069   COMPLEX_ASSIGN (n, rt, it);
1070   COMPLEX_ASSIGN (d, 1, - (rt * it));
1071
1072   return n / d;
1073 }
1074 #endif
1075
1076 #if !defined(HAVE_CTANH)
1077 #define HAVE_CTANH 1
1078 double complex
1079 ctanh (double complex a)
1080 {
1081   double rt, it;
1082   double complex n, d;
1083
1084   rt = tanh (REALPART (a));
1085   it = tan (IMAGPART (a));
1086   COMPLEX_ASSIGN (n, rt, it);
1087   COMPLEX_ASSIGN (d, 1, - (rt * it));
1088
1089   return n / d;
1090 }
1091 #endif
1092
1093 #if !defined(HAVE_CTANHL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1094 #define HAVE_CTANHL 1
1095 long double complex
1096 ctanhl (long double complex a)
1097 {
1098   long double rt, it;
1099   long double complex n, d;
1100
1101   rt = tanhl (REALPART (a));
1102   it = tanl (IMAGPART (a));
1103   COMPLEX_ASSIGN (n, rt, it);
1104   COMPLEX_ASSIGN (d, 1, - (rt * it));
1105
1106   return n / d;
1107 }
1108 #endif
1109
1110
1111 /* sin(a + i b) = sin(a) cosh(b) + i cos(a) sinh(b)  */
1112 #if !defined(HAVE_CSINF)
1113 #define HAVE_CSINF 1
1114 float complex
1115 csinf (float complex a)
1116 {
1117   float r, i;
1118   float complex v;
1119
1120   r = REALPART (a);
1121   i = IMAGPART (a);
1122   COMPLEX_ASSIGN (v, sinf (r) * coshf (i), cosf (r) * sinhf (i));
1123   return v;
1124 }
1125 #endif
1126
1127 #if !defined(HAVE_CSIN)
1128 #define HAVE_CSIN 1
1129 double complex
1130 csin (double complex a)
1131 {
1132   double r, i;
1133   double complex v;
1134
1135   r = REALPART (a);
1136   i = IMAGPART (a);
1137   COMPLEX_ASSIGN (v, sin (r) * cosh (i), cos (r) * sinh (i));
1138   return v;
1139 }
1140 #endif
1141
1142 #if !defined(HAVE_CSINL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1143 #define HAVE_CSINL 1
1144 long double complex
1145 csinl (long double complex a)
1146 {
1147   long double r, i;
1148   long double complex v;
1149
1150   r = REALPART (a);
1151   i = IMAGPART (a);
1152   COMPLEX_ASSIGN (v, sinl (r) * coshl (i), cosl (r) * sinhl (i));
1153   return v;
1154 }
1155 #endif
1156
1157
1158 /* cos(a + i b) = cos(a) cosh(b) - i sin(a) sinh(b)  */
1159 #if !defined(HAVE_CCOSF)
1160 #define HAVE_CCOSF 1
1161 float complex
1162 ccosf (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, cosf (r) * coshf (i), - (sinf (r) * sinhf (i)));
1170   return v;
1171 }
1172 #endif
1173
1174 #if !defined(HAVE_CCOS)
1175 #define HAVE_CCOS 1
1176 double complex
1177 ccos (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, cos (r) * cosh (i), - (sin (r) * sinh (i)));
1185   return v;
1186 }
1187 #endif
1188
1189 #if !defined(HAVE_CCOSL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1190 #define HAVE_CCOSL 1
1191 long double complex
1192 ccosl (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, cosl (r) * coshl (i), - (sinl (r) * sinhl (i)));
1200   return v;
1201 }
1202 #endif
1203
1204
1205 /* tan(a + i b) = (tan(a) + i tanh(b)) / (1 - i tan(a) tanh(b))  */
1206 #if !defined(HAVE_CTANF)
1207 #define HAVE_CTANF 1
1208 float complex
1209 ctanf (float complex a)
1210 {
1211   float rt, it;
1212   float complex n, d;
1213
1214   rt = tanf (REALPART (a));
1215   it = tanhf (IMAGPART (a));
1216   COMPLEX_ASSIGN (n, rt, it);
1217   COMPLEX_ASSIGN (d, 1, - (rt * it));
1218
1219   return n / d;
1220 }
1221 #endif
1222
1223 #if !defined(HAVE_CTAN)
1224 #define HAVE_CTAN 1
1225 double complex
1226 ctan (double complex a)
1227 {
1228   double rt, it;
1229   double complex n, d;
1230
1231   rt = tan (REALPART (a));
1232   it = tanh (IMAGPART (a));
1233   COMPLEX_ASSIGN (n, rt, it);
1234   COMPLEX_ASSIGN (d, 1, - (rt * it));
1235
1236   return n / d;
1237 }
1238 #endif
1239
1240 #if !defined(HAVE_CTANL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1241 #define HAVE_CTANL 1
1242 long double complex
1243 ctanl (long double complex a)
1244 {
1245   long double rt, it;
1246   long double complex n, d;
1247
1248   rt = tanl (REALPART (a));
1249   it = tanhl (IMAGPART (a));
1250   COMPLEX_ASSIGN (n, rt, it);
1251   COMPLEX_ASSIGN (d, 1, - (rt * it));
1252
1253   return n / d;
1254 }
1255 #endif
1256