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