1 /* Implementation of various C99 functions
2 Copyright (C) 2004, 2009, 2010 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 /* IRIX's <math.h> declares a non-C99 compliant implementation of cabs,
31 which takes two floating point arguments instead of a single complex.
32 If <complex.h> is missing this prevents building of c99_functions.c.
33 To work around this we redirect cabs{,f,l} calls to __gfc_cabs{,f,l}. */
35 #if defined(__sgi__) && !defined(HAVE_COMPLEX_H)
39 #define cabs __gfc_cabs
40 #define cabsf __gfc_cabsf
41 #define cabsl __gfc_cabsl
44 /* Tru64's <math.h> declares a non-C99 compliant implementation of cabs,
45 which takes two floating point arguments instead of a single complex.
46 To work around this we redirect cabs{,f,l} calls to __gfc_cabs{,f,l}. */
52 #define cabs __gfc_cabs
53 #define cabsf __gfc_cabsf
54 #define cabsl __gfc_cabsl
57 /* On a C99 system "I" (with I*I = -1) should be defined in complex.h;
58 if not, we define a fallback version here. */
60 # if defined(_Imaginary_I)
61 # define I _Imaginary_I
62 # elif defined(_Complex_I)
69 /* Prototypes are included to silence -Wstrict-prototypes
70 -Wmissing-prototypes. */
73 /* Wrappers for systems without the various C99 single precision Bessel
76 #if defined(HAVE_J0) && ! defined(HAVE_J0F)
83 return (float) j0 ((double) x
);
87 #if defined(HAVE_J1) && !defined(HAVE_J1F)
93 return (float) j1 ((double) x
);
97 #if defined(HAVE_JN) && !defined(HAVE_JNF)
99 float jnf (int, float);
104 return (float) jn (n
, (double) x
);
108 #if defined(HAVE_Y0) && !defined(HAVE_Y0F)
115 return (float) y0 ((double) x
);
119 #if defined(HAVE_Y1) && !defined(HAVE_Y1F)
126 return (float) y1 ((double) x
);
130 #if defined(HAVE_YN) && !defined(HAVE_YNF)
132 float ynf (int, float);
137 return (float) yn (n
, (double) x
);
142 /* Wrappers for systems without the C99 erff() and erfcf() functions. */
144 #if defined(HAVE_ERF) && !defined(HAVE_ERFF)
151 return (float) erf ((double) x
);
155 #if defined(HAVE_ERFC) && !defined(HAVE_ERFCF)
162 return (float) erfc ((double) x
);
169 float acosf (float x
);
174 return (float) acos (x
);
178 #if HAVE_ACOSH && !HAVE_ACOSHF
179 float acoshf (float x
);
184 return (float) acosh ((double) x
);
190 float asinf (float x
);
195 return (float) asin (x
);
199 #if HAVE_ASINH && !HAVE_ASINHF
200 float asinhf (float x
);
205 return (float) asinh ((double) x
);
210 #define HAVE_ATAN2F 1
211 float atan2f (float y
, float x
);
214 atan2f (float y
, float x
)
216 return (float) atan2 (y
, x
);
222 float atanf (float x
);
227 return (float) atan (x
);
231 #if HAVE_ATANH && !HAVE_ATANHF
232 float atanhf (float x
);
237 return (float) atanh ((double) x
);
243 float ceilf (float x
);
248 return (float) ceil (x
);
252 #ifndef HAVE_COPYSIGNF
253 #define HAVE_COPYSIGNF 1
254 float copysignf (float x
, float y
);
257 copysignf (float x
, float y
)
259 return (float) copysign (x
, y
);
265 float cosf (float x
);
270 return (float) cos (x
);
276 float coshf (float x
);
281 return (float) cosh (x
);
287 float expf (float x
);
292 return (float) exp (x
);
298 float fabsf (float x
);
303 return (float) fabs (x
);
308 #define HAVE_FLOORF 1
309 float floorf (float x
);
314 return (float) floor (x
);
320 float fmodf (float x
, float y
);
323 fmodf (float x
, float y
)
325 return (float) fmod (x
, y
);
330 #define HAVE_FREXPF 1
331 float frexpf (float x
, int *exp
);
334 frexpf (float x
, int *exp
)
336 return (float) frexp (x
, exp
);
341 #define HAVE_HYPOTF 1
342 float hypotf (float x
, float y
);
345 hypotf (float x
, float y
)
347 return (float) hypot (x
, y
);
353 float logf (float x
);
358 return (float) log (x
);
363 #define HAVE_LOG10F 1
364 float log10f (float x
);
369 return (float) log10 (x
);
374 #define HAVE_SCALBN 1
375 double scalbn (double x
, int y
);
378 scalbn (double x
, int y
)
380 #if (FLT_RADIX == 2) && defined(HAVE_LDEXP)
383 return x
* pow (FLT_RADIX
, y
);
389 #define HAVE_SCALBNF 1
390 float scalbnf (float x
, int y
);
393 scalbnf (float x
, int y
)
395 return (float) scalbn (x
, y
);
401 float sinf (float x
);
406 return (float) sin (x
);
412 float sinhf (float x
);
417 return (float) sinh (x
);
423 float sqrtf (float x
);
428 return (float) sqrt (x
);
434 float tanf (float x
);
439 return (float) tan (x
);
445 float tanhf (float x
);
450 return (float) tanh (x
);
456 double trunc (double x
);
472 #define HAVE_TRUNCF 1
473 float truncf (float x
);
478 return (float) trunc (x
);
482 #ifndef HAVE_NEXTAFTERF
483 #define HAVE_NEXTAFTERF 1
484 /* This is a portable implementation of nextafterf that is intended to be
485 independent of the floating point format or its in memory representation.
486 This implementation works correctly with denormalized values. */
487 float nextafterf (float x
, float y
);
490 nextafterf (float x
, float y
)
492 /* This variable is marked volatile to avoid excess precision problems
493 on some platforms, including IA-32. */
494 volatile float delta
;
495 float absx
, denorm_min
;
497 if (isnan (x
) || isnan (y
))
502 return x
> 0 ? __FLT_MAX__
: - __FLT_MAX__
;
504 /* absx = fabsf (x); */
505 absx
= (x
< 0.0) ? -x
: x
;
507 /* __FLT_DENORM_MIN__ is non-zero iff the target supports denormals. */
508 if (__FLT_DENORM_MIN__
== 0.0f
)
509 denorm_min
= __FLT_MIN__
;
511 denorm_min
= __FLT_DENORM_MIN__
;
513 if (absx
< __FLT_MIN__
)
520 /* Discard the fraction from x. */
521 frac
= frexpf (absx
, &exp
);
522 delta
= scalbnf (0.5f
, exp
);
524 /* Scale x by the epsilon of the representation. By rights we should
525 have been able to combine this with scalbnf, but some targets don't
526 get that correct with denormals. */
527 delta
*= __FLT_EPSILON__
;
529 /* If we're going to be reducing the absolute value of X, and doing so
530 would reduce the exponent of X, then the delta to be applied is
531 one exponent smaller. */
532 if (frac
== 0.5f
&& (y
< x
) == (x
> 0))
535 /* If that underflows to zero, then we're back to the minimum. */
548 #if !defined(HAVE_POWF) || defined(HAVE_BROKEN_POWF)
552 float powf (float x
, float y
);
555 powf (float x
, float y
)
557 return (float) pow (x
, y
);
562 /* Algorithm by Steven G. Kargl. */
564 #if !defined(HAVE_ROUNDL)
565 #define HAVE_ROUNDL 1
566 long double roundl (long double x
);
568 #if defined(HAVE_CEILL)
569 /* Round to nearest integral value. If the argument is halfway between two
570 integral values then round away from zero. */
573 roundl (long double x
)
596 /* Poor version of roundl for system that don't have ceill. */
598 roundl (long double x
)
600 if (x
> DBL_MAX
|| x
< -DBL_MAX
)
602 #ifdef HAVE_NEXTAFTERL
603 static long double prechalf
= nexafterl (0.5L, LDBL_MAX
);
605 static long double prechalf
= 0.5L;
607 return (GFC_INTEGER_LARGEST
) (x
+ (x
> 0 ? prechalf
: -prechalf
));
611 return round ((double) x
);
619 /* Round to nearest integral value. If the argument is halfway between two
620 integral values then round away from zero. */
621 double round (double x
);
648 #define HAVE_ROUNDF 1
649 /* Round to nearest integral value. If the argument is halfway between two
650 integral values then round away from zero. */
651 float roundf (float x
);
678 /* lround{f,,l} and llround{f,,l} functions. */
680 #if !defined(HAVE_LROUNDF) && defined(HAVE_ROUNDF)
681 #define HAVE_LROUNDF 1
682 long int lroundf (float x
);
687 return (long int) roundf (x
);
691 #if !defined(HAVE_LROUND) && defined(HAVE_ROUND)
692 #define HAVE_LROUND 1
693 long int lround (double x
);
698 return (long int) round (x
);
702 #if !defined(HAVE_LROUNDL) && defined(HAVE_ROUNDL)
703 #define HAVE_LROUNDL 1
704 long int lroundl (long double x
);
707 lroundl (long double x
)
709 return (long long int) roundl (x
);
713 #if !defined(HAVE_LLROUNDF) && defined(HAVE_ROUNDF)
714 #define HAVE_LLROUNDF 1
715 long long int llroundf (float x
);
720 return (long long int) roundf (x
);
724 #if !defined(HAVE_LLROUND) && defined(HAVE_ROUND)
725 #define HAVE_LLROUND 1
726 long long int llround (double x
);
731 return (long long int) round (x
);
735 #if !defined(HAVE_LLROUNDL) && defined(HAVE_ROUNDL)
736 #define HAVE_LLROUNDL 1
737 long long int llroundl (long double x
);
740 llroundl (long double x
)
742 return (long long int) roundl (x
);
748 #define HAVE_LOG10L 1
749 /* log10 function for long double variables. The version provided here
750 reduces the argument until it fits into a double, then use log10. */
751 long double log10l (long double x
);
754 log10l (long double x
)
756 #if LDBL_MAX_EXP > DBL_MAX_EXP
761 if (x
> 0x1p
16383L) { p2_result
+= 16383; x
/= 0x1p
16383L; }
762 if (x
> 0x1p
8191L) { p2_result
+= 8191; x
/= 0x1p
8191L; }
763 if (x
> 0x1p
4095L) { p2_result
+= 4095; x
/= 0x1p
4095L; }
764 if (x
> 0x1p
2047L) { p2_result
+= 2047; x
/= 0x1p
2047L; }
765 if (x
> 0x1p
1023L) { p2_result
+= 1023; x
/= 0x1p
1023L; }
766 val
= log10 ((double) x
);
767 return (val
+ p2_result
* .30102999566398119521373889472449302L);
770 #if LDBL_MIN_EXP < DBL_MIN_EXP
775 if (x
< 0x1p
-16380L) { p2_result
+= 16380; x
/= 0x1p
-16380L; }
776 if (x
< 0x1p
-8189L) { p2_result
+= 8189; x
/= 0x1p
-8189L; }
777 if (x
< 0x1p
-4093L) { p2_result
+= 4093; x
/= 0x1p
-4093L; }
778 if (x
< 0x1p
-2045L) { p2_result
+= 2045; x
/= 0x1p
-2045L; }
779 if (x
< 0x1p
-1021L) { p2_result
+= 1021; x
/= 0x1p
-1021L; }
780 val
= fabs (log10 ((double) x
));
781 return (- val
- p2_result
* .30102999566398119521373889472449302L);
790 #define HAVE_FLOORL 1
791 long double floorl (long double x
);
794 floorl (long double x
)
796 /* Zero, possibly signed. */
800 /* Large magnitude. */
801 if (x
> DBL_MAX
|| x
< (-DBL_MAX
))
804 /* Small positive values. */
805 if (x
>= 0 && x
< DBL_MIN
)
808 /* Small negative values. */
809 if (x
< 0 && x
> (-DBL_MIN
))
819 long double fmodl (long double x
, long double y
);
822 fmodl (long double x
, long double y
)
827 /* Need to check that the result has the same sign as x and magnitude
828 less than the magnitude of y. */
829 return x
- floorl (x
/ y
) * y
;
834 #if !defined(HAVE_CABSF)
836 float cabsf (float complex z
);
839 cabsf (float complex z
)
841 return hypotf (REALPART (z
), IMAGPART (z
));
845 #if !defined(HAVE_CABS)
847 double cabs (double complex z
);
850 cabs (double complex z
)
852 return hypot (REALPART (z
), IMAGPART (z
));
856 #if !defined(HAVE_CABSL) && defined(HAVE_HYPOTL)
858 long double cabsl (long double complex z
);
861 cabsl (long double complex z
)
863 return hypotl (REALPART (z
), IMAGPART (z
));
868 #if !defined(HAVE_CARGF)
870 float cargf (float complex z
);
873 cargf (float complex z
)
875 return atan2f (IMAGPART (z
), REALPART (z
));
879 #if !defined(HAVE_CARG)
881 double carg (double complex z
);
884 carg (double complex z
)
886 return atan2 (IMAGPART (z
), REALPART (z
));
890 #if !defined(HAVE_CARGL) && defined(HAVE_ATAN2L)
892 long double cargl (long double complex z
);
895 cargl (long double complex z
)
897 return atan2l (IMAGPART (z
), REALPART (z
));
902 /* exp(z) = exp(a)*(cos(b) + i sin(b)) */
903 #if !defined(HAVE_CEXPF)
905 float complex cexpf (float complex z
);
908 cexpf (float complex z
)
915 COMPLEX_ASSIGN (v
, cosf (b
), sinf (b
));
920 #if !defined(HAVE_CEXP)
922 double complex cexp (double complex z
);
925 cexp (double complex z
)
932 COMPLEX_ASSIGN (v
, cos (b
), sin (b
));
937 #if !defined(HAVE_CEXPL) && defined(HAVE_COSL) && defined(HAVE_SINL) && defined(EXPL)
939 long double complex cexpl (long double complex z
);
942 cexpl (long double complex z
)
945 long double complex v
;
949 COMPLEX_ASSIGN (v
, cosl (b
), sinl (b
));
955 /* log(z) = log (cabs(z)) + i*carg(z) */
956 #if !defined(HAVE_CLOGF)
958 float complex clogf (float complex z
);
961 clogf (float complex z
)
965 COMPLEX_ASSIGN (v
, logf (cabsf (z
)), cargf (z
));
970 #if !defined(HAVE_CLOG)
972 double complex clog (double complex z
);
975 clog (double complex z
)
979 COMPLEX_ASSIGN (v
, log (cabs (z
)), carg (z
));
984 #if !defined(HAVE_CLOGL) && defined(HAVE_LOGL) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
986 long double complex clogl (long double complex z
);
989 clogl (long double complex z
)
991 long double complex v
;
993 COMPLEX_ASSIGN (v
, logl (cabsl (z
)), cargl (z
));
999 /* log10(z) = log10 (cabs(z)) + i*carg(z) */
1000 #if !defined(HAVE_CLOG10F)
1001 #define HAVE_CLOG10F 1
1002 float complex clog10f (float complex z
);
1005 clog10f (float complex z
)
1009 COMPLEX_ASSIGN (v
, log10f (cabsf (z
)), cargf (z
));
1014 #if !defined(HAVE_CLOG10)
1015 #define HAVE_CLOG10 1
1016 double complex clog10 (double complex z
);
1019 clog10 (double complex z
)
1023 COMPLEX_ASSIGN (v
, log10 (cabs (z
)), carg (z
));
1028 #if !defined(HAVE_CLOG10L) && defined(HAVE_LOG10L) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
1029 #define HAVE_CLOG10L 1
1030 long double complex clog10l (long double complex z
);
1033 clog10l (long double complex z
)
1035 long double complex v
;
1037 COMPLEX_ASSIGN (v
, log10l (cabsl (z
)), cargl (z
));
1043 /* pow(base, power) = cexp (power * clog (base)) */
1044 #if !defined(HAVE_CPOWF)
1045 #define HAVE_CPOWF 1
1046 float complex cpowf (float complex base
, float complex power
);
1049 cpowf (float complex base
, float complex power
)
1051 return cexpf (power
* clogf (base
));
1055 #if !defined(HAVE_CPOW)
1057 double complex cpow (double complex base
, double complex power
);
1060 cpow (double complex base
, double complex power
)
1062 return cexp (power
* clog (base
));
1066 #if !defined(HAVE_CPOWL) && defined(HAVE_CEXPL) && defined(HAVE_CLOGL)
1067 #define HAVE_CPOWL 1
1068 long double complex cpowl (long double complex base
, long double complex power
);
1071 cpowl (long double complex base
, long double complex power
)
1073 return cexpl (power
* clogl (base
));
1078 /* sqrt(z). Algorithm pulled from glibc. */
1079 #if !defined(HAVE_CSQRTF)
1080 #define HAVE_CSQRTF 1
1081 float complex csqrtf (float complex z
);
1084 csqrtf (float complex z
)
1095 COMPLEX_ASSIGN (v
, 0, copysignf (sqrtf (-re
), im
));
1099 COMPLEX_ASSIGN (v
, fabsf (sqrtf (re
)), copysignf (0, im
));
1106 r
= sqrtf (0.5 * fabsf (im
));
1108 COMPLEX_ASSIGN (v
, r
, copysignf (r
, im
));
1114 d
= hypotf (re
, im
);
1115 /* Use the identity 2 Re res Im res = Im x
1116 to avoid cancellation error in d +/- Re x. */
1119 r
= sqrtf (0.5 * d
+ 0.5 * re
);
1124 s
= sqrtf (0.5 * d
- 0.5 * re
);
1125 r
= fabsf ((0.5 * im
) / s
);
1128 COMPLEX_ASSIGN (v
, r
, copysignf (s
, im
));
1134 #if !defined(HAVE_CSQRT)
1135 #define HAVE_CSQRT 1
1136 double complex csqrt (double complex z
);
1139 csqrt (double complex z
)
1150 COMPLEX_ASSIGN (v
, 0, copysign (sqrt (-re
), im
));
1154 COMPLEX_ASSIGN (v
, fabs (sqrt (re
)), copysign (0, im
));
1161 r
= sqrt (0.5 * fabs (im
));
1163 COMPLEX_ASSIGN (v
, r
, copysign (r
, im
));
1170 /* Use the identity 2 Re res Im res = Im x
1171 to avoid cancellation error in d +/- Re x. */
1174 r
= sqrt (0.5 * d
+ 0.5 * re
);
1179 s
= sqrt (0.5 * d
- 0.5 * re
);
1180 r
= fabs ((0.5 * im
) / s
);
1183 COMPLEX_ASSIGN (v
, r
, copysign (s
, im
));
1189 #if !defined(HAVE_CSQRTL) && defined(HAVE_COPYSIGNL) && defined(HAVE_SQRTL) && defined(HAVE_FABSL) && defined(HAVE_HYPOTL)
1190 #define HAVE_CSQRTL 1
1191 long double complex csqrtl (long double complex z
);
1194 csqrtl (long double complex z
)
1197 long double complex v
;
1205 COMPLEX_ASSIGN (v
, 0, copysignl (sqrtl (-re
), im
));
1209 COMPLEX_ASSIGN (v
, fabsl (sqrtl (re
)), copysignl (0, im
));
1216 r
= sqrtl (0.5 * fabsl (im
));
1218 COMPLEX_ASSIGN (v
, copysignl (r
, im
), r
);
1222 long double d
, r
, s
;
1224 d
= hypotl (re
, im
);
1225 /* Use the identity 2 Re res Im res = Im x
1226 to avoid cancellation error in d +/- Re x. */
1229 r
= sqrtl (0.5 * d
+ 0.5 * re
);
1234 s
= sqrtl (0.5 * d
- 0.5 * re
);
1235 r
= fabsl ((0.5 * im
) / s
);
1238 COMPLEX_ASSIGN (v
, r
, copysignl (s
, im
));
1245 /* sinh(a + i b) = sinh(a) cos(b) + i cosh(a) sin(b) */
1246 #if !defined(HAVE_CSINHF)
1247 #define HAVE_CSINHF 1
1248 float complex csinhf (float complex a
);
1251 csinhf (float complex a
)
1258 COMPLEX_ASSIGN (v
, sinhf (r
) * cosf (i
), coshf (r
) * sinf (i
));
1263 #if !defined(HAVE_CSINH)
1264 #define HAVE_CSINH 1
1265 double complex csinh (double complex a
);
1268 csinh (double complex a
)
1275 COMPLEX_ASSIGN (v
, sinh (r
) * cos (i
), cosh (r
) * sin (i
));
1280 #if !defined(HAVE_CSINHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1281 #define HAVE_CSINHL 1
1282 long double complex csinhl (long double complex a
);
1285 csinhl (long double complex a
)
1288 long double complex v
;
1292 COMPLEX_ASSIGN (v
, sinhl (r
) * cosl (i
), coshl (r
) * sinl (i
));
1298 /* cosh(a + i b) = cosh(a) cos(b) + i sinh(a) sin(b) */
1299 #if !defined(HAVE_CCOSHF)
1300 #define HAVE_CCOSHF 1
1301 float complex ccoshf (float complex a
);
1304 ccoshf (float complex a
)
1311 COMPLEX_ASSIGN (v
, coshf (r
) * cosf (i
), sinhf (r
) * sinf (i
));
1316 #if !defined(HAVE_CCOSH)
1317 #define HAVE_CCOSH 1
1318 double complex ccosh (double complex a
);
1321 ccosh (double complex a
)
1328 COMPLEX_ASSIGN (v
, cosh (r
) * cos (i
), sinh (r
) * sin (i
));
1333 #if !defined(HAVE_CCOSHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1334 #define HAVE_CCOSHL 1
1335 long double complex ccoshl (long double complex a
);
1338 ccoshl (long double complex a
)
1341 long double complex v
;
1345 COMPLEX_ASSIGN (v
, coshl (r
) * cosl (i
), sinhl (r
) * sinl (i
));
1351 /* tanh(a + i b) = (tanh(a) + i tan(b)) / (1 + i tanh(a) tan(b)) */
1352 #if !defined(HAVE_CTANHF)
1353 #define HAVE_CTANHF 1
1354 float complex ctanhf (float complex a
);
1357 ctanhf (float complex a
)
1362 rt
= tanhf (REALPART (a
));
1363 it
= tanf (IMAGPART (a
));
1364 COMPLEX_ASSIGN (n
, rt
, it
);
1365 COMPLEX_ASSIGN (d
, 1, rt
* it
);
1371 #if !defined(HAVE_CTANH)
1372 #define HAVE_CTANH 1
1373 double complex ctanh (double complex a
);
1375 ctanh (double complex a
)
1378 double complex n
, d
;
1380 rt
= tanh (REALPART (a
));
1381 it
= tan (IMAGPART (a
));
1382 COMPLEX_ASSIGN (n
, rt
, it
);
1383 COMPLEX_ASSIGN (d
, 1, rt
* it
);
1389 #if !defined(HAVE_CTANHL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1390 #define HAVE_CTANHL 1
1391 long double complex ctanhl (long double complex a
);
1394 ctanhl (long double complex a
)
1397 long double complex n
, d
;
1399 rt
= tanhl (REALPART (a
));
1400 it
= tanl (IMAGPART (a
));
1401 COMPLEX_ASSIGN (n
, rt
, it
);
1402 COMPLEX_ASSIGN (d
, 1, rt
* it
);
1409 /* sin(a + i b) = sin(a) cosh(b) + i cos(a) sinh(b) */
1410 #if !defined(HAVE_CSINF)
1411 #define HAVE_CSINF 1
1412 float complex csinf (float complex a
);
1415 csinf (float complex a
)
1422 COMPLEX_ASSIGN (v
, sinf (r
) * coshf (i
), cosf (r
) * sinhf (i
));
1427 #if !defined(HAVE_CSIN)
1429 double complex csin (double complex a
);
1432 csin (double complex a
)
1439 COMPLEX_ASSIGN (v
, sin (r
) * cosh (i
), cos (r
) * sinh (i
));
1444 #if !defined(HAVE_CSINL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1445 #define HAVE_CSINL 1
1446 long double complex csinl (long double complex a
);
1449 csinl (long double complex a
)
1452 long double complex v
;
1456 COMPLEX_ASSIGN (v
, sinl (r
) * coshl (i
), cosl (r
) * sinhl (i
));
1462 /* cos(a + i b) = cos(a) cosh(b) - i sin(a) sinh(b) */
1463 #if !defined(HAVE_CCOSF)
1464 #define HAVE_CCOSF 1
1465 float complex ccosf (float complex a
);
1468 ccosf (float complex a
)
1475 COMPLEX_ASSIGN (v
, cosf (r
) * coshf (i
), - (sinf (r
) * sinhf (i
)));
1480 #if !defined(HAVE_CCOS)
1482 double complex ccos (double complex a
);
1485 ccos (double complex a
)
1492 COMPLEX_ASSIGN (v
, cos (r
) * cosh (i
), - (sin (r
) * sinh (i
)));
1497 #if !defined(HAVE_CCOSL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1498 #define HAVE_CCOSL 1
1499 long double complex ccosl (long double complex a
);
1502 ccosl (long double complex a
)
1505 long double complex v
;
1509 COMPLEX_ASSIGN (v
, cosl (r
) * coshl (i
), - (sinl (r
) * sinhl (i
)));
1515 /* tan(a + i b) = (tan(a) + i tanh(b)) / (1 - i tan(a) tanh(b)) */
1516 #if !defined(HAVE_CTANF)
1517 #define HAVE_CTANF 1
1518 float complex ctanf (float complex a
);
1521 ctanf (float complex a
)
1526 rt
= tanf (REALPART (a
));
1527 it
= tanhf (IMAGPART (a
));
1528 COMPLEX_ASSIGN (n
, rt
, it
);
1529 COMPLEX_ASSIGN (d
, 1, - (rt
* it
));
1535 #if !defined(HAVE_CTAN)
1537 double complex ctan (double complex a
);
1540 ctan (double complex a
)
1543 double complex n
, d
;
1545 rt
= tan (REALPART (a
));
1546 it
= tanh (IMAGPART (a
));
1547 COMPLEX_ASSIGN (n
, rt
, it
);
1548 COMPLEX_ASSIGN (d
, 1, - (rt
* it
));
1554 #if !defined(HAVE_CTANL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1555 #define HAVE_CTANL 1
1556 long double complex ctanl (long double complex a
);
1559 ctanl (long double complex a
)
1562 long double complex n
, d
;
1564 rt
= tanl (REALPART (a
));
1565 it
= tanhl (IMAGPART (a
));
1566 COMPLEX_ASSIGN (n
, rt
, it
);
1567 COMPLEX_ASSIGN (d
, 1, - (rt
* it
));
1574 /* Complex ASIN. Returns wrongly NaN for infinite arguments.
1575 Algorithm taken from Abramowitz & Stegun. */
1577 #if !defined(HAVE_CASINF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1578 #define HAVE_CASINF 1
1579 complex float casinf (complex float z
);
1582 casinf (complex float z
)
1584 return -I
*clogf (I
*z
+ csqrtf (1.0f
-z
*z
));
1589 #if !defined(HAVE_CASIN) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1590 #define HAVE_CASIN 1
1591 complex double casin (complex double z
);
1594 casin (complex double z
)
1596 return -I
*clog (I
*z
+ csqrt (1.0-z
*z
));
1601 #if !defined(HAVE_CASINL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1602 #define HAVE_CASINL 1
1603 complex long double casinl (complex long double z
);
1606 casinl (complex long double z
)
1608 return -I
*clogl (I
*z
+ csqrtl (1.0L-z
*z
));
1613 /* Complex ACOS. Returns wrongly NaN for infinite arguments.
1614 Algorithm taken from Abramowitz & Stegun. */
1616 #if !defined(HAVE_CACOSF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1617 #define HAVE_CACOSF 1
1618 complex float cacosf (complex float z
);
1621 cacosf (complex float z
)
1623 return -I
*clogf (z
+ I
*csqrtf (1.0f
-z
*z
));
1628 #if !defined(HAVE_CACOS) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1629 #define HAVE_CACOS 1
1630 complex double cacos (complex double z
);
1633 cacos (complex double z
)
1635 return -I
*clog (z
+ I
*csqrt (1.0-z
*z
));
1640 #if !defined(HAVE_CACOSL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1641 #define HAVE_CACOSL 1
1642 complex long double cacosl (complex long double z
);
1645 cacosl (complex long double z
)
1647 return -I
*clogl (z
+ I
*csqrtl (1.0L-z
*z
));
1652 /* Complex ATAN. Returns wrongly NaN for infinite arguments.
1653 Algorithm taken from Abramowitz & Stegun. */
1655 #if !defined(HAVE_CATANF) && defined(HAVE_CLOGF)
1656 #define HAVE_CACOSF 1
1657 complex float catanf (complex float z
);
1660 catanf (complex float z
)
1662 return I
*clogf ((I
+z
)/(I
-z
))/2.0f
;
1667 #if !defined(HAVE_CATAN) && defined(HAVE_CLOG)
1668 #define HAVE_CACOS 1
1669 complex double catan (complex double z
);
1672 catan (complex double z
)
1674 return I
*clog ((I
+z
)/(I
-z
))/2.0;
1679 #if !defined(HAVE_CATANL) && defined(HAVE_CLOGL)
1680 #define HAVE_CACOSL 1
1681 complex long double catanl (complex long double z
);
1684 catanl (complex long double z
)
1686 return I
*clogl ((I
+z
)/(I
-z
))/2.0L;
1691 /* Complex ASINH. Returns wrongly NaN for infinite arguments.
1692 Algorithm taken from Abramowitz & Stegun. */
1694 #if !defined(HAVE_CASINHF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1695 #define HAVE_CASINHF 1
1696 complex float casinhf (complex float z
);
1699 casinhf (complex float z
)
1701 return clogf (z
+ csqrtf (z
*z
+1.0f
));
1706 #if !defined(HAVE_CASINH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1707 #define HAVE_CASINH 1
1708 complex double casinh (complex double z
);
1711 casinh (complex double z
)
1713 return clog (z
+ csqrt (z
*z
+1.0));
1718 #if !defined(HAVE_CASINHL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1719 #define HAVE_CASINHL 1
1720 complex long double casinhl (complex long double z
);
1723 casinhl (complex long double z
)
1725 return clogl (z
+ csqrtl (z
*z
+1.0L));
1730 /* Complex ACOSH. Returns wrongly NaN for infinite arguments.
1731 Algorithm taken from Abramowitz & Stegun. */
1733 #if !defined(HAVE_CACOSHF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1734 #define HAVE_CACOSHF 1
1735 complex float cacoshf (complex float z
);
1738 cacoshf (complex float z
)
1740 return clogf (z
+ csqrtf (z
-1.0f
) * csqrtf (z
+1.0f
));
1745 #if !defined(HAVE_CACOSH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1746 #define HAVE_CACOSH 1
1747 complex double cacosh (complex double z
);
1750 cacosh (complex double z
)
1752 return clog (z
+ csqrt (z
-1.0) * csqrt (z
+1.0));
1757 #if !defined(HAVE_CACOSHL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1758 #define HAVE_CACOSHL 1
1759 complex long double cacoshl (complex long double z
);
1762 cacoshl (complex long double z
)
1764 return clogl (z
+ csqrtl (z
-1.0L) * csqrtl (z
+1.0L));
1769 /* Complex ATANH. Returns wrongly NaN for infinite arguments.
1770 Algorithm taken from Abramowitz & Stegun. */
1772 #if !defined(HAVE_CATANHF) && defined(HAVE_CLOGF)
1773 #define HAVE_CATANHF 1
1774 complex float catanhf (complex float z
);
1777 catanhf (complex float z
)
1779 return clogf ((1.0f
+z
)/(1.0f
-z
))/2.0f
;
1784 #if !defined(HAVE_CATANH) && defined(HAVE_CLOG)
1785 #define HAVE_CATANH 1
1786 complex double catanh (complex double z
);
1789 catanh (complex double z
)
1791 return clog ((1.0+z
)/(1.0-z
))/2.0;
1795 #if !defined(HAVE_CATANHL) && defined(HAVE_CLOGL)
1796 #define HAVE_CATANHL 1
1797 complex long double catanhl (complex long double z
);
1800 catanhl (complex long double z
)
1802 return clogl ((1.0L+z
)/(1.0L-z
))/2.0L;
1807 #if !defined(HAVE_TGAMMA)
1808 #define HAVE_TGAMMA 1
1809 double tgamma (double);
1811 /* Fallback tgamma() function. Uses the algorithm from
1812 http://www.netlib.org/specfun/gamma and references therein. */
1815 #define SQRTPI 0.9189385332046727417803297
1818 #define PI 3.1415926535897932384626434
1824 double fact
, res
, sum
, xden
, xnum
, y
, y1
, ysq
, z
;
1826 static double p
[8] = {
1827 -1.71618513886549492533811e0
, 2.47656508055759199108314e1
,
1828 -3.79804256470945635097577e2
, 6.29331155312818442661052e2
,
1829 8.66966202790413211295064e2
, -3.14512729688483675254357e4
,
1830 -3.61444134186911729807069e4
, 6.64561438202405440627855e4
};
1832 static double q
[8] = {
1833 -3.08402300119738975254353e1
, 3.15350626979604161529144e2
,
1834 -1.01515636749021914166146e3
, -3.10777167157231109440444e3
,
1835 2.25381184209801510330112e4
, 4.75584627752788110767815e3
,
1836 -1.34659959864969306392456e5
, -1.15132259675553483497211e5
};
1838 static double c
[7] = { -1.910444077728e-03,
1839 8.4171387781295e-04, -5.952379913043012e-04,
1840 7.93650793500350248e-04, -2.777777777777681622553e-03,
1841 8.333333333333333331554247e-02, 5.7083835261e-03 };
1843 static const double xminin
= 2.23e-308;
1844 static const double xbig
= 171.624;
1845 static const double xnan
= __builtin_nan ("0x0"), xinf
= __builtin_inf ();
1846 static double eps
= 0;
1849 eps
= nextafter (1., 2.) - 1.;
1867 if (y1
!= trunc (y1
*0.5l)*2)
1869 fact
= -PI
/ sin (PI
*res
);
1873 return x
== 0 ? copysign (xinf
, x
) : xnan
;
1900 for (i
= 0; i
< 8; i
++)
1902 xnum
= (xnum
+ p
[i
]) * z
;
1903 xden
= xden
* z
+ q
[i
];
1906 res
= xnum
/ xden
+ 1;
1911 for (i
= 1; i
<= n
; i
++)
1923 for (i
= 0; i
< 6; i
++)
1924 sum
= sum
/ ysq
+ c
[i
];
1926 sum
= sum
/y
- y
+ SQRTPI
;
1927 sum
= sum
+ (y
- 0.5) * log (y
);
1931 return x
< 0 ? xnan
: xinf
;
1945 #if !defined(HAVE_LGAMMA)
1946 #define HAVE_LGAMMA 1
1947 double lgamma (double);
1949 /* Fallback lgamma() function. Uses the algorithm from
1950 http://www.netlib.org/specfun/algama and references therein,
1951 except for negative arguments (where netlib would return +Inf)
1952 where we use the following identity:
1953 lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y)
1961 #define SQRTPI 0.9189385332046727417803297
1964 #define PI 3.1415926535897932384626434
1966 #define PNT68 0.6796875
1967 #define D1 -0.5772156649015328605195174
1968 #define D2 0.4227843350984671393993777
1969 #define D4 1.791759469228055000094023
1971 static double p1
[8] = {
1972 4.945235359296727046734888e0
, 2.018112620856775083915565e2
,
1973 2.290838373831346393026739e3
, 1.131967205903380828685045e4
,
1974 2.855724635671635335736389e4
, 3.848496228443793359990269e4
,
1975 2.637748787624195437963534e4
, 7.225813979700288197698961e3
};
1976 static double q1
[8] = {
1977 6.748212550303777196073036e1
, 1.113332393857199323513008e3
,
1978 7.738757056935398733233834e3
, 2.763987074403340708898585e4
,
1979 5.499310206226157329794414e4
, 6.161122180066002127833352e4
,
1980 3.635127591501940507276287e4
, 8.785536302431013170870835e3
};
1981 static double p2
[8] = {
1982 4.974607845568932035012064e0
, 5.424138599891070494101986e2
,
1983 1.550693864978364947665077e4
, 1.847932904445632425417223e5
,
1984 1.088204769468828767498470e6
, 3.338152967987029735917223e6
,
1985 5.106661678927352456275255e6
, 3.074109054850539556250927e6
};
1986 static double q2
[8] = {
1987 1.830328399370592604055942e2
, 7.765049321445005871323047e3
,
1988 1.331903827966074194402448e5
, 1.136705821321969608938755e6
,
1989 5.267964117437946917577538e6
, 1.346701454311101692290052e7
,
1990 1.782736530353274213975932e7
, 9.533095591844353613395747e6
};
1991 static double p4
[8] = {
1992 1.474502166059939948905062e4
, 2.426813369486704502836312e6
,
1993 1.214755574045093227939592e8
, 2.663432449630976949898078e9
,
1994 2.940378956634553899906876e10
, 1.702665737765398868392998e11
,
1995 4.926125793377430887588120e11
, 5.606251856223951465078242e11
};
1996 static double q4
[8] = {
1997 2.690530175870899333379843e3
, 6.393885654300092398984238e5
,
1998 4.135599930241388052042842e7
, 1.120872109616147941376570e9
,
1999 1.488613728678813811542398e10
, 1.016803586272438228077304e11
,
2000 3.417476345507377132798597e11
, 4.463158187419713286462081e11
};
2001 static double c
[7] = {
2002 -1.910444077728e-03, 8.4171387781295e-04,
2003 -5.952379913043012e-04, 7.93650793500350248e-04,
2004 -2.777777777777681622553e-03, 8.333333333333333331554247e-02,
2007 static double xbig
= 2.55e305
, xinf
= __builtin_inf (), eps
= 0,
2011 double corr
, res
, xden
, xm1
, xm2
, xm4
, xnum
, ysq
;
2014 eps
= __builtin_nextafter (1., 2.) - 1.;
2016 if ((y
> 0) && (y
<= xbig
))
2030 xm1
= (y
- 0.5) - 0.5;
2033 if ((y
<= 0.5) || (y
>= PNT68
))
2037 for (i
= 0; i
< 8; i
++)
2039 xnum
= xnum
*xm1
+ p1
[i
];
2040 xden
= xden
*xm1
+ q1
[i
];
2042 res
= corr
+ (xm1
* (D1
+ xm1
*(xnum
/xden
)));
2046 xm2
= (y
- 0.5) - 0.5;
2049 for (i
= 0; i
< 8; i
++)
2051 xnum
= xnum
*xm2
+ p2
[i
];
2052 xden
= xden
*xm2
+ q2
[i
];
2054 res
= corr
+ xm2
* (D2
+ xm2
*(xnum
/xden
));
2062 for (i
= 0; i
< 8; i
++)
2064 xnum
= xnum
*xm2
+ p2
[i
];
2065 xden
= xden
*xm2
+ q2
[i
];
2067 res
= xm2
* (D2
+ xm2
*(xnum
/xden
));
2074 for (i
= 0; i
< 8; i
++)
2076 xnum
= xnum
*xm4
+ p4
[i
];
2077 xden
= xden
*xm4
+ q4
[i
];
2079 res
= D4
+ xm4
*(xnum
/xden
);
2088 for (i
= 0; i
< 6; i
++)
2089 res
= res
/ ysq
+ c
[i
];
2093 res
= res
+ SQRTPI
- 0.5*corr
;
2094 res
= res
+ y
*(corr
-1);
2097 else if (y
< 0 && __builtin_floor (y
) != y
)
2099 /* lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y)
2100 For abs(y) very close to zero, we use a series expansion to
2101 the first order in y to avoid overflow. */
2103 res
= -2 * log (fabs (y
)) - lgamma (-y
);
2105 res
= log (PI
/ fabs (y
* sin (PI
* y
))) - lgamma (-y
);
2115 #if defined(HAVE_TGAMMA) && !defined(HAVE_TGAMMAF)
2116 #define HAVE_TGAMMAF 1
2117 float tgammaf (float);
2122 return (float) tgamma ((double) x
);
2126 #if defined(HAVE_LGAMMA) && !defined(HAVE_LGAMMAF)
2127 #define HAVE_LGAMMAF 1
2128 float lgammaf (float);
2133 return (float) lgamma ((double) x
);