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. */
32 #define C99_PROTOS_H WE_DONT_WANT_PROTOS_NOW
33 #include "libgfortran.h"
35 /* IRIX's <math.h> declares a non-C99 compliant implementation of cabs,
36 which takes two floating point arguments instead of a single complex.
37 If <complex.h> is missing this prevents building of c99_functions.c.
38 To work around this we redirect cabs{,f,l} calls to __gfc_cabs{,f,l}. */
40 #if defined(__sgi__) && !defined(HAVE_COMPLEX_H)
44 #define cabs __gfc_cabs
45 #define cabsf __gfc_cabsf
46 #define cabsl __gfc_cabsl
49 /* Tru64's <math.h> declares a non-C99 compliant implementation of cabs,
50 which takes two floating point arguments instead of a single complex.
51 To work around this we redirect cabs{,f,l} calls to __gfc_cabs{,f,l}. */
57 #define cabs __gfc_cabs
58 #define cabsf __gfc_cabsf
59 #define cabsl __gfc_cabsl
62 /* Prototypes to silence -Wstrict-prototypes -Wmissing-prototypes. */
64 float cabsf(float complex);
65 double cabs(double complex);
66 long double cabsl(long double complex);
68 float cargf(float complex);
69 double carg(double complex);
70 long double cargl(long double complex);
72 float complex clog10f(float complex);
73 double complex clog10(double complex);
74 long double complex clog10l(long double complex);
77 /* Wrappers for systems without the various C99 single precision Bessel
80 #if defined(HAVE_J0) && ! defined(HAVE_J0F)
82 extern float j0f (float);
87 return (float) j0 ((double) x);
91 #if defined(HAVE_J1) && !defined(HAVE_J1F)
93 extern float j1f (float);
97 return (float) j1 ((double) x);
101 #if defined(HAVE_JN) && !defined(HAVE_JNF)
103 extern float jnf (int, float);
108 return (float) jn (n, (double) x);
112 #if defined(HAVE_Y0) && !defined(HAVE_Y0F)
114 extern float y0f (float);
119 return (float) y0 ((double) x);
123 #if defined(HAVE_Y1) && !defined(HAVE_Y1F)
125 extern float y1f (float);
130 return (float) y1 ((double) x);
134 #if defined(HAVE_YN) && !defined(HAVE_YNF)
136 extern float ynf (int, float);
141 return (float) yn (n, (double) x);
146 /* Wrappers for systems without the C99 erff() and erfcf() functions. */
148 #if defined(HAVE_ERF) && !defined(HAVE_ERFF)
150 extern float erff (float);
155 return (float) erf ((double) x);
159 #if defined(HAVE_ERFC) && !defined(HAVE_ERFCF)
161 extern float erfcf (float);
166 return (float) erfc ((double) x);
176 return (float) acos(x);
180 #if HAVE_ACOSH && !HAVE_ACOSHF
184 return (float) acosh ((double) x);
193 return (float) asin(x);
197 #if HAVE_ASINH && !HAVE_ASINHF
201 return (float) asinh ((double) x);
206 #define HAVE_ATAN2F 1
208 atan2f(float y, float x)
210 return (float) atan2(y, x);
219 return (float) atan(x);
223 #if HAVE_ATANH && !HAVE_ATANHF
227 return (float) atanh ((double) x);
236 return (float) ceil(x);
240 #ifndef HAVE_COPYSIGNF
241 #define HAVE_COPYSIGNF 1
243 copysignf(float x, float y)
245 return (float) copysign(x, y);
254 return (float) cos(x);
263 return (float) cosh(x);
272 return (float) exp(x);
281 return (float) fabs(x);
286 #define HAVE_FLOORF 1
290 return (float) floor(x);
297 fmodf (float x, float y)
299 return (float) fmod (x, y);
304 #define HAVE_FREXPF 1
306 frexpf(float x, int *exp)
308 return (float) frexp(x, exp);
313 #define HAVE_HYPOTF 1
315 hypotf(float x, float y)
317 return (float) hypot(x, y);
326 return (float) log(x);
331 #define HAVE_LOG10F 1
335 return (float) log10(x);
340 #define HAVE_SCALBN 1
342 scalbn(double x, int y)
344 #if (FLT_RADIX == 2) && defined(HAVE_LDEXP)
347 return x * pow(FLT_RADIX, y);
353 #define HAVE_SCALBNF 1
355 scalbnf(float x, int y)
357 return (float) scalbn(x, y);
366 return (float) sin(x);
375 return (float) sinh(x);
384 return (float) sqrt(x);
393 return (float) tan(x);
402 return (float) tanh(x);
422 #define HAVE_TRUNCF 1
426 return (float) trunc (x);
430 #ifndef HAVE_NEXTAFTERF
431 #define HAVE_NEXTAFTERF 1
432 /* This is a portable implementation of nextafterf that is intended to be
433 independent of the floating point format or its in memory representation.
434 This implementation works correctly with denormalized values. */
436 nextafterf(float x, float y)
438 /* This variable is marked volatile to avoid excess precision problems
439 on some platforms, including IA-32. */
440 volatile float delta;
441 float absx, denorm_min;
443 if (isnan(x) || isnan(y))
448 return x > 0 ? __FLT_MAX__ : - __FLT_MAX__;
450 /* absx = fabsf (x); */
451 absx = (x < 0.0) ? -x : x;
453 /* __FLT_DENORM_MIN__ is non-zero iff the target supports denormals. */
454 if (__FLT_DENORM_MIN__ == 0.0f)
455 denorm_min = __FLT_MIN__;
457 denorm_min = __FLT_DENORM_MIN__;
459 if (absx < __FLT_MIN__)
466 /* Discard the fraction from x. */
467 frac = frexpf (absx, &exp);
468 delta = scalbnf (0.5f, exp);
470 /* Scale x by the epsilon of the representation. By rights we should
471 have been able to combine this with scalbnf, but some targets don't
472 get that correct with denormals. */
473 delta *= __FLT_EPSILON__;
475 /* If we're going to be reducing the absolute value of X, and doing so
476 would reduce the exponent of X, then the delta to be applied is
477 one exponent smaller. */
478 if (frac == 0.5f && (y < x) == (x > 0))
481 /* If that underflows to zero, then we're back to the minimum. */
497 powf(float x, float y)
499 return (float) pow(x, y);
503 /* Note that if fpclassify is not defined, then NaN is not handled */
505 /* Algorithm by Steven G. Kargl. */
507 #if !defined(HAVE_ROUNDL)
508 #define HAVE_ROUNDL 1
509 #if defined(HAVE_CEILL)
510 /* Round to nearest integral value. If the argument is halfway between two
511 integral values then round away from zero. */
514 roundl(long double x)
537 /* Poor version of roundl for system that don't have ceill. */
539 roundl(long double x)
541 if (x > DBL_MAX || x < -DBL_MAX)
543 #ifdef HAVE_NEXTAFTERL
544 static long double prechalf = nexafterl (0.5L, LDBL_MAX);
546 static long double prechalf = 0.5L;
548 return (GFC_INTEGER_LARGEST) (x + (x > 0 ? prechalf : -prechalf));
552 return round((double) x);
560 /* Round to nearest integral value. If the argument is halfway between two
561 integral values then round away from zero. */
588 #define HAVE_ROUNDF 1
589 /* Round to nearest integral value. If the argument is halfway between two
590 integral values then round away from zero. */
617 /* lround{f,,l} and llround{f,,l} functions. */
619 #if !defined(HAVE_LROUNDF) && defined(HAVE_ROUNDF)
620 #define HAVE_LROUNDF 1
624 return (long int) roundf (x);
628 #if !defined(HAVE_LROUND) && defined(HAVE_ROUND)
629 #define HAVE_LROUND 1
633 return (long int) round (x);
637 #if !defined(HAVE_LROUNDL) && defined(HAVE_ROUNDL)
638 #define HAVE_LROUNDL 1
640 lroundl (long double x)
642 return (long long int) roundl (x);
646 #if !defined(HAVE_LLROUNDF) && defined(HAVE_ROUNDF)
647 #define HAVE_LLROUNDF 1
651 return (long long int) roundf (x);
655 #if !defined(HAVE_LLROUND) && defined(HAVE_ROUND)
656 #define HAVE_LLROUND 1
660 return (long long int) round (x);
664 #if !defined(HAVE_LLROUNDL) && defined(HAVE_ROUNDL)
665 #define HAVE_LLROUNDL 1
667 llroundl (long double x)
669 return (long long int) roundl (x);
675 #define HAVE_LOG10L 1
676 /* log10 function for long double variables. The version provided here
677 reduces the argument until it fits into a double, then use log10. */
679 log10l(long double x)
681 #if LDBL_MAX_EXP > DBL_MAX_EXP
686 if (x > 0x1p16383L) { p2_result += 16383; x /= 0x1p16383L; }
687 if (x > 0x1p8191L) { p2_result += 8191; x /= 0x1p8191L; }
688 if (x > 0x1p4095L) { p2_result += 4095; x /= 0x1p4095L; }
689 if (x > 0x1p2047L) { p2_result += 2047; x /= 0x1p2047L; }
690 if (x > 0x1p1023L) { p2_result += 1023; x /= 0x1p1023L; }
691 val = log10 ((double) x);
692 return (val + p2_result * .30102999566398119521373889472449302L);
695 #if LDBL_MIN_EXP < DBL_MIN_EXP
700 if (x < 0x1p-16380L) { p2_result += 16380; x /= 0x1p-16380L; }
701 if (x < 0x1p-8189L) { p2_result += 8189; x /= 0x1p-8189L; }
702 if (x < 0x1p-4093L) { p2_result += 4093; x /= 0x1p-4093L; }
703 if (x < 0x1p-2045L) { p2_result += 2045; x /= 0x1p-2045L; }
704 if (x < 0x1p-1021L) { p2_result += 1021; x /= 0x1p-1021L; }
705 val = fabs(log10 ((double) x));
706 return (- val - p2_result * .30102999566398119521373889472449302L);
715 #define HAVE_FLOORL 1
717 floorl (long double x)
719 /* Zero, possibly signed. */
723 /* Large magnitude. */
724 if (x > DBL_MAX || x < (-DBL_MAX))
727 /* Small positive values. */
728 if (x >= 0 && x < DBL_MIN)
731 /* Small negative values. */
732 if (x < 0 && x > (-DBL_MIN))
743 fmodl (long double x, long double y)
748 /* Need to check that the result has the same sign as x and magnitude
749 less than the magnitude of y. */
750 return x - floorl (x / y) * y;
755 #if !defined(HAVE_CABSF)
758 cabsf (float complex z)
760 return hypotf (REALPART (z), IMAGPART (z));
764 #if !defined(HAVE_CABS)
767 cabs (double complex z)
769 return hypot (REALPART (z), IMAGPART (z));
773 #if !defined(HAVE_CABSL) && defined(HAVE_HYPOTL)
776 cabsl (long double complex z)
778 return hypotl (REALPART (z), IMAGPART (z));
783 #if !defined(HAVE_CARGF)
786 cargf (float complex z)
788 return atan2f (IMAGPART (z), REALPART (z));
792 #if !defined(HAVE_CARG)
795 carg (double complex z)
797 return atan2 (IMAGPART (z), REALPART (z));
801 #if !defined(HAVE_CARGL) && defined(HAVE_ATAN2L)
804 cargl (long double complex z)
806 return atan2l (IMAGPART (z), REALPART (z));
811 /* exp(z) = exp(a)*(cos(b) + i sin(b)) */
812 #if !defined(HAVE_CEXPF)
815 cexpf (float complex z)
822 COMPLEX_ASSIGN (v, cosf (b), sinf (b));
827 #if !defined(HAVE_CEXP)
830 cexp (double complex z)
837 COMPLEX_ASSIGN (v, cos (b), sin (b));
842 #if !defined(HAVE_CEXPL) && defined(HAVE_COSL) && defined(HAVE_SINL) && defined(EXPL)
845 cexpl (long double complex z)
848 long double complex v;
852 COMPLEX_ASSIGN (v, cosl (b), sinl (b));
858 /* log(z) = log (cabs(z)) + i*carg(z) */
859 #if !defined(HAVE_CLOGF)
862 clogf (float complex z)
866 COMPLEX_ASSIGN (v, logf (cabsf (z)), cargf (z));
871 #if !defined(HAVE_CLOG)
874 clog (double complex z)
878 COMPLEX_ASSIGN (v, log (cabs (z)), carg (z));
883 #if !defined(HAVE_CLOGL) && defined(HAVE_LOGL) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
886 clogl (long double complex z)
888 long double complex v;
890 COMPLEX_ASSIGN (v, logl (cabsl (z)), cargl (z));
896 /* log10(z) = log10 (cabs(z)) + i*carg(z) */
897 #if !defined(HAVE_CLOG10F)
898 #define HAVE_CLOG10F 1
900 clog10f (float complex z)
904 COMPLEX_ASSIGN (v, log10f (cabsf (z)), cargf (z));
909 #if !defined(HAVE_CLOG10)
910 #define HAVE_CLOG10 1
912 clog10 (double complex z)
916 COMPLEX_ASSIGN (v, log10 (cabs (z)), carg (z));
921 #if !defined(HAVE_CLOG10L) && defined(HAVE_LOG10L) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
922 #define HAVE_CLOG10L 1
924 clog10l (long double complex z)
926 long double complex v;
928 COMPLEX_ASSIGN (v, log10l (cabsl (z)), cargl (z));
934 /* pow(base, power) = cexp (power * clog (base)) */
935 #if !defined(HAVE_CPOWF)
938 cpowf (float complex base, float complex power)
940 return cexpf (power * clogf (base));
944 #if !defined(HAVE_CPOW)
947 cpow (double complex base, double complex power)
949 return cexp (power * clog (base));
953 #if !defined(HAVE_CPOWL) && defined(HAVE_CEXPL) && defined(HAVE_CLOGL)
956 cpowl (long double complex base, long double complex power)
958 return cexpl (power * clogl (base));
963 /* sqrt(z). Algorithm pulled from glibc. */
964 #if !defined(HAVE_CSQRTF)
965 #define HAVE_CSQRTF 1
967 csqrtf (float complex z)
978 COMPLEX_ASSIGN (v, 0, copysignf (sqrtf (-re), im));
982 COMPLEX_ASSIGN (v, fabsf (sqrtf (re)), copysignf (0, im));
989 r = sqrtf (0.5 * fabsf (im));
991 COMPLEX_ASSIGN (v, r, copysignf (r, im));
998 /* Use the identity 2 Re res Im res = Im x
999 to avoid cancellation error in d +/- Re x. */
1002 r = sqrtf (0.5 * d + 0.5 * re);
1007 s = sqrtf (0.5 * d - 0.5 * re);
1008 r = fabsf ((0.5 * im) / s);
1011 COMPLEX_ASSIGN (v, r, copysignf (s, im));
1017 #if !defined(HAVE_CSQRT)
1018 #define HAVE_CSQRT 1
1020 csqrt (double complex z)
1031 COMPLEX_ASSIGN (v, 0, copysign (sqrt (-re), im));
1035 COMPLEX_ASSIGN (v, fabs (sqrt (re)), copysign (0, im));
1042 r = sqrt (0.5 * fabs (im));
1044 COMPLEX_ASSIGN (v, r, copysign (r, im));
1051 /* Use the identity 2 Re res Im res = Im x
1052 to avoid cancellation error in d +/- Re x. */
1055 r = sqrt (0.5 * d + 0.5 * re);
1060 s = sqrt (0.5 * d - 0.5 * re);
1061 r = fabs ((0.5 * im) / s);
1064 COMPLEX_ASSIGN (v, r, copysign (s, im));
1070 #if !defined(HAVE_CSQRTL) && defined(HAVE_COPYSIGNL) && defined(HAVE_SQRTL) && defined(HAVE_FABSL) && defined(HAVE_HYPOTL)
1071 #define HAVE_CSQRTL 1
1073 csqrtl (long double complex z)
1076 long double complex v;
1084 COMPLEX_ASSIGN (v, 0, copysignl (sqrtl (-re), im));
1088 COMPLEX_ASSIGN (v, fabsl (sqrtl (re)), copysignl (0, im));
1095 r = sqrtl (0.5 * fabsl (im));
1097 COMPLEX_ASSIGN (v, copysignl (r, im), r);
1101 long double d, r, s;
1103 d = hypotl (re, im);
1104 /* Use the identity 2 Re res Im res = Im x
1105 to avoid cancellation error in d +/- Re x. */
1108 r = sqrtl (0.5 * d + 0.5 * re);
1113 s = sqrtl (0.5 * d - 0.5 * re);
1114 r = fabsl ((0.5 * im) / s);
1117 COMPLEX_ASSIGN (v, r, copysignl (s, im));
1124 /* sinh(a + i b) = sinh(a) cos(b) + i cosh(a) sin(b) */
1125 #if !defined(HAVE_CSINHF)
1126 #define HAVE_CSINHF 1
1128 csinhf (float complex a)
1135 COMPLEX_ASSIGN (v, sinhf (r) * cosf (i), coshf (r) * sinf (i));
1140 #if !defined(HAVE_CSINH)
1141 #define HAVE_CSINH 1
1143 csinh (double complex a)
1150 COMPLEX_ASSIGN (v, sinh (r) * cos (i), cosh (r) * sin (i));
1155 #if !defined(HAVE_CSINHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1156 #define HAVE_CSINHL 1
1158 csinhl (long double complex a)
1161 long double complex v;
1165 COMPLEX_ASSIGN (v, sinhl (r) * cosl (i), coshl (r) * sinl (i));
1171 /* cosh(a + i b) = cosh(a) cos(b) - i sinh(a) sin(b) */
1172 #if !defined(HAVE_CCOSHF)
1173 #define HAVE_CCOSHF 1
1175 ccoshf (float complex a)
1182 COMPLEX_ASSIGN (v, coshf (r) * cosf (i), - (sinhf (r) * sinf (i)));
1187 #if !defined(HAVE_CCOSH)
1188 #define HAVE_CCOSH 1
1190 ccosh (double complex a)
1197 COMPLEX_ASSIGN (v, cosh (r) * cos (i), - (sinh (r) * sin (i)));
1202 #if !defined(HAVE_CCOSHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1203 #define HAVE_CCOSHL 1
1205 ccoshl (long double complex a)
1208 long double complex v;
1212 COMPLEX_ASSIGN (v, coshl (r) * cosl (i), - (sinhl (r) * sinl (i)));
1218 /* tanh(a + i b) = (tanh(a) + i tan(b)) / (1 - i tanh(a) tan(b)) */
1219 #if !defined(HAVE_CTANHF)
1220 #define HAVE_CTANHF 1
1222 ctanhf (float complex a)
1227 rt = tanhf (REALPART (a));
1228 it = tanf (IMAGPART (a));
1229 COMPLEX_ASSIGN (n, rt, it);
1230 COMPLEX_ASSIGN (d, 1, - (rt * it));
1236 #if !defined(HAVE_CTANH)
1237 #define HAVE_CTANH 1
1239 ctanh (double complex a)
1242 double complex n, d;
1244 rt = tanh (REALPART (a));
1245 it = tan (IMAGPART (a));
1246 COMPLEX_ASSIGN (n, rt, it);
1247 COMPLEX_ASSIGN (d, 1, - (rt * it));
1253 #if !defined(HAVE_CTANHL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1254 #define HAVE_CTANHL 1
1256 ctanhl (long double complex a)
1259 long double complex n, d;
1261 rt = tanhl (REALPART (a));
1262 it = tanl (IMAGPART (a));
1263 COMPLEX_ASSIGN (n, rt, it);
1264 COMPLEX_ASSIGN (d, 1, - (rt * it));
1271 /* sin(a + i b) = sin(a) cosh(b) + i cos(a) sinh(b) */
1272 #if !defined(HAVE_CSINF)
1273 #define HAVE_CSINF 1
1275 csinf (float complex a)
1282 COMPLEX_ASSIGN (v, sinf (r) * coshf (i), cosf (r) * sinhf (i));
1287 #if !defined(HAVE_CSIN)
1290 csin (double complex a)
1297 COMPLEX_ASSIGN (v, sin (r) * cosh (i), cos (r) * sinh (i));
1302 #if !defined(HAVE_CSINL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1303 #define HAVE_CSINL 1
1305 csinl (long double complex a)
1308 long double complex v;
1312 COMPLEX_ASSIGN (v, sinl (r) * coshl (i), cosl (r) * sinhl (i));
1318 /* cos(a + i b) = cos(a) cosh(b) - i sin(a) sinh(b) */
1319 #if !defined(HAVE_CCOSF)
1320 #define HAVE_CCOSF 1
1322 ccosf (float complex a)
1329 COMPLEX_ASSIGN (v, cosf (r) * coshf (i), - (sinf (r) * sinhf (i)));
1334 #if !defined(HAVE_CCOS)
1337 ccos (double complex a)
1344 COMPLEX_ASSIGN (v, cos (r) * cosh (i), - (sin (r) * sinh (i)));
1349 #if !defined(HAVE_CCOSL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1350 #define HAVE_CCOSL 1
1352 ccosl (long double complex a)
1355 long double complex v;
1359 COMPLEX_ASSIGN (v, cosl (r) * coshl (i), - (sinl (r) * sinhl (i)));
1365 /* tan(a + i b) = (tan(a) + i tanh(b)) / (1 - i tan(a) tanh(b)) */
1366 #if !defined(HAVE_CTANF)
1367 #define HAVE_CTANF 1
1369 ctanf (float complex a)
1374 rt = tanf (REALPART (a));
1375 it = tanhf (IMAGPART (a));
1376 COMPLEX_ASSIGN (n, rt, it);
1377 COMPLEX_ASSIGN (d, 1, - (rt * it));
1383 #if !defined(HAVE_CTAN)
1386 ctan (double complex a)
1389 double complex n, d;
1391 rt = tan (REALPART (a));
1392 it = tanh (IMAGPART (a));
1393 COMPLEX_ASSIGN (n, rt, it);
1394 COMPLEX_ASSIGN (d, 1, - (rt * it));
1400 #if !defined(HAVE_CTANL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1401 #define HAVE_CTANL 1
1403 ctanl (long double complex a)
1406 long double complex n, d;
1408 rt = tanl (REALPART (a));
1409 it = tanhl (IMAGPART (a));
1410 COMPLEX_ASSIGN (n, rt, it);
1411 COMPLEX_ASSIGN (d, 1, - (rt * it));