OSDN Git Service

96b5ef88851a39aaa6d6bc837a2ad32a46c574be
[pf3gnuchains/gcc-fork.git] / 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