1 /* Implementation of various C99 functions
2 Copyright (C) 2004-2024 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 3 of the License, or (at your option) any later version.
11 Libgfortran is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
16 Under Section 7 of GPL version 3, you are granted additional
17 permissions described in the GCC Runtime Library Exception, version
18 3.1, as published by the Free Software Foundation.
20 You should have received a copy of the GNU General Public License and
21 a copy of the GCC Runtime Library Exception along with this program;
22 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
23 <http://www.gnu.org/licenses/>. */
27 #define C99_PROTOS_H WE_DONT_WANT_PROTOS_NOW
28 #include "libgfortran.h"
30 /* On a C99 system "I" (with I*I = -1) should be defined in complex.h;
31 if not, we define a fallback version here. */
33 # if defined(_Imaginary_I)
34 # define I _Imaginary_I
35 # elif defined(_Complex_I)
42 /* Macros to get real and imaginary parts of a complex, and set
44 #define REALPART(z) (__real__(z))
45 #define IMAGPART(z) (__imag__(z))
46 #define COMPLEX_ASSIGN(z_, r_, i_) {__real__(z_) = (r_); __imag__(z_) = (i_);}
49 /* Prototypes are included to silence -Wstrict-prototypes
50 -Wmissing-prototypes. */
53 /* Wrappers for systems without the various C99 single precision Bessel
56 #if defined(HAVE_J0) && ! defined(HAVE_J0F)
63 return (float) j0 ((double) x
);
67 #if defined(HAVE_J1) && !defined(HAVE_J1F)
73 return (float) j1 ((double) x
);
77 #if defined(HAVE_JN) && !defined(HAVE_JNF)
79 float jnf (int, float);
84 return (float) jn (n
, (double) x
);
88 #if defined(HAVE_Y0) && !defined(HAVE_Y0F)
95 return (float) y0 ((double) x
);
99 #if defined(HAVE_Y1) && !defined(HAVE_Y1F)
106 return (float) y1 ((double) x
);
110 #if defined(HAVE_YN) && !defined(HAVE_YNF)
112 float ynf (int, float);
117 return (float) yn (n
, (double) x
);
122 /* Wrappers for systems without the C99 erff() and erfcf() functions. */
124 #if defined(HAVE_ERF) && !defined(HAVE_ERFF)
131 return (float) erf ((double) x
);
135 #if defined(HAVE_ERFC) && !defined(HAVE_ERFCF)
142 return (float) erfc ((double) x
);
149 float acosf (float x
);
154 return (float) acos (x
);
158 #if HAVE_ACOSH && !HAVE_ACOSHF
159 float acoshf (float x
);
164 return (float) acosh ((double) x
);
170 float asinf (float x
);
175 return (float) asin (x
);
179 #if HAVE_ASINH && !HAVE_ASINHF
180 float asinhf (float x
);
185 return (float) asinh ((double) x
);
190 #define HAVE_ATAN2F 1
191 float atan2f (float y
, float x
);
194 atan2f (float y
, float x
)
196 return (float) atan2 (y
, x
);
202 float atanf (float x
);
207 return (float) atan (x
);
211 #if HAVE_ATANH && !HAVE_ATANHF
212 float atanhf (float x
);
217 return (float) atanh ((double) x
);
223 float ceilf (float x
);
228 return (float) ceil (x
);
232 #if !defined(HAVE_COPYSIGN) && defined(HAVE_INLINE_BUILTIN_COPYSIGN)
233 #define HAVE_COPYSIGN 1
234 double copysign (double x
, double y
);
237 copysign (double x
, double y
)
239 return __builtin_copysign (x
, y
);
243 #ifndef HAVE_COPYSIGNF
244 #define HAVE_COPYSIGNF 1
245 float copysignf (float x
, float y
);
248 copysignf (float x
, float y
)
250 return (float) copysign (x
, y
);
254 #if !defined(HAVE_COPYSIGNL) && defined(HAVE_INLINE_BUILTIN_COPYSIGNL)
255 #define HAVE_COPYSIGNL 1
256 long double copysignl (long double x
, long double y
);
259 copysignl (long double x
, long double y
)
261 return __builtin_copysignl (x
, y
);
267 float cosf (float x
);
272 return (float) cos (x
);
278 float coshf (float x
);
283 return (float) cosh (x
);
289 float expf (float x
);
294 return (float) exp (x
);
298 #if !defined(HAVE_FABS) && defined(HAVE_INLINE_BUILTIN_FABS)
300 double fabs (double x
);
305 return __builtin_fabs (x
);
311 float fabsf (float x
);
316 return (float) fabs (x
);
320 #if !defined(HAVE_FABSL) && defined(HAVE_INLINE_BUILTIN_FABSL)
322 long double fabsl (long double x
);
325 fabsl (long double x
)
327 return __builtin_fabsl (x
);
332 #define HAVE_FLOORF 1
333 float floorf (float x
);
338 return (float) floor (x
);
344 float fmodf (float x
, float y
);
347 fmodf (float x
, float y
)
349 return (float) fmod (x
, y
);
354 #define HAVE_FREXPF 1
355 float frexpf (float x
, int *exp
);
358 frexpf (float x
, int *exp
)
360 return (float) frexp (x
, exp
);
365 #define HAVE_HYPOTF 1
366 float hypotf (float x
, float y
);
369 hypotf (float x
, float y
)
371 return (float) hypot (x
, y
);
377 float logf (float x
);
382 return (float) log (x
);
387 #define HAVE_LOG10F 1
388 float log10f (float x
);
393 return (float) log10 (x
);
398 #define HAVE_SCALBN 1
399 double scalbn (double x
, int y
);
402 scalbn (double x
, int y
)
404 #if (FLT_RADIX == 2) && defined(HAVE_LDEXP)
407 return x
* pow (FLT_RADIX
, y
);
413 #define HAVE_SCALBNF 1
414 float scalbnf (float x
, int y
);
417 scalbnf (float x
, int y
)
419 return (float) scalbn (x
, y
);
425 float sinf (float x
);
430 return (float) sin (x
);
436 float sinhf (float x
);
441 return (float) sinh (x
);
447 float sqrtf (float x
);
452 return (float) sqrt (x
);
458 float tanf (float x
);
463 return (float) tan (x
);
469 float tanhf (float x
);
474 return (float) tanh (x
);
480 double trunc (double x
);
496 #define HAVE_TRUNCF 1
497 float truncf (float x
);
502 return (float) trunc (x
);
506 #ifndef HAVE_NEXTAFTERF
507 #define HAVE_NEXTAFTERF 1
508 /* This is a portable implementation of nextafterf that is intended to be
509 independent of the floating point format or its in memory representation.
510 This implementation works correctly with denormalized values. */
511 float nextafterf (float x
, float y
);
514 nextafterf (float x
, float y
)
516 /* This variable is marked volatile to avoid excess precision problems
517 on some platforms, including IA-32. */
518 volatile float delta
;
519 float absx
, denorm_min
;
521 if (isnan (x
) || isnan (y
))
526 return x
> 0 ? __FLT_MAX__
: - __FLT_MAX__
;
528 /* absx = fabsf (x); */
529 absx
= (x
< 0.0) ? -x
: x
;
531 /* __FLT_DENORM_MIN__ is non-zero iff the target supports denormals. */
532 if (__FLT_DENORM_MIN__
== 0.0f
)
533 denorm_min
= __FLT_MIN__
;
535 denorm_min
= __FLT_DENORM_MIN__
;
537 if (absx
< __FLT_MIN__
)
544 /* Discard the fraction from x. */
545 frac
= frexpf (absx
, &exp
);
546 delta
= scalbnf (0.5f
, exp
);
548 /* Scale x by the epsilon of the representation. By rights we should
549 have been able to combine this with scalbnf, but some targets don't
550 get that correct with denormals. */
551 delta
*= __FLT_EPSILON__
;
553 /* If we're going to be reducing the absolute value of X, and doing so
554 would reduce the exponent of X, then the delta to be applied is
555 one exponent smaller. */
556 if (frac
== 0.5f
&& (y
< x
) == (x
> 0))
559 /* If that underflows to zero, then we're back to the minimum. */
574 float powf (float x
, float y
);
577 powf (float x
, float y
)
579 return (float) pow (x
, y
);
586 /* Round to nearest integral value. If the argument is halfway between two
587 integral values then round away from zero. */
588 double round (double x
);
615 /* Algorithm by Steven G. Kargl. */
617 #if !defined(HAVE_ROUNDL)
618 #define HAVE_ROUNDL 1
619 long double roundl (long double x
);
621 #if defined(HAVE_CEILL)
622 /* Round to nearest integral value. If the argument is halfway between two
623 integral values then round away from zero. */
626 roundl (long double x
)
649 /* Poor version of roundl for system that don't have ceill. */
651 roundl (long double x
)
653 if (x
> DBL_MAX
|| x
< -DBL_MAX
)
655 #ifdef HAVE_NEXTAFTERL
656 long double prechalf
= nextafterl (0.5L, LDBL_MAX
);
658 static long double prechalf
= 0.5L;
660 return (GFC_INTEGER_LARGEST
) (x
+ (x
> 0 ? prechalf
: -prechalf
));
664 return round ((double) x
);
671 #define HAVE_ROUNDF 1
672 /* Round to nearest integral value. If the argument is halfway between two
673 integral values then round away from zero. */
674 float roundf (float x
);
701 /* lround{f,,l} and llround{f,,l} functions. */
703 #if !defined(HAVE_LROUNDF) && defined(HAVE_ROUNDF)
704 #define HAVE_LROUNDF 1
705 long int lroundf (float x
);
710 return (long int) roundf (x
);
714 #if !defined(HAVE_LROUND) && defined(HAVE_ROUND)
715 #define HAVE_LROUND 1
716 long int lround (double x
);
721 return (long int) round (x
);
725 #if !defined(HAVE_LROUNDL) && defined(HAVE_ROUNDL)
726 #define HAVE_LROUNDL 1
727 long int lroundl (long double x
);
730 lroundl (long double x
)
732 return (long long int) roundl (x
);
736 #if !defined(HAVE_LLROUNDF) && defined(HAVE_ROUNDF)
737 #define HAVE_LLROUNDF 1
738 long long int llroundf (float x
);
743 return (long long int) roundf (x
);
747 #if !defined(HAVE_LLROUND) && defined(HAVE_ROUND)
748 #define HAVE_LLROUND 1
749 long long int llround (double x
);
754 return (long long int) round (x
);
758 #if !defined(HAVE_LLROUNDL) && defined(HAVE_ROUNDL)
759 #define HAVE_LLROUNDL 1
760 long long int llroundl (long double x
);
763 llroundl (long double x
)
765 return (long long int) roundl (x
);
771 #define HAVE_LOG10L 1
772 /* log10 function for long double variables. The version provided here
773 reduces the argument until it fits into a double, then use log10. */
774 long double log10l (long double x
);
777 log10l (long double x
)
779 #if LDBL_MAX_EXP > DBL_MAX_EXP
784 if (x
> 0x1p
16383L) { p2_result
+= 16383; x
/= 0x1p
16383L; }
785 if (x
> 0x1p
8191L) { p2_result
+= 8191; x
/= 0x1p
8191L; }
786 if (x
> 0x1p
4095L) { p2_result
+= 4095; x
/= 0x1p
4095L; }
787 if (x
> 0x1p
2047L) { p2_result
+= 2047; x
/= 0x1p
2047L; }
788 if (x
> 0x1p
1023L) { p2_result
+= 1023; x
/= 0x1p
1023L; }
789 val
= log10 ((double) x
);
790 return (val
+ p2_result
* .30102999566398119521373889472449302L);
793 #if LDBL_MIN_EXP < DBL_MIN_EXP
798 if (x
< 0x1p
-16380L) { p2_result
+= 16380; x
/= 0x1p
-16380L; }
799 if (x
< 0x1p
-8189L) { p2_result
+= 8189; x
/= 0x1p
-8189L; }
800 if (x
< 0x1p
-4093L) { p2_result
+= 4093; x
/= 0x1p
-4093L; }
801 if (x
< 0x1p
-2045L) { p2_result
+= 2045; x
/= 0x1p
-2045L; }
802 if (x
< 0x1p
-1021L) { p2_result
+= 1021; x
/= 0x1p
-1021L; }
803 val
= fabs (log10 ((double) x
));
804 return (- val
- p2_result
* .30102999566398119521373889472449302L);
813 #define HAVE_FLOORL 1
814 long double floorl (long double x
);
817 floorl (long double x
)
819 /* Zero, possibly signed. */
823 /* Large magnitude. */
824 if (x
> DBL_MAX
|| x
< (-DBL_MAX
))
827 /* Small positive values. */
828 if (x
>= 0 && x
< DBL_MIN
)
831 /* Small negative values. */
832 if (x
< 0 && x
> (-DBL_MIN
))
842 long double fmodl (long double x
, long double y
);
845 fmodl (long double x
, long double y
)
850 /* Need to check that the result has the same sign as x and magnitude
851 less than the magnitude of y. */
852 return x
- floorl (x
/ y
) * y
;
857 #if !defined(HAVE_CABSF)
859 float cabsf (float complex z
);
862 cabsf (float complex z
)
864 return hypotf (REALPART (z
), IMAGPART (z
));
868 #if !defined(HAVE_CABS)
870 double cabs (double complex z
);
873 cabs (double complex z
)
875 return hypot (REALPART (z
), IMAGPART (z
));
879 #if !defined(HAVE_CABSL) && defined(HAVE_HYPOTL)
881 long double cabsl (long double complex z
);
884 cabsl (long double complex z
)
886 return hypotl (REALPART (z
), IMAGPART (z
));
891 #if !defined(HAVE_CARGF)
893 float cargf (float complex z
);
896 cargf (float complex z
)
898 return atan2f (IMAGPART (z
), REALPART (z
));
902 #if !defined(HAVE_CARG)
904 double carg (double complex z
);
907 carg (double complex z
)
909 return atan2 (IMAGPART (z
), REALPART (z
));
913 #if !defined(HAVE_CARGL) && defined(HAVE_ATAN2L)
915 long double cargl (long double complex z
);
918 cargl (long double complex z
)
920 return atan2l (IMAGPART (z
), REALPART (z
));
925 /* exp(z) = exp(a)*(cos(b) + i sin(b)) */
926 #if !defined(HAVE_CEXPF)
928 float complex cexpf (float complex z
);
931 cexpf (float complex z
)
938 COMPLEX_ASSIGN (v
, cosf (b
), sinf (b
));
943 #if !defined(HAVE_CEXP)
945 double complex cexp (double complex z
);
948 cexp (double complex z
)
955 COMPLEX_ASSIGN (v
, cos (b
), sin (b
));
960 #if !defined(HAVE_CEXPL) && defined(HAVE_COSL) && defined(HAVE_SINL) && defined(HAVE_EXPL)
962 long double complex cexpl (long double complex z
);
965 cexpl (long double complex z
)
968 long double complex v
;
972 COMPLEX_ASSIGN (v
, cosl (b
), sinl (b
));
978 /* log(z) = log (cabs(z)) + i*carg(z) */
979 #if !defined(HAVE_CLOGF)
981 float complex clogf (float complex z
);
984 clogf (float complex z
)
988 COMPLEX_ASSIGN (v
, logf (cabsf (z
)), cargf (z
));
993 #if !defined(HAVE_CLOG)
995 double complex clog (double complex z
);
998 clog (double complex z
)
1002 COMPLEX_ASSIGN (v
, log (cabs (z
)), carg (z
));
1007 #if !defined(HAVE_CLOGL) && defined(HAVE_LOGL) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
1008 #define HAVE_CLOGL 1
1009 long double complex clogl (long double complex z
);
1012 clogl (long double complex z
)
1014 long double complex v
;
1016 COMPLEX_ASSIGN (v
, logl (cabsl (z
)), cargl (z
));
1022 /* log10(z) = log10 (cabs(z)) + i*carg(z) */
1023 #if !defined(HAVE_CLOG10F)
1024 #define HAVE_CLOG10F 1
1025 float complex clog10f (float complex z
);
1028 clog10f (float complex z
)
1032 COMPLEX_ASSIGN (v
, log10f (cabsf (z
)), cargf (z
));
1037 #if !defined(HAVE_CLOG10)
1038 #define HAVE_CLOG10 1
1039 double complex clog10 (double complex z
);
1042 clog10 (double complex z
)
1046 COMPLEX_ASSIGN (v
, log10 (cabs (z
)), carg (z
));
1051 #if !defined(HAVE_CLOG10L) && defined(HAVE_LOG10L) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
1052 #define HAVE_CLOG10L 1
1053 long double complex clog10l (long double complex z
);
1056 clog10l (long double complex z
)
1058 long double complex v
;
1060 COMPLEX_ASSIGN (v
, log10l (cabsl (z
)), cargl (z
));
1066 /* pow(base, power) = cexp (power * clog (base)) */
1067 #if !defined(HAVE_CPOWF)
1068 #define HAVE_CPOWF 1
1069 float complex cpowf (float complex base
, float complex power
);
1072 cpowf (float complex base
, float complex power
)
1074 return cexpf (power
* clogf (base
));
1078 #if !defined(HAVE_CPOW)
1080 double complex cpow (double complex base
, double complex power
);
1083 cpow (double complex base
, double complex power
)
1085 return cexp (power
* clog (base
));
1089 #if !defined(HAVE_CPOWL) && defined(HAVE_CEXPL) && defined(HAVE_CLOGL)
1090 #define HAVE_CPOWL 1
1091 long double complex cpowl (long double complex base
, long double complex power
);
1094 cpowl (long double complex base
, long double complex power
)
1096 return cexpl (power
* clogl (base
));
1101 /* sqrt(z). Algorithm pulled from glibc. */
1102 #if !defined(HAVE_CSQRTF)
1103 #define HAVE_CSQRTF 1
1104 float complex csqrtf (float complex z
);
1107 csqrtf (float complex z
)
1118 COMPLEX_ASSIGN (v
, 0, copysignf (sqrtf (-re
), im
));
1122 COMPLEX_ASSIGN (v
, fabsf (sqrtf (re
)), copysignf (0, im
));
1129 r
= sqrtf (0.5 * fabsf (im
));
1131 COMPLEX_ASSIGN (v
, r
, copysignf (r
, im
));
1137 d
= hypotf (re
, im
);
1138 /* Use the identity 2 Re res Im res = Im x
1139 to avoid cancellation error in d +/- Re x. */
1142 r
= sqrtf (0.5 * d
+ 0.5 * re
);
1147 s
= sqrtf (0.5 * d
- 0.5 * re
);
1148 r
= fabsf ((0.5 * im
) / s
);
1151 COMPLEX_ASSIGN (v
, r
, copysignf (s
, im
));
1157 #if !defined(HAVE_CSQRT)
1158 #define HAVE_CSQRT 1
1159 double complex csqrt (double complex z
);
1162 csqrt (double complex z
)
1173 COMPLEX_ASSIGN (v
, 0, copysign (sqrt (-re
), im
));
1177 COMPLEX_ASSIGN (v
, fabs (sqrt (re
)), copysign (0, im
));
1184 r
= sqrt (0.5 * fabs (im
));
1186 COMPLEX_ASSIGN (v
, r
, copysign (r
, im
));
1193 /* Use the identity 2 Re res Im res = Im x
1194 to avoid cancellation error in d +/- Re x. */
1197 r
= sqrt (0.5 * d
+ 0.5 * re
);
1202 s
= sqrt (0.5 * d
- 0.5 * re
);
1203 r
= fabs ((0.5 * im
) / s
);
1206 COMPLEX_ASSIGN (v
, r
, copysign (s
, im
));
1212 #if !defined(HAVE_CSQRTL) && defined(HAVE_COPYSIGNL) && defined(HAVE_SQRTL) && defined(HAVE_FABSL) && defined(HAVE_HYPOTL)
1213 #define HAVE_CSQRTL 1
1214 long double complex csqrtl (long double complex z
);
1217 csqrtl (long double complex z
)
1220 long double complex v
;
1228 COMPLEX_ASSIGN (v
, 0, copysignl (sqrtl (-re
), im
));
1232 COMPLEX_ASSIGN (v
, fabsl (sqrtl (re
)), copysignl (0, im
));
1239 r
= sqrtl (0.5 * fabsl (im
));
1241 COMPLEX_ASSIGN (v
, copysignl (r
, im
), r
);
1245 long double d
, r
, s
;
1247 d
= hypotl (re
, im
);
1248 /* Use the identity 2 Re res Im res = Im x
1249 to avoid cancellation error in d +/- Re x. */
1252 r
= sqrtl (0.5 * d
+ 0.5 * re
);
1257 s
= sqrtl (0.5 * d
- 0.5 * re
);
1258 r
= fabsl ((0.5 * im
) / s
);
1261 COMPLEX_ASSIGN (v
, r
, copysignl (s
, im
));
1268 /* sinh(a + i b) = sinh(a) cos(b) + i cosh(a) sin(b) */
1269 #if !defined(HAVE_CSINHF)
1270 #define HAVE_CSINHF 1
1271 float complex csinhf (float complex a
);
1274 csinhf (float complex a
)
1281 COMPLEX_ASSIGN (v
, sinhf (r
) * cosf (i
), coshf (r
) * sinf (i
));
1286 #if !defined(HAVE_CSINH)
1287 #define HAVE_CSINH 1
1288 double complex csinh (double complex a
);
1291 csinh (double complex a
)
1298 COMPLEX_ASSIGN (v
, sinh (r
) * cos (i
), cosh (r
) * sin (i
));
1303 #if !defined(HAVE_CSINHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1304 #define HAVE_CSINHL 1
1305 long double complex csinhl (long double complex a
);
1308 csinhl (long double complex a
)
1311 long double complex v
;
1315 COMPLEX_ASSIGN (v
, sinhl (r
) * cosl (i
), coshl (r
) * sinl (i
));
1321 /* cosh(a + i b) = cosh(a) cos(b) + i sinh(a) sin(b) */
1322 #if !defined(HAVE_CCOSHF)
1323 #define HAVE_CCOSHF 1
1324 float complex ccoshf (float complex a
);
1327 ccoshf (float complex a
)
1334 COMPLEX_ASSIGN (v
, coshf (r
) * cosf (i
), sinhf (r
) * sinf (i
));
1339 #if !defined(HAVE_CCOSH)
1340 #define HAVE_CCOSH 1
1341 double complex ccosh (double complex a
);
1344 ccosh (double complex a
)
1351 COMPLEX_ASSIGN (v
, cosh (r
) * cos (i
), sinh (r
) * sin (i
));
1356 #if !defined(HAVE_CCOSHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1357 #define HAVE_CCOSHL 1
1358 long double complex ccoshl (long double complex a
);
1361 ccoshl (long double complex a
)
1364 long double complex v
;
1368 COMPLEX_ASSIGN (v
, coshl (r
) * cosl (i
), sinhl (r
) * sinl (i
));
1374 /* tanh(a + i b) = (tanh(a) + i tan(b)) / (1 + i tanh(a) tan(b)) */
1375 #if !defined(HAVE_CTANHF)
1376 #define HAVE_CTANHF 1
1377 float complex ctanhf (float complex a
);
1380 ctanhf (float complex a
)
1385 rt
= tanhf (REALPART (a
));
1386 it
= tanf (IMAGPART (a
));
1387 COMPLEX_ASSIGN (n
, rt
, it
);
1388 COMPLEX_ASSIGN (d
, 1, rt
* it
);
1394 #if !defined(HAVE_CTANH)
1395 #define HAVE_CTANH 1
1396 double complex ctanh (double complex a
);
1398 ctanh (double complex a
)
1401 double complex n
, d
;
1403 rt
= tanh (REALPART (a
));
1404 it
= tan (IMAGPART (a
));
1405 COMPLEX_ASSIGN (n
, rt
, it
);
1406 COMPLEX_ASSIGN (d
, 1, rt
* it
);
1412 #if !defined(HAVE_CTANHL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1413 #define HAVE_CTANHL 1
1414 long double complex ctanhl (long double complex a
);
1417 ctanhl (long double complex a
)
1420 long double complex n
, d
;
1422 rt
= tanhl (REALPART (a
));
1423 it
= tanl (IMAGPART (a
));
1424 COMPLEX_ASSIGN (n
, rt
, it
);
1425 COMPLEX_ASSIGN (d
, 1, rt
* it
);
1432 /* sin(a + i b) = sin(a) cosh(b) + i cos(a) sinh(b) */
1433 #if !defined(HAVE_CSINF)
1434 #define HAVE_CSINF 1
1435 float complex csinf (float complex a
);
1438 csinf (float complex a
)
1445 COMPLEX_ASSIGN (v
, sinf (r
) * coshf (i
), cosf (r
) * sinhf (i
));
1450 #if !defined(HAVE_CSIN)
1452 double complex csin (double complex a
);
1455 csin (double complex a
)
1462 COMPLEX_ASSIGN (v
, sin (r
) * cosh (i
), cos (r
) * sinh (i
));
1467 #if !defined(HAVE_CSINL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1468 #define HAVE_CSINL 1
1469 long double complex csinl (long double complex a
);
1472 csinl (long double complex a
)
1475 long double complex v
;
1479 COMPLEX_ASSIGN (v
, sinl (r
) * coshl (i
), cosl (r
) * sinhl (i
));
1485 /* cos(a + i b) = cos(a) cosh(b) - i sin(a) sinh(b) */
1486 #if !defined(HAVE_CCOSF)
1487 #define HAVE_CCOSF 1
1488 float complex ccosf (float complex a
);
1491 ccosf (float complex a
)
1498 COMPLEX_ASSIGN (v
, cosf (r
) * coshf (i
), - (sinf (r
) * sinhf (i
)));
1503 #if !defined(HAVE_CCOS)
1505 double complex ccos (double complex a
);
1508 ccos (double complex a
)
1515 COMPLEX_ASSIGN (v
, cos (r
) * cosh (i
), - (sin (r
) * sinh (i
)));
1520 #if !defined(HAVE_CCOSL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1521 #define HAVE_CCOSL 1
1522 long double complex ccosl (long double complex a
);
1525 ccosl (long double complex a
)
1528 long double complex v
;
1532 COMPLEX_ASSIGN (v
, cosl (r
) * coshl (i
), - (sinl (r
) * sinhl (i
)));
1538 /* tan(a + i b) = (tan(a) + i tanh(b)) / (1 - i tan(a) tanh(b)) */
1539 #if !defined(HAVE_CTANF)
1540 #define HAVE_CTANF 1
1541 float complex ctanf (float complex a
);
1544 ctanf (float complex a
)
1549 rt
= tanf (REALPART (a
));
1550 it
= tanhf (IMAGPART (a
));
1551 COMPLEX_ASSIGN (n
, rt
, it
);
1552 COMPLEX_ASSIGN (d
, 1, - (rt
* it
));
1558 #if !defined(HAVE_CTAN)
1560 double complex ctan (double complex a
);
1563 ctan (double complex a
)
1566 double complex n
, d
;
1568 rt
= tan (REALPART (a
));
1569 it
= tanh (IMAGPART (a
));
1570 COMPLEX_ASSIGN (n
, rt
, it
);
1571 COMPLEX_ASSIGN (d
, 1, - (rt
* it
));
1577 #if !defined(HAVE_CTANL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1578 #define HAVE_CTANL 1
1579 long double complex ctanl (long double complex a
);
1582 ctanl (long double complex a
)
1585 long double complex n
, d
;
1587 rt
= tanl (REALPART (a
));
1588 it
= tanhl (IMAGPART (a
));
1589 COMPLEX_ASSIGN (n
, rt
, it
);
1590 COMPLEX_ASSIGN (d
, 1, - (rt
* it
));
1597 /* Complex ASIN. Returns wrongly NaN for infinite arguments.
1598 Algorithm taken from Abramowitz & Stegun. */
1600 #if !defined(HAVE_CASINF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1601 #define HAVE_CASINF 1
1602 complex float casinf (complex float z
);
1605 casinf (complex float z
)
1607 return -I
*clogf (I
*z
+ csqrtf (1.0f
-z
*z
));
1612 #if !defined(HAVE_CASIN) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1613 #define HAVE_CASIN 1
1614 complex double casin (complex double z
);
1617 casin (complex double z
)
1619 return -I
*clog (I
*z
+ csqrt (1.0-z
*z
));
1624 #if !defined(HAVE_CASINL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1625 #define HAVE_CASINL 1
1626 complex long double casinl (complex long double z
);
1629 casinl (complex long double z
)
1631 return -I
*clogl (I
*z
+ csqrtl (1.0L-z
*z
));
1636 /* Complex ACOS. Returns wrongly NaN for infinite arguments.
1637 Algorithm taken from Abramowitz & Stegun. */
1639 #if !defined(HAVE_CACOSF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1640 #define HAVE_CACOSF 1
1641 complex float cacosf (complex float z
);
1644 cacosf (complex float z
)
1646 return -I
*clogf (z
+ I
*csqrtf (1.0f
-z
*z
));
1651 #if !defined(HAVE_CACOS) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1652 #define HAVE_CACOS 1
1653 complex double cacos (complex double z
);
1656 cacos (complex double z
)
1658 return -I
*clog (z
+ I
*csqrt (1.0-z
*z
));
1663 #if !defined(HAVE_CACOSL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1664 #define HAVE_CACOSL 1
1665 complex long double cacosl (complex long double z
);
1668 cacosl (complex long double z
)
1670 return -I
*clogl (z
+ I
*csqrtl (1.0L-z
*z
));
1675 /* Complex ATAN. Returns wrongly NaN for infinite arguments.
1676 Algorithm taken from Abramowitz & Stegun. */
1678 #if !defined(HAVE_CATANF) && defined(HAVE_CLOGF)
1679 #define HAVE_CACOSF 1
1680 complex float catanf (complex float z
);
1683 catanf (complex float z
)
1685 return I
*clogf ((I
+z
)/(I
-z
))/2.0f
;
1690 #if !defined(HAVE_CATAN) && defined(HAVE_CLOG)
1691 #define HAVE_CACOS 1
1692 complex double catan (complex double z
);
1695 catan (complex double z
)
1697 return I
*clog ((I
+z
)/(I
-z
))/2.0;
1702 #if !defined(HAVE_CATANL) && defined(HAVE_CLOGL)
1703 #define HAVE_CACOSL 1
1704 complex long double catanl (complex long double z
);
1707 catanl (complex long double z
)
1709 return I
*clogl ((I
+z
)/(I
-z
))/2.0L;
1714 /* Complex ASINH. Returns wrongly NaN for infinite arguments.
1715 Algorithm taken from Abramowitz & Stegun. */
1717 #if !defined(HAVE_CASINHF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1718 #define HAVE_CASINHF 1
1719 complex float casinhf (complex float z
);
1722 casinhf (complex float z
)
1724 return clogf (z
+ csqrtf (z
*z
+1.0f
));
1729 #if !defined(HAVE_CASINH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1730 #define HAVE_CASINH 1
1731 complex double casinh (complex double z
);
1734 casinh (complex double z
)
1736 return clog (z
+ csqrt (z
*z
+1.0));
1741 #if !defined(HAVE_CASINHL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1742 #define HAVE_CASINHL 1
1743 complex long double casinhl (complex long double z
);
1746 casinhl (complex long double z
)
1748 return clogl (z
+ csqrtl (z
*z
+1.0L));
1753 /* Complex ACOSH. Returns wrongly NaN for infinite arguments.
1754 Algorithm taken from Abramowitz & Stegun. */
1756 #if !defined(HAVE_CACOSHF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1757 #define HAVE_CACOSHF 1
1758 complex float cacoshf (complex float z
);
1761 cacoshf (complex float z
)
1763 return clogf (z
+ csqrtf (z
-1.0f
) * csqrtf (z
+1.0f
));
1768 #if !defined(HAVE_CACOSH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1769 #define HAVE_CACOSH 1
1770 complex double cacosh (complex double z
);
1773 cacosh (complex double z
)
1775 return clog (z
+ csqrt (z
-1.0) * csqrt (z
+1.0));
1780 #if !defined(HAVE_CACOSHL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1781 #define HAVE_CACOSHL 1
1782 complex long double cacoshl (complex long double z
);
1785 cacoshl (complex long double z
)
1787 return clogl (z
+ csqrtl (z
-1.0L) * csqrtl (z
+1.0L));
1792 /* Complex ATANH. Returns wrongly NaN for infinite arguments.
1793 Algorithm taken from Abramowitz & Stegun. */
1795 #if !defined(HAVE_CATANHF) && defined(HAVE_CLOGF)
1796 #define HAVE_CATANHF 1
1797 complex float catanhf (complex float z
);
1800 catanhf (complex float z
)
1802 return clogf ((1.0f
+z
)/(1.0f
-z
))/2.0f
;
1807 #if !defined(HAVE_CATANH) && defined(HAVE_CLOG)
1808 #define HAVE_CATANH 1
1809 complex double catanh (complex double z
);
1812 catanh (complex double z
)
1814 return clog ((1.0+z
)/(1.0-z
))/2.0;
1818 #if !defined(HAVE_CATANHL) && defined(HAVE_CLOGL)
1819 #define HAVE_CATANHL 1
1820 complex long double catanhl (complex long double z
);
1823 catanhl (complex long double z
)
1825 return clogl ((1.0L+z
)/(1.0L-z
))/2.0L;
1830 #if !defined(HAVE_TGAMMA)
1831 #define HAVE_TGAMMA 1
1832 double tgamma (double);
1834 /* Fallback tgamma() function. Uses the algorithm from
1835 http://www.netlib.org/specfun/gamma and references therein. */
1838 #define SQRTPI 0.9189385332046727417803297
1841 #define PI 3.1415926535897932384626434
1847 double fact
, res
, sum
, xden
, xnum
, y
, y1
, ysq
, z
;
1849 static double p
[8] = {
1850 -1.71618513886549492533811e0
, 2.47656508055759199108314e1
,
1851 -3.79804256470945635097577e2
, 6.29331155312818442661052e2
,
1852 8.66966202790413211295064e2
, -3.14512729688483675254357e4
,
1853 -3.61444134186911729807069e4
, 6.64561438202405440627855e4
};
1855 static double q
[8] = {
1856 -3.08402300119738975254353e1
, 3.15350626979604161529144e2
,
1857 -1.01515636749021914166146e3
, -3.10777167157231109440444e3
,
1858 2.25381184209801510330112e4
, 4.75584627752788110767815e3
,
1859 -1.34659959864969306392456e5
, -1.15132259675553483497211e5
};
1861 static double c
[7] = { -1.910444077728e-03,
1862 8.4171387781295e-04, -5.952379913043012e-04,
1863 7.93650793500350248e-04, -2.777777777777681622553e-03,
1864 8.333333333333333331554247e-02, 5.7083835261e-03 };
1866 static const double xminin
= 2.23e-308;
1867 static const double xbig
= 171.624;
1868 static const double xnan
= __builtin_nan ("0x0"), xinf
= __builtin_inf ();
1869 static double eps
= 0;
1872 eps
= nextafter (1., 2.) - 1.;
1890 if (y1
!= trunc (y1
*0.5l)*2)
1892 fact
= -PI
/ sin (PI
*res
);
1896 return x
== 0 ? copysign (xinf
, x
) : xnan
;
1923 for (i
= 0; i
< 8; i
++)
1925 xnum
= (xnum
+ p
[i
]) * z
;
1926 xden
= xden
* z
+ q
[i
];
1929 res
= xnum
/ xden
+ 1;
1934 for (i
= 1; i
<= n
; i
++)
1946 for (i
= 0; i
< 6; i
++)
1947 sum
= sum
/ ysq
+ c
[i
];
1949 sum
= sum
/y
- y
+ SQRTPI
;
1950 sum
= sum
+ (y
- 0.5) * log (y
);
1954 return x
< 0 ? xnan
: xinf
;
1968 #if !defined(HAVE_LGAMMA)
1969 #define HAVE_LGAMMA 1
1970 double lgamma (double);
1972 /* Fallback lgamma() function. Uses the algorithm from
1973 http://www.netlib.org/specfun/algama and references therein,
1974 except for negative arguments (where netlib would return +Inf)
1975 where we use the following identity:
1976 lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y)
1984 #define SQRTPI 0.9189385332046727417803297
1987 #define PI 3.1415926535897932384626434
1989 #define PNT68 0.6796875
1990 #define D1 -0.5772156649015328605195174
1991 #define D2 0.4227843350984671393993777
1992 #define D4 1.791759469228055000094023
1994 static double p1
[8] = {
1995 4.945235359296727046734888e0
, 2.018112620856775083915565e2
,
1996 2.290838373831346393026739e3
, 1.131967205903380828685045e4
,
1997 2.855724635671635335736389e4
, 3.848496228443793359990269e4
,
1998 2.637748787624195437963534e4
, 7.225813979700288197698961e3
};
1999 static double q1
[8] = {
2000 6.748212550303777196073036e1
, 1.113332393857199323513008e3
,
2001 7.738757056935398733233834e3
, 2.763987074403340708898585e4
,
2002 5.499310206226157329794414e4
, 6.161122180066002127833352e4
,
2003 3.635127591501940507276287e4
, 8.785536302431013170870835e3
};
2004 static double p2
[8] = {
2005 4.974607845568932035012064e0
, 5.424138599891070494101986e2
,
2006 1.550693864978364947665077e4
, 1.847932904445632425417223e5
,
2007 1.088204769468828767498470e6
, 3.338152967987029735917223e6
,
2008 5.106661678927352456275255e6
, 3.074109054850539556250927e6
};
2009 static double q2
[8] = {
2010 1.830328399370592604055942e2
, 7.765049321445005871323047e3
,
2011 1.331903827966074194402448e5
, 1.136705821321969608938755e6
,
2012 5.267964117437946917577538e6
, 1.346701454311101692290052e7
,
2013 1.782736530353274213975932e7
, 9.533095591844353613395747e6
};
2014 static double p4
[8] = {
2015 1.474502166059939948905062e4
, 2.426813369486704502836312e6
,
2016 1.214755574045093227939592e8
, 2.663432449630976949898078e9
,
2017 2.940378956634553899906876e10
, 1.702665737765398868392998e11
,
2018 4.926125793377430887588120e11
, 5.606251856223951465078242e11
};
2019 static double q4
[8] = {
2020 2.690530175870899333379843e3
, 6.393885654300092398984238e5
,
2021 4.135599930241388052042842e7
, 1.120872109616147941376570e9
,
2022 1.488613728678813811542398e10
, 1.016803586272438228077304e11
,
2023 3.417476345507377132798597e11
, 4.463158187419713286462081e11
};
2024 static double c
[7] = {
2025 -1.910444077728e-03, 8.4171387781295e-04,
2026 -5.952379913043012e-04, 7.93650793500350248e-04,
2027 -2.777777777777681622553e-03, 8.333333333333333331554247e-02,
2030 static double xbig
= 2.55e305
, xinf
= __builtin_inf (), eps
= 0,
2034 double corr
, res
, xden
, xm1
, xm2
, xm4
, xnum
, ysq
;
2037 eps
= __builtin_nextafter (1., 2.) - 1.;
2039 if ((y
> 0) && (y
<= xbig
))
2053 xm1
= (y
- 0.5) - 0.5;
2056 if ((y
<= 0.5) || (y
>= PNT68
))
2060 for (i
= 0; i
< 8; i
++)
2062 xnum
= xnum
*xm1
+ p1
[i
];
2063 xden
= xden
*xm1
+ q1
[i
];
2065 res
= corr
+ (xm1
* (D1
+ xm1
*(xnum
/xden
)));
2069 xm2
= (y
- 0.5) - 0.5;
2072 for (i
= 0; i
< 8; i
++)
2074 xnum
= xnum
*xm2
+ p2
[i
];
2075 xden
= xden
*xm2
+ q2
[i
];
2077 res
= corr
+ xm2
* (D2
+ xm2
*(xnum
/xden
));
2085 for (i
= 0; i
< 8; i
++)
2087 xnum
= xnum
*xm2
+ p2
[i
];
2088 xden
= xden
*xm2
+ q2
[i
];
2090 res
= xm2
* (D2
+ xm2
*(xnum
/xden
));
2097 for (i
= 0; i
< 8; i
++)
2099 xnum
= xnum
*xm4
+ p4
[i
];
2100 xden
= xden
*xm4
+ q4
[i
];
2102 res
= D4
+ xm4
*(xnum
/xden
);
2111 for (i
= 0; i
< 6; i
++)
2112 res
= res
/ ysq
+ c
[i
];
2116 res
= res
+ SQRTPI
- 0.5*corr
;
2117 res
= res
+ y
*(corr
-1);
2120 else if (y
< 0 && __builtin_floor (y
) != y
)
2122 /* lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y)
2123 For abs(y) very close to zero, we use a series expansion to
2124 the first order in y to avoid overflow. */
2126 res
= -2 * log (fabs (y
)) - lgamma (-y
);
2128 res
= log (PI
/ fabs (y
* sin (PI
* y
))) - lgamma (-y
);
2138 #if defined(HAVE_TGAMMA) && !defined(HAVE_TGAMMAF)
2139 #define HAVE_TGAMMAF 1
2140 float tgammaf (float);
2145 return (float) tgamma ((double) x
);
2149 #if defined(HAVE_LGAMMA) && !defined(HAVE_LGAMMAF)
2150 #define HAVE_LGAMMAF 1
2151 float lgammaf (float);
2156 return (float) lgamma ((double) x
);
2162 double fma (double, double, double);
2165 fma (double x
, double y
, double z
)
2173 float fmaf (float, float, float);
2176 fmaf (float x
, float y
, float z
)
2178 return fma (x
, y
, z
);
2184 long double fmal (long double, long double, long double);
2187 fmal (long double x
, long double y
, long double z
)