1 /* Implementation of various C99 functions
2 Copyright (C) 2004 Free Software Foundation, Inc.
4 This file is part of the GNU Fortran 95 runtime library (libgfortran).
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.
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
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.
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. */
31 #include <sys/types.h>
35 #define C99_PROTOS_H WE_DONT_WANT_PROTOS_NOW
36 #include "libgfortran.h"
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}. */
43 #if defined(__sgi__) && !defined(HAVE_COMPLEX_H)
47 #define cabs __gfc_cabs
48 #define cabsf __gfc_cabsf
49 #define cabsl __gfc_cabsl
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}. */
60 #define cabs __gfc_cabs
61 #define cabsf __gfc_cabsf
62 #define cabsl __gfc_cabsl
65 /* Prototypes to silence -Wstrict-prototypes -Wmissing-prototypes. */
67 float cabsf(float complex);
68 double cabs(double complex);
69 long double cabsl(long double complex);
71 float cargf(float complex);
72 double carg(double complex);
73 long double cargl(long double complex);
75 float complex clog10f(float complex);
76 double complex clog10(double complex);
77 long double complex clog10l(long double complex);
85 return (float) acos(x
);
94 return (float) asin(x
);
101 atan2f(float y
, float x
)
103 return (float) atan2(y
, x
);
112 return (float) atan(x
);
121 return (float) ceil(x
);
125 #ifndef HAVE_COPYSIGNF
126 #define HAVE_COPYSIGNF 1
128 copysignf(float x
, float y
)
130 return (float) copysign(x
, y
);
139 return (float) cos(x
);
148 return (float) cosh(x
);
157 return (float) exp(x
);
166 return (float) fabs(x
);
171 #define HAVE_FLOORF 1
175 return (float) floor(x
);
180 #define HAVE_FREXPF 1
182 frexpf(float x
, int *exp
)
184 return (float) frexp(x
, exp
);
189 #define HAVE_HYPOTF 1
191 hypotf(float x
, float y
)
193 return (float) hypot(x
, y
);
202 return (float) log(x
);
207 #define HAVE_LOG10F 1
211 return (float) log10(x
);
216 #define HAVE_SCALBN 1
218 scalbn(double x
, int y
)
220 return x
* pow(FLT_RADIX
, y
);
225 #define HAVE_SCALBNF 1
227 scalbnf(float x
, int y
)
229 return (float) scalbn(x
, y
);
238 return (float) sin(x
);
247 return (float) sinh(x
);
256 return (float) sqrt(x
);
265 return (float) tan(x
);
274 return (float) tanh(x
);
294 #define HAVE_TRUNCF 1
298 return (float) trunc (x
);
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. */
308 nextafterf(float x
, float y
)
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
;
315 if (isnan(x
) || isnan(y
))
320 return x
> 0 ? __FLT_MAX__
: - __FLT_MAX__
;
322 /* absx = fabsf (x); */
323 absx
= (x
< 0.0) ? -x
: x
;
325 /* __FLT_DENORM_MIN__ is non-zero iff the target supports denormals. */
326 if (__FLT_DENORM_MIN__
== 0.0f
)
327 denorm_min
= __FLT_MIN__
;
329 denorm_min
= __FLT_DENORM_MIN__
;
331 if (absx
< __FLT_MIN__
)
338 /* Discard the fraction from x. */
339 frac
= frexpf (absx
, &exp
);
340 delta
= scalbnf (0.5f
, exp
);
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__
;
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))
353 /* If that underflows to zero, then we're back to the minimum. */
369 powf(float x
, float y
)
371 return (float) pow(x
, y
);
375 /* Note that if fpclassify is not defined, then NaN is not handled */
377 /* Algorithm by Steven G. Kargl. */
381 /* Round to nearest integral value. If the argument is halfway between two
382 integral values then round away from zero. */
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. */
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. */
442 log10l(long double x
)
444 #if LDBL_MAX_EXP > DBL_MAX_EXP
449 if (x
> 0x1p
16383L) { p2_result
+= 16383; x
/= 0x1p
16383L; }
450 if (x
> 0x1p
8191L) { p2_result
+= 8191; x
/= 0x1p
8191L; }
451 if (x
> 0x1p
4095L) { p2_result
+= 4095; x
/= 0x1p
4095L; }
452 if (x
> 0x1p
2047L) { p2_result
+= 2047; x
/= 0x1p
2047L; }
453 if (x
> 0x1p
1023L) { p2_result
+= 1023; x
/= 0x1p
1023L; }
454 val
= log10 ((double) x
);
455 return (val
+ p2_result
* .30102999566398119521373889472449302L);
458 #if LDBL_MIN_EXP < DBL_MIN_EXP
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);
477 #if !defined(HAVE_CABSF)
480 cabsf (float complex z
)
482 return hypotf (REALPART (z
), IMAGPART (z
));
486 #if !defined(HAVE_CABS)
489 cabs (double complex z
)
491 return hypot (REALPART (z
), IMAGPART (z
));
495 #if !defined(HAVE_CABSL) && defined(HAVE_HYPOTL)
498 cabsl (long double complex z
)
500 return hypotl (REALPART (z
), IMAGPART (z
));
505 #if !defined(HAVE_CARGF)
508 cargf (float complex z
)
510 return atan2f (IMAGPART (z
), REALPART (z
));
514 #if !defined(HAVE_CARG)
517 carg (double complex z
)
519 return atan2 (IMAGPART (z
), REALPART (z
));
523 #if !defined(HAVE_CARGL) && defined(HAVE_ATAN2L)
526 cargl (long double complex z
)
528 return atan2l (IMAGPART (z
), REALPART (z
));
533 /* exp(z) = exp(a)*(cos(b) + i sin(b)) */
534 #if !defined(HAVE_CEXPF)
537 cexpf (float complex z
)
544 COMPLEX_ASSIGN (v
, cosf (b
), sinf (b
));
549 #if !defined(HAVE_CEXP)
552 cexp (double complex z
)
559 COMPLEX_ASSIGN (v
, cos (b
), sin (b
));
564 #if !defined(HAVE_CEXPL) && defined(HAVE_COSL) && defined(HAVE_SINL) && defined(EXPL)
567 cexpl (long double complex z
)
570 long double complex v
;
574 COMPLEX_ASSIGN (v
, cosl (b
), sinl (b
));
580 /* log(z) = log (cabs(z)) + i*carg(z) */
581 #if !defined(HAVE_CLOGF)
584 clogf (float complex z
)
588 COMPLEX_ASSIGN (v
, logf (cabsf (z
)), cargf (z
));
593 #if !defined(HAVE_CLOG)
596 clog (double complex z
)
600 COMPLEX_ASSIGN (v
, log (cabs (z
)), carg (z
));
605 #if !defined(HAVE_CLOGL) && defined(HAVE_LOGL) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
608 clogl (long double complex z
)
610 long double complex v
;
612 COMPLEX_ASSIGN (v
, logl (cabsl (z
)), cargl (z
));
618 /* log10(z) = log10 (cabs(z)) + i*carg(z) */
619 #if !defined(HAVE_CLOG10F)
620 #define HAVE_CLOG10F 1
622 clog10f (float complex z
)
626 COMPLEX_ASSIGN (v
, log10f (cabsf (z
)), cargf (z
));
631 #if !defined(HAVE_CLOG10)
632 #define HAVE_CLOG10 1
634 clog10 (double complex z
)
638 COMPLEX_ASSIGN (v
, log10 (cabs (z
)), carg (z
));
643 #if !defined(HAVE_CLOG10L) && defined(HAVE_LOG10L) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
644 #define HAVE_CLOG10L 1
646 clog10l (long double complex z
)
648 long double complex v
;
650 COMPLEX_ASSIGN (v
, log10l (cabsl (z
)), cargl (z
));
656 /* pow(base, power) = cexp (power * clog (base)) */
657 #if !defined(HAVE_CPOWF)
660 cpowf (float complex base
, float complex power
)
662 return cexpf (power
* clogf (base
));
666 #if !defined(HAVE_CPOW)
669 cpow (double complex base
, double complex power
)
671 return cexp (power
* clog (base
));
675 #if !defined(HAVE_CPOWL) && defined(HAVE_CEXPL) && defined(HAVE_CLOGL)
678 cpowl (long double complex base
, long double complex power
)
680 return cexpl (power
* clogl (base
));
685 /* sqrt(z). Algorithm pulled from glibc. */
686 #if !defined(HAVE_CSQRTF)
687 #define HAVE_CSQRTF 1
689 csqrtf (float complex z
)
700 COMPLEX_ASSIGN (v
, 0, copysignf (sqrtf (-re
), im
));
704 COMPLEX_ASSIGN (v
, fabsf (sqrtf (re
)), copysignf (0, im
));
711 r
= sqrtf (0.5 * fabsf (im
));
713 COMPLEX_ASSIGN (v
, r
, copysignf (r
, im
));
720 /* Use the identity 2 Re res Im res = Im x
721 to avoid cancellation error in d +/- Re x. */
724 r
= sqrtf (0.5 * d
+ 0.5 * re
);
729 s
= sqrtf (0.5 * d
- 0.5 * re
);
730 r
= fabsf ((0.5 * im
) / s
);
733 COMPLEX_ASSIGN (v
, r
, copysignf (s
, im
));
739 #if !defined(HAVE_CSQRT)
742 csqrt (double complex z
)
753 COMPLEX_ASSIGN (v
, 0, copysign (sqrt (-re
), im
));
757 COMPLEX_ASSIGN (v
, fabs (sqrt (re
)), copysign (0, im
));
764 r
= sqrt (0.5 * fabs (im
));
766 COMPLEX_ASSIGN (v
, r
, copysign (r
, im
));
773 /* Use the identity 2 Re res Im res = Im x
774 to avoid cancellation error in d +/- Re x. */
777 r
= sqrt (0.5 * d
+ 0.5 * re
);
782 s
= sqrt (0.5 * d
- 0.5 * re
);
783 r
= fabs ((0.5 * im
) / s
);
786 COMPLEX_ASSIGN (v
, r
, copysign (s
, im
));
792 #if !defined(HAVE_CSQRTL) && defined(HAVE_COPYSIGNL) && defined(HAVE_SQRTL) && defined(HAVE_FABSL) && defined(HAVE_HYPOTL)
793 #define HAVE_CSQRTL 1
795 csqrtl (long double complex z
)
798 long double complex v
;
806 COMPLEX_ASSIGN (v
, 0, copysignl (sqrtl (-re
), im
));
810 COMPLEX_ASSIGN (v
, fabsl (sqrtl (re
)), copysignl (0, im
));
817 r
= sqrtl (0.5 * fabsl (im
));
819 COMPLEX_ASSIGN (v
, copysignl (r
, im
), r
);
826 /* Use the identity 2 Re res Im res = Im x
827 to avoid cancellation error in d +/- Re x. */
830 r
= sqrtl (0.5 * d
+ 0.5 * re
);
835 s
= sqrtl (0.5 * d
- 0.5 * re
);
836 r
= fabsl ((0.5 * im
) / s
);
839 COMPLEX_ASSIGN (v
, r
, copysignl (s
, im
));
846 /* sinh(a + i b) = sinh(a) cos(b) + i cosh(a) sin(b) */
847 #if !defined(HAVE_CSINHF)
848 #define HAVE_CSINHF 1
850 csinhf (float complex a
)
857 COMPLEX_ASSIGN (v
, sinhf (r
) * cosf (i
), coshf (r
) * sinf (i
));
862 #if !defined(HAVE_CSINH)
865 csinh (double complex a
)
872 COMPLEX_ASSIGN (v
, sinh (r
) * cos (i
), cosh (r
) * sin (i
));
877 #if !defined(HAVE_CSINHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
878 #define HAVE_CSINHL 1
880 csinhl (long double complex a
)
883 long double complex v
;
887 COMPLEX_ASSIGN (v
, sinhl (r
) * cosl (i
), coshl (r
) * sinl (i
));
893 /* cosh(a + i b) = cosh(a) cos(b) - i sinh(a) sin(b) */
894 #if !defined(HAVE_CCOSHF)
895 #define HAVE_CCOSHF 1
897 ccoshf (float complex a
)
904 COMPLEX_ASSIGN (v
, coshf (r
) * cosf (i
), - (sinhf (r
) * sinf (i
)));
909 #if !defined(HAVE_CCOSH)
912 ccosh (double complex a
)
919 COMPLEX_ASSIGN (v
, cosh (r
) * cos (i
), - (sinh (r
) * sin (i
)));
924 #if !defined(HAVE_CCOSHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
925 #define HAVE_CCOSHL 1
927 ccoshl (long double complex a
)
930 long double complex v
;
934 COMPLEX_ASSIGN (v
, coshl (r
) * cosl (i
), - (sinhl (r
) * sinl (i
)));
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
944 ctanhf (float complex a
)
949 rt
= tanhf (REALPART (a
));
950 it
= tanf (IMAGPART (a
));
951 COMPLEX_ASSIGN (n
, rt
, it
);
952 COMPLEX_ASSIGN (d
, 1, - (rt
* it
));
958 #if !defined(HAVE_CTANH)
961 ctanh (double complex a
)
966 rt
= tanh (REALPART (a
));
967 it
= tan (IMAGPART (a
));
968 COMPLEX_ASSIGN (n
, rt
, it
);
969 COMPLEX_ASSIGN (d
, 1, - (rt
* it
));
975 #if !defined(HAVE_CTANHL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
976 #define HAVE_CTANHL 1
978 ctanhl (long double complex a
)
981 long double complex n
, d
;
983 rt
= tanhl (REALPART (a
));
984 it
= tanl (IMAGPART (a
));
985 COMPLEX_ASSIGN (n
, rt
, it
);
986 COMPLEX_ASSIGN (d
, 1, - (rt
* it
));
993 /* sin(a + i b) = sin(a) cosh(b) + i cos(a) sinh(b) */
994 #if !defined(HAVE_CSINF)
997 csinf (float complex a
)
1004 COMPLEX_ASSIGN (v
, sinf (r
) * coshf (i
), cosf (r
) * sinhf (i
));
1009 #if !defined(HAVE_CSIN)
1012 csin (double complex a
)
1019 COMPLEX_ASSIGN (v
, sin (r
) * cosh (i
), cos (r
) * sinh (i
));
1024 #if !defined(HAVE_CSINL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1025 #define HAVE_CSINL 1
1027 csinl (long double complex a
)
1030 long double complex v
;
1034 COMPLEX_ASSIGN (v
, sinl (r
) * coshl (i
), cosl (r
) * sinhl (i
));
1040 /* cos(a + i b) = cos(a) cosh(b) - i sin(a) sinh(b) */
1041 #if !defined(HAVE_CCOSF)
1042 #define HAVE_CCOSF 1
1044 ccosf (float complex a
)
1051 COMPLEX_ASSIGN (v
, cosf (r
) * coshf (i
), - (sinf (r
) * sinhf (i
)));
1056 #if !defined(HAVE_CCOS)
1059 ccos (double complex a
)
1066 COMPLEX_ASSIGN (v
, cos (r
) * cosh (i
), - (sin (r
) * sinh (i
)));
1071 #if !defined(HAVE_CCOSL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1072 #define HAVE_CCOSL 1
1074 ccosl (long double complex a
)
1077 long double complex v
;
1081 COMPLEX_ASSIGN (v
, cosl (r
) * coshl (i
), - (sinl (r
) * sinhl (i
)));
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
1091 ctanf (float complex a
)
1096 rt
= tanf (REALPART (a
));
1097 it
= tanhf (IMAGPART (a
));
1098 COMPLEX_ASSIGN (n
, rt
, it
);
1099 COMPLEX_ASSIGN (d
, 1, - (rt
* it
));
1105 #if !defined(HAVE_CTAN)
1108 ctan (double complex a
)
1111 double complex n
, d
;
1113 rt
= tan (REALPART (a
));
1114 it
= tanh (IMAGPART (a
));
1115 COMPLEX_ASSIGN (n
, rt
, it
);
1116 COMPLEX_ASSIGN (d
, 1, - (rt
* it
));
1122 #if !defined(HAVE_CTANL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1123 #define HAVE_CTANL 1
1125 ctanl (long double complex a
)
1128 long double complex n
, d
;
1130 rt
= tanl (REALPART (a
));
1131 it
= tanhl (IMAGPART (a
));
1132 COMPLEX_ASSIGN (n
, rt
, it
);
1133 COMPLEX_ASSIGN (d
, 1, - (rt
* it
));