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"
44 return (float) acos(x);
53 return (float) asin(x);
60 atan2f(float y, float x)
62 return (float) atan2(y, x);
71 return (float) atan(x);
80 return (float) ceil(x);
84 #ifndef HAVE_COPYSIGNF
85 #define HAVE_COPYSIGNF 1
87 copysignf(float x, float y)
89 return (float) copysign(x, y);
98 return (float) cos(x);
107 return (float) cosh(x);
116 return (float) exp(x);
125 return (float) fabs(x);
130 #define HAVE_FLOORF 1
134 return (float) floor(x);
139 #define HAVE_FREXPF 1
141 frexpf(float x, int *exp)
143 return (float) frexp(x, exp);
148 #define HAVE_HYPOTF 1
150 hypotf(float x, float y)
152 return (float) hypot(x, y);
161 return (float) log(x);
166 #define HAVE_LOG10F 1
170 return (float) log10(x);
175 #define HAVE_SCALBN 1
177 scalbn(double x, int y)
179 return x * pow(FLT_RADIX, y);
184 #define HAVE_SCALBNF 1
186 scalbnf(float x, int y)
188 return (float) scalbn(x, y);
197 return (float) sin(x);
206 return (float) sinh(x);
215 return (float) sqrt(x);
224 return (float) tan(x);
233 return (float) tanh(x);
253 #define HAVE_TRUNCF 1
257 return (float) trunc (x);
261 #ifndef HAVE_NEXTAFTERF
262 #define HAVE_NEXTAFTERF 1
263 /* This is a portable implementation of nextafterf that is intended to be
264 independent of the floating point format or its in memory representation.
265 This implementation works correctly with denormalized values. */
267 nextafterf(float x, float y)
269 /* This variable is marked volatile to avoid excess precision problems
270 on some platforms, including IA-32. */
271 volatile float delta;
272 float absx, denorm_min;
274 if (isnan(x) || isnan(y))
279 return x > 0 ? __FLT_MAX__ : - __FLT_MAX__;
281 /* absx = fabsf (x); */
282 absx = (x < 0.0) ? -x : x;
284 /* __FLT_DENORM_MIN__ is non-zero iff the target supports denormals. */
285 if (__FLT_DENORM_MIN__ == 0.0f)
286 denorm_min = __FLT_MIN__;
288 denorm_min = __FLT_DENORM_MIN__;
290 if (absx < __FLT_MIN__)
297 /* Discard the fraction from x. */
298 frac = frexpf (absx, &exp);
299 delta = scalbnf (0.5f, exp);
301 /* Scale x by the epsilon of the representation. By rights we should
302 have been able to combine this with scalbnf, but some targets don't
303 get that correct with denormals. */
304 delta *= __FLT_EPSILON__;
306 /* If we're going to be reducing the absolute value of X, and doing so
307 would reduce the exponent of X, then the delta to be applied is
308 one exponent smaller. */
309 if (frac == 0.5f && (y < x) == (x > 0))
312 /* If that underflows to zero, then we're back to the minimum. */
328 powf(float x, float y)
330 return (float) pow(x, y);
334 /* Note that if fpclassify is not defined, then NaN is not handled */
336 /* Algorithm by Steven G. Kargl. */
340 /* Round to nearest integral value. If the argument is halfway between two
341 integral values then round away from zero. */
368 #define HAVE_ROUNDF 1
369 /* Round to nearest integral value. If the argument is halfway between two
370 integral values then round away from zero. */
397 #define HAVE_LOG10L 1
398 /* log10 function for long double variables. The version provided here
399 reduces the argument until it fits into a double, then use log10. */
401 log10l(long double x)
403 #if LDBL_MAX_EXP > DBL_MAX_EXP
408 if (x > 0x1p16383L) { p2_result += 16383; x /= 0x1p16383L; }
409 if (x > 0x1p8191L) { p2_result += 8191; x /= 0x1p8191L; }
410 if (x > 0x1p4095L) { p2_result += 4095; x /= 0x1p4095L; }
411 if (x > 0x1p2047L) { p2_result += 2047; x /= 0x1p2047L; }
412 if (x > 0x1p1023L) { p2_result += 1023; x /= 0x1p1023L; }
413 val = log10 ((double) x);
414 return (val + p2_result * .30102999566398119521373889472449302L);
417 #if LDBL_MIN_EXP < DBL_MIN_EXP
422 if (x < 0x1p-16380L) { p2_result += 16380; x /= 0x1p-16380L; }
423 if (x < 0x1p-8189L) { p2_result += 8189; x /= 0x1p-8189L; }
424 if (x < 0x1p-4093L) { p2_result += 4093; x /= 0x1p-4093L; }
425 if (x < 0x1p-2045L) { p2_result += 2045; x /= 0x1p-2045L; }
426 if (x < 0x1p-1021L) { p2_result += 1021; x /= 0x1p-1021L; }
427 val = fabs(log10 ((double) x));
428 return (- val - p2_result * .30102999566398119521373889472449302L);
436 #if !defined(HAVE_CABSF)
439 cabsf (float complex z)
441 return hypotf (REALPART (z), IMAGPART (z));
445 #if !defined(HAVE_CABS)
448 cabs (double complex z)
450 return hypot (REALPART (z), IMAGPART (z));
454 #if !defined(HAVE_CABSL) && defined(HAVE_HYPOTL)
457 cabsl (long double complex z)
459 return hypotl (REALPART (z), IMAGPART (z));
464 #if !defined(HAVE_CARGF)
467 cargf (float complex z)
469 return atan2f (IMAGPART (z), REALPART (z));
473 #if !defined(HAVE_CARG)
476 carg (double complex z)
478 return atan2 (IMAGPART (z), REALPART (z));
482 #if !defined(HAVE_CARGL) && defined(HAVE_ATAN2L)
485 cargl (long double complex z)
487 return atan2l (IMAGPART (z), REALPART (z));
492 /* exp(z) = exp(a)*(cos(b) + i sin(b)) */
493 #if !defined(HAVE_CEXPF)
496 cexpf (float complex z)
503 COMPLEX_ASSIGN (v, cosf (b), sinf (b));
508 #if !defined(HAVE_CEXP)
511 cexp (double complex z)
518 COMPLEX_ASSIGN (v, cos (b), sin (b));
523 #if !defined(HAVE_CEXPL) && defined(HAVE_COSL) && defined(HAVE_SINL) && defined(EXPL)
526 cexpl (long double complex z)
529 long double complex v;
533 COMPLEX_ASSIGN (v, cosl (b), sinl (b));
539 /* log(z) = log (cabs(z)) + i*carg(z) */
540 #if !defined(HAVE_CLOGF)
543 clogf (float complex z)
547 COMPLEX_ASSIGN (v, logf (cabsf (z)), cargf (z));
552 #if !defined(HAVE_CLOG)
555 clog (double complex z)
559 COMPLEX_ASSIGN (v, log (cabs (z)), carg (z));
564 #if !defined(HAVE_CLOGL) && defined(HAVE_LOGL) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
567 clogl (long double complex z)
569 long double complex v;
571 COMPLEX_ASSIGN (v, logl (cabsl (z)), cargl (z));
577 /* log10(z) = log10 (cabs(z)) + i*carg(z) */
578 #if !defined(HAVE_CLOG10F)
579 #define HAVE_CLOG10F 1
581 clog10f (float complex z)
585 COMPLEX_ASSIGN (v, log10f (cabsf (z)), cargf (z));
590 #if !defined(HAVE_CLOG10)
591 #define HAVE_CLOG10 1
593 clog10 (double complex z)
597 COMPLEX_ASSIGN (v, log10 (cabs (z)), carg (z));
602 #if !defined(HAVE_CLOG10L) && defined(HAVE_LOG10L) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
603 #define HAVE_CLOG10L 1
605 clog10l (long double complex z)
607 long double complex v;
609 COMPLEX_ASSIGN (v, log10l (cabsl (z)), cargl (z));
615 /* pow(base, power) = cexp (power * clog (base)) */
616 #if !defined(HAVE_CPOWF)
619 cpowf (float complex base, float complex power)
621 return cexpf (power * clogf (base));
625 #if !defined(HAVE_CPOW)
628 cpow (double complex base, double complex power)
630 return cexp (power * clog (base));
634 #if !defined(HAVE_CPOWL) && defined(HAVE_CEXPL) && defined(HAVE_CLOGL)
637 cpowl (long double complex base, long double complex power)
639 return cexpl (power * clogl (base));
644 /* sqrt(z). Algorithm pulled from glibc. */
645 #if !defined(HAVE_CSQRTF)
646 #define HAVE_CSQRTF 1
648 csqrtf (float complex z)
659 COMPLEX_ASSIGN (v, 0, copysignf (sqrtf (-re), im));
663 COMPLEX_ASSIGN (v, fabsf (sqrtf (re)), copysignf (0, im));
670 r = sqrtf (0.5 * fabsf (im));
672 COMPLEX_ASSIGN (v, r, copysignf (r, im));
679 /* Use the identity 2 Re res Im res = Im x
680 to avoid cancellation error in d +/- Re x. */
683 r = sqrtf (0.5 * d + 0.5 * re);
688 s = sqrtf (0.5 * d - 0.5 * re);
689 r = fabsf ((0.5 * im) / s);
692 COMPLEX_ASSIGN (v, r, copysignf (s, im));
698 #if !defined(HAVE_CSQRT)
701 csqrt (double complex z)
712 COMPLEX_ASSIGN (v, 0, copysign (sqrt (-re), im));
716 COMPLEX_ASSIGN (v, fabs (sqrt (re)), copysign (0, im));
723 r = sqrt (0.5 * fabs (im));
725 COMPLEX_ASSIGN (v, r, copysign (r, im));
732 /* Use the identity 2 Re res Im res = Im x
733 to avoid cancellation error in d +/- Re x. */
736 r = sqrt (0.5 * d + 0.5 * re);
741 s = sqrt (0.5 * d - 0.5 * re);
742 r = fabs ((0.5 * im) / s);
745 COMPLEX_ASSIGN (v, r, copysign (s, im));
751 #if !defined(HAVE_CSQRTL) && defined(HAVE_COPYSIGNL) && defined(HAVE_SQRTL) && defined(HAVE_FABSL) && defined(HAVE_HYPOTL)
752 #define HAVE_CSQRTL 1
754 csqrtl (long double complex z)
757 long double complex v;
765 COMPLEX_ASSIGN (v, 0, copysignl (sqrtl (-re), im));
769 COMPLEX_ASSIGN (v, fabsl (sqrtl (re)), copysignl (0, im));
776 r = sqrtl (0.5 * fabsl (im));
778 COMPLEX_ASSIGN (v, copysignl (r, im), r);
785 /* Use the identity 2 Re res Im res = Im x
786 to avoid cancellation error in d +/- Re x. */
789 r = sqrtl (0.5 * d + 0.5 * re);
794 s = sqrtl (0.5 * d - 0.5 * re);
795 r = fabsl ((0.5 * im) / s);
798 COMPLEX_ASSIGN (v, r, copysignl (s, im));
805 /* sinh(a + i b) = sinh(a) cos(b) + i cosh(a) sin(b) */
806 #if !defined(HAVE_CSINHF)
807 #define HAVE_CSINHF 1
809 csinhf (float complex a)
816 COMPLEX_ASSIGN (v, sinhf (r) * cosf (i), coshf (r) * sinf (i));
821 #if !defined(HAVE_CSINH)
824 csinh (double complex a)
831 COMPLEX_ASSIGN (v, sinh (r) * cos (i), cosh (r) * sin (i));
836 #if !defined(HAVE_CSINHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
837 #define HAVE_CSINHL 1
839 csinhl (long double complex a)
842 long double complex v;
846 COMPLEX_ASSIGN (v, sinhl (r) * cosl (i), coshl (r) * sinl (i));
852 /* cosh(a + i b) = cosh(a) cos(b) - i sinh(a) sin(b) */
853 #if !defined(HAVE_CCOSHF)
854 #define HAVE_CCOSHF 1
856 ccoshf (float complex a)
863 COMPLEX_ASSIGN (v, coshf (r) * cosf (i), - (sinhf (r) * sinf (i)));
868 #if !defined(HAVE_CCOSH)
871 ccosh (double complex a)
878 COMPLEX_ASSIGN (v, cosh (r) * cos (i), - (sinh (r) * sin (i)));
883 #if !defined(HAVE_CCOSHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
884 #define HAVE_CCOSHL 1
886 ccoshl (long double complex a)
889 long double complex v;
893 COMPLEX_ASSIGN (v, coshl (r) * cosl (i), - (sinhl (r) * sinl (i)));
899 /* tanh(a + i b) = (tanh(a) + i tan(b)) / (1 - i tanh(a) tan(b)) */
900 #if !defined(HAVE_CTANHF)
901 #define HAVE_CTANHF 1
903 ctanhf (float complex a)
908 rt = tanhf (REALPART (a));
909 it = tanf (IMAGPART (a));
910 COMPLEX_ASSIGN (n, rt, it);
911 COMPLEX_ASSIGN (d, 1, - (rt * it));
917 #if !defined(HAVE_CTANH)
920 ctanh (double complex a)
925 rt = tanh (REALPART (a));
926 it = tan (IMAGPART (a));
927 COMPLEX_ASSIGN (n, rt, it);
928 COMPLEX_ASSIGN (d, 1, - (rt * it));
934 #if !defined(HAVE_CTANHL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
935 #define HAVE_CTANHL 1
937 ctanhl (long double complex a)
940 long double complex n, d;
942 rt = tanhl (REALPART (a));
943 it = tanl (IMAGPART (a));
944 COMPLEX_ASSIGN (n, rt, it);
945 COMPLEX_ASSIGN (d, 1, - (rt * it));
952 /* sin(a + i b) = sin(a) cosh(b) + i cos(a) sinh(b) */
953 #if !defined(HAVE_CSINF)
956 csinf (float complex a)
963 COMPLEX_ASSIGN (v, sinf (r) * coshf (i), cosf (r) * sinhf (i));
968 #if !defined(HAVE_CSIN)
971 csin (double complex a)
978 COMPLEX_ASSIGN (v, sin (r) * cosh (i), cos (r) * sinh (i));
983 #if !defined(HAVE_CSINL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
986 csinl (long double complex a)
989 long double complex v;
993 COMPLEX_ASSIGN (v, sinl (r) * coshl (i), cosl (r) * sinhl (i));
999 /* cos(a + i b) = cos(a) cosh(b) - i sin(a) sinh(b) */
1000 #if !defined(HAVE_CCOSF)
1001 #define HAVE_CCOSF 1
1003 ccosf (float complex a)
1010 COMPLEX_ASSIGN (v, cosf (r) * coshf (i), - (sinf (r) * sinhf (i)));
1015 #if !defined(HAVE_CCOS)
1018 ccos (double complex a)
1025 COMPLEX_ASSIGN (v, cos (r) * cosh (i), - (sin (r) * sinh (i)));
1030 #if !defined(HAVE_CCOSL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1031 #define HAVE_CCOSL 1
1033 ccosl (long double complex a)
1036 long double complex v;
1040 COMPLEX_ASSIGN (v, cosl (r) * coshl (i), - (sinl (r) * sinhl (i)));
1046 /* tan(a + i b) = (tan(a) + i tanh(b)) / (1 - i tan(a) tanh(b)) */
1047 #if !defined(HAVE_CTANF)
1048 #define HAVE_CTANF 1
1050 ctanf (float complex a)
1055 rt = tanf (REALPART (a));
1056 it = tanhf (IMAGPART (a));
1057 COMPLEX_ASSIGN (n, rt, it);
1058 COMPLEX_ASSIGN (d, 1, - (rt * it));
1064 #if !defined(HAVE_CTAN)
1067 ctan (double complex a)
1070 double complex n, d;
1072 rt = tan (REALPART (a));
1073 it = tanh (IMAGPART (a));
1074 COMPLEX_ASSIGN (n, rt, it);
1075 COMPLEX_ASSIGN (d, 1, - (rt * it));
1081 #if !defined(HAVE_CTANL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1082 #define HAVE_CTANL 1
1084 ctanl (long double complex a)
1087 long double complex n, d;
1089 rt = tanl (REALPART (a));
1090 it = tanhl (IMAGPART (a));
1091 COMPLEX_ASSIGN (n, rt, it);
1092 COMPLEX_ASSIGN (d, 1, - (rt * it));