make __stl_prime_list in comdat
[official-gcc.git] / libgfortran / intrinsics / c99_functions.c
blob9ba5544a02b7ea5ea2fb6f3ea301e156b81d4de1
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/>. */
25 #include "config.h"
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)
36 #undef HAVE_CABS
37 #undef HAVE_CABSF
38 #undef HAVE_CABSL
39 #define cabs __gfc_cabs
40 #define cabsf __gfc_cabsf
41 #define cabsl __gfc_cabsl
42 #endif
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}. */
48 #ifdef __osf__
49 #undef HAVE_CABS
50 #undef HAVE_CABSF
51 #undef HAVE_CABSL
52 #define cabs __gfc_cabs
53 #define cabsf __gfc_cabsf
54 #define cabsl __gfc_cabsl
55 #endif
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. */
59 #ifndef I
60 # if defined(_Imaginary_I)
61 # define I _Imaginary_I
62 # elif defined(_Complex_I)
63 # define I _Complex_I
64 # else
65 # define I (1.0fi)
66 # endif
67 #endif
69 /* Prototypes are included to silence -Wstrict-prototypes
70 -Wmissing-prototypes. */
73 /* Wrappers for systems without the various C99 single precision Bessel
74 functions. */
76 #if defined(HAVE_J0) && ! defined(HAVE_J0F)
77 #define HAVE_J0F 1
78 float j0f (float);
80 float
81 j0f (float x)
83 return (float) j0 ((double) x);
85 #endif
87 #if defined(HAVE_J1) && !defined(HAVE_J1F)
88 #define HAVE_J1F 1
89 float j1f (float);
91 float j1f (float x)
93 return (float) j1 ((double) x);
95 #endif
97 #if defined(HAVE_JN) && !defined(HAVE_JNF)
98 #define HAVE_JNF 1
99 float jnf (int, float);
101 float
102 jnf (int n, float x)
104 return (float) jn (n, (double) x);
106 #endif
108 #if defined(HAVE_Y0) && !defined(HAVE_Y0F)
109 #define HAVE_Y0F 1
110 float y0f (float);
112 float
113 y0f (float x)
115 return (float) y0 ((double) x);
117 #endif
119 #if defined(HAVE_Y1) && !defined(HAVE_Y1F)
120 #define HAVE_Y1F 1
121 float y1f (float);
123 float
124 y1f (float x)
126 return (float) y1 ((double) x);
128 #endif
130 #if defined(HAVE_YN) && !defined(HAVE_YNF)
131 #define HAVE_YNF 1
132 float ynf (int, float);
134 float
135 ynf (int n, float x)
137 return (float) yn (n, (double) x);
139 #endif
142 /* Wrappers for systems without the C99 erff() and erfcf() functions. */
144 #if defined(HAVE_ERF) && !defined(HAVE_ERFF)
145 #define HAVE_ERFF 1
146 float erff (float);
148 float
149 erff (float x)
151 return (float) erf ((double) x);
153 #endif
155 #if defined(HAVE_ERFC) && !defined(HAVE_ERFCF)
156 #define HAVE_ERFCF 1
157 float erfcf (float);
159 float
160 erfcf (float x)
162 return (float) erfc ((double) x);
164 #endif
167 #ifndef HAVE_ACOSF
168 #define HAVE_ACOSF 1
169 float acosf (float x);
171 float
172 acosf (float x)
174 return (float) acos (x);
176 #endif
178 #if HAVE_ACOSH && !HAVE_ACOSHF
179 float acoshf (float x);
181 float
182 acoshf (float x)
184 return (float) acosh ((double) x);
186 #endif
188 #ifndef HAVE_ASINF
189 #define HAVE_ASINF 1
190 float asinf (float x);
192 float
193 asinf (float x)
195 return (float) asin (x);
197 #endif
199 #if HAVE_ASINH && !HAVE_ASINHF
200 float asinhf (float x);
202 float
203 asinhf (float x)
205 return (float) asinh ((double) x);
207 #endif
209 #ifndef HAVE_ATAN2F
210 #define HAVE_ATAN2F 1
211 float atan2f (float y, float x);
213 float
214 atan2f (float y, float x)
216 return (float) atan2 (y, x);
218 #endif
220 #ifndef HAVE_ATANF
221 #define HAVE_ATANF 1
222 float atanf (float x);
224 float
225 atanf (float x)
227 return (float) atan (x);
229 #endif
231 #if HAVE_ATANH && !HAVE_ATANHF
232 float atanhf (float x);
234 float
235 atanhf (float x)
237 return (float) atanh ((double) x);
239 #endif
241 #ifndef HAVE_CEILF
242 #define HAVE_CEILF 1
243 float ceilf (float x);
245 float
246 ceilf (float x)
248 return (float) ceil (x);
250 #endif
252 #ifndef HAVE_COPYSIGNF
253 #define HAVE_COPYSIGNF 1
254 float copysignf (float x, float y);
256 float
257 copysignf (float x, float y)
259 return (float) copysign (x, y);
261 #endif
263 #ifndef HAVE_COSF
264 #define HAVE_COSF 1
265 float cosf (float x);
267 float
268 cosf (float x)
270 return (float) cos (x);
272 #endif
274 #ifndef HAVE_COSHF
275 #define HAVE_COSHF 1
276 float coshf (float x);
278 float
279 coshf (float x)
281 return (float) cosh (x);
283 #endif
285 #ifndef HAVE_EXPF
286 #define HAVE_EXPF 1
287 float expf (float x);
289 float
290 expf (float x)
292 return (float) exp (x);
294 #endif
296 #ifndef HAVE_FABSF
297 #define HAVE_FABSF 1
298 float fabsf (float x);
300 float
301 fabsf (float x)
303 return (float) fabs (x);
305 #endif
307 #ifndef HAVE_FLOORF
308 #define HAVE_FLOORF 1
309 float floorf (float x);
311 float
312 floorf (float x)
314 return (float) floor (x);
316 #endif
318 #ifndef HAVE_FMODF
319 #define HAVE_FMODF 1
320 float fmodf (float x, float y);
322 float
323 fmodf (float x, float y)
325 return (float) fmod (x, y);
327 #endif
329 #ifndef HAVE_FREXPF
330 #define HAVE_FREXPF 1
331 float frexpf (float x, int *exp);
333 float
334 frexpf (float x, int *exp)
336 return (float) frexp (x, exp);
338 #endif
340 #ifndef HAVE_HYPOTF
341 #define HAVE_HYPOTF 1
342 float hypotf (float x, float y);
344 float
345 hypotf (float x, float y)
347 return (float) hypot (x, y);
349 #endif
351 #ifndef HAVE_LOGF
352 #define HAVE_LOGF 1
353 float logf (float x);
355 float
356 logf (float x)
358 return (float) log (x);
360 #endif
362 #ifndef HAVE_LOG10F
363 #define HAVE_LOG10F 1
364 float log10f (float x);
366 float
367 log10f (float x)
369 return (float) log10 (x);
371 #endif
373 #ifndef HAVE_SCALBN
374 #define HAVE_SCALBN 1
375 double scalbn (double x, int y);
377 double
378 scalbn (double x, int y)
380 #if (FLT_RADIX == 2) && defined(HAVE_LDEXP)
381 return ldexp (x, y);
382 #else
383 return x * pow (FLT_RADIX, y);
384 #endif
386 #endif
388 #ifndef HAVE_SCALBNF
389 #define HAVE_SCALBNF 1
390 float scalbnf (float x, int y);
392 float
393 scalbnf (float x, int y)
395 return (float) scalbn (x, y);
397 #endif
399 #ifndef HAVE_SINF
400 #define HAVE_SINF 1
401 float sinf (float x);
403 float
404 sinf (float x)
406 return (float) sin (x);
408 #endif
410 #ifndef HAVE_SINHF
411 #define HAVE_SINHF 1
412 float sinhf (float x);
414 float
415 sinhf (float x)
417 return (float) sinh (x);
419 #endif
421 #ifndef HAVE_SQRTF
422 #define HAVE_SQRTF 1
423 float sqrtf (float x);
425 float
426 sqrtf (float x)
428 return (float) sqrt (x);
430 #endif
432 #ifndef HAVE_TANF
433 #define HAVE_TANF 1
434 float tanf (float x);
436 float
437 tanf (float x)
439 return (float) tan (x);
441 #endif
443 #ifndef HAVE_TANHF
444 #define HAVE_TANHF 1
445 float tanhf (float x);
447 float
448 tanhf (float x)
450 return (float) tanh (x);
452 #endif
454 #ifndef HAVE_TRUNC
455 #define HAVE_TRUNC 1
456 double trunc (double x);
458 double
459 trunc (double x)
461 if (!isfinite (x))
462 return x;
464 if (x < 0.0)
465 return - floor (-x);
466 else
467 return floor (x);
469 #endif
471 #ifndef HAVE_TRUNCF
472 #define HAVE_TRUNCF 1
473 float truncf (float x);
475 float
476 truncf (float x)
478 return (float) trunc (x);
480 #endif
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);
489 float
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))
498 return x + y;
499 if (x == y)
500 return x;
501 if (!isfinite (x))
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__;
510 else
511 denorm_min = __FLT_DENORM_MIN__;
513 if (absx < __FLT_MIN__)
514 delta = denorm_min;
515 else
517 float frac;
518 int exp;
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))
533 delta *= 0.5f;
535 /* If that underflows to zero, then we're back to the minimum. */
536 if (delta == 0.0f)
537 delta = denorm_min;
540 if (y < x)
541 delta = -delta;
543 return x + delta;
545 #endif
548 #if !defined(HAVE_POWF) || defined(HAVE_BROKEN_POWF)
549 #ifndef HAVE_POWF
550 #define HAVE_POWF 1
551 #endif
552 float powf (float x, float y);
554 float
555 powf (float x, float y)
557 return (float) pow (x, y);
559 #endif
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. */
572 long double
573 roundl (long double x)
575 long double t;
576 if (!isfinite (x))
577 return (x);
579 if (x >= 0.0)
581 t = ceill (x);
582 if (t - x > 0.5)
583 t -= 1.0;
584 return (t);
586 else
588 t = ceill (-x);
589 if (t + x > 0.5)
590 t -= 1.0;
591 return (-t);
594 #else
596 /* Poor version of roundl for system that don't have ceill. */
597 long double
598 roundl (long double x)
600 if (x > DBL_MAX || x < -DBL_MAX)
602 #ifdef HAVE_NEXTAFTERL
603 long double prechalf = nextafterl (0.5L, LDBL_MAX);
604 #else
605 static long double prechalf = 0.5L;
606 #endif
607 return (GFC_INTEGER_LARGEST) (x + (x > 0 ? prechalf : -prechalf));
609 else
610 /* Use round(). */
611 return round ((double) x);
614 #endif
615 #endif
617 #ifndef HAVE_ROUND
618 #define HAVE_ROUND 1
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);
623 double
624 round (double x)
626 double t;
627 if (!isfinite (x))
628 return (x);
630 if (x >= 0.0)
632 t = floor (x);
633 if (t - x <= -0.5)
634 t += 1.0;
635 return (t);
637 else
639 t = floor (-x);
640 if (t + x <= -0.5)
641 t += 1.0;
642 return (-t);
645 #endif
647 #ifndef HAVE_ROUNDF
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);
653 float
654 roundf (float x)
656 float t;
657 if (!isfinite (x))
658 return (x);
660 if (x >= 0.0)
662 t = floorf (x);
663 if (t - x <= -0.5)
664 t += 1.0;
665 return (t);
667 else
669 t = floorf (-x);
670 if (t + x <= -0.5)
671 t += 1.0;
672 return (-t);
675 #endif
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);
684 long int
685 lroundf (float x)
687 return (long int) roundf (x);
689 #endif
691 #if !defined(HAVE_LROUND) && defined(HAVE_ROUND)
692 #define HAVE_LROUND 1
693 long int lround (double x);
695 long int
696 lround (double x)
698 return (long int) round (x);
700 #endif
702 #if !defined(HAVE_LROUNDL) && defined(HAVE_ROUNDL)
703 #define HAVE_LROUNDL 1
704 long int lroundl (long double x);
706 long int
707 lroundl (long double x)
709 return (long long int) roundl (x);
711 #endif
713 #if !defined(HAVE_LLROUNDF) && defined(HAVE_ROUNDF)
714 #define HAVE_LLROUNDF 1
715 long long int llroundf (float x);
717 long long int
718 llroundf (float x)
720 return (long long int) roundf (x);
722 #endif
724 #if !defined(HAVE_LLROUND) && defined(HAVE_ROUND)
725 #define HAVE_LLROUND 1
726 long long int llround (double x);
728 long long int
729 llround (double x)
731 return (long long int) round (x);
733 #endif
735 #if !defined(HAVE_LLROUNDL) && defined(HAVE_ROUNDL)
736 #define HAVE_LLROUNDL 1
737 long long int llroundl (long double x);
739 long long int
740 llroundl (long double x)
742 return (long long int) roundl (x);
744 #endif
747 #ifndef HAVE_LOG10L
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);
753 long double
754 log10l (long double x)
756 #if LDBL_MAX_EXP > DBL_MAX_EXP
757 if (x > DBL_MAX)
759 double val;
760 int p2_result = 0;
761 if (x > 0x1p16383L) { p2_result += 16383; x /= 0x1p16383L; }
762 if (x > 0x1p8191L) { p2_result += 8191; x /= 0x1p8191L; }
763 if (x > 0x1p4095L) { p2_result += 4095; x /= 0x1p4095L; }
764 if (x > 0x1p2047L) { p2_result += 2047; x /= 0x1p2047L; }
765 if (x > 0x1p1023L) { p2_result += 1023; x /= 0x1p1023L; }
766 val = log10 ((double) x);
767 return (val + p2_result * .30102999566398119521373889472449302L);
769 #endif
770 #if LDBL_MIN_EXP < DBL_MIN_EXP
771 if (x < DBL_MIN)
773 double val;
774 int p2_result = 0;
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);
783 #endif
784 return log10 (x);
786 #endif
789 #ifndef HAVE_FLOORL
790 #define HAVE_FLOORL 1
791 long double floorl (long double x);
793 long double
794 floorl (long double x)
796 /* Zero, possibly signed. */
797 if (x == 0)
798 return x;
800 /* Large magnitude. */
801 if (x > DBL_MAX || x < (-DBL_MAX))
802 return x;
804 /* Small positive values. */
805 if (x >= 0 && x < DBL_MIN)
806 return 0;
808 /* Small negative values. */
809 if (x < 0 && x > (-DBL_MIN))
810 return -1;
812 return floor (x);
814 #endif
817 #ifndef HAVE_FMODL
818 #define HAVE_FMODL 1
819 long double fmodl (long double x, long double y);
821 long double
822 fmodl (long double x, long double y)
824 if (y == 0.0L)
825 return 0.0L;
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;
831 #endif
834 #if !defined(HAVE_CABSF)
835 #define HAVE_CABSF 1
836 float cabsf (float complex z);
838 float
839 cabsf (float complex z)
841 return hypotf (REALPART (z), IMAGPART (z));
843 #endif
845 #if !defined(HAVE_CABS)
846 #define HAVE_CABS 1
847 double cabs (double complex z);
849 double
850 cabs (double complex z)
852 return hypot (REALPART (z), IMAGPART (z));
854 #endif
856 #if !defined(HAVE_CABSL) && defined(HAVE_HYPOTL)
857 #define HAVE_CABSL 1
858 long double cabsl (long double complex z);
860 long double
861 cabsl (long double complex z)
863 return hypotl (REALPART (z), IMAGPART (z));
865 #endif
868 #if !defined(HAVE_CARGF)
869 #define HAVE_CARGF 1
870 float cargf (float complex z);
872 float
873 cargf (float complex z)
875 return atan2f (IMAGPART (z), REALPART (z));
877 #endif
879 #if !defined(HAVE_CARG)
880 #define HAVE_CARG 1
881 double carg (double complex z);
883 double
884 carg (double complex z)
886 return atan2 (IMAGPART (z), REALPART (z));
888 #endif
890 #if !defined(HAVE_CARGL) && defined(HAVE_ATAN2L)
891 #define HAVE_CARGL 1
892 long double cargl (long double complex z);
894 long double
895 cargl (long double complex z)
897 return atan2l (IMAGPART (z), REALPART (z));
899 #endif
902 /* exp(z) = exp(a)*(cos(b) + i sin(b)) */
903 #if !defined(HAVE_CEXPF)
904 #define HAVE_CEXPF 1
905 float complex cexpf (float complex z);
907 float complex
908 cexpf (float complex z)
910 float a, b;
911 float complex v;
913 a = REALPART (z);
914 b = IMAGPART (z);
915 COMPLEX_ASSIGN (v, cosf (b), sinf (b));
916 return expf (a) * v;
918 #endif
920 #if !defined(HAVE_CEXP)
921 #define HAVE_CEXP 1
922 double complex cexp (double complex z);
924 double complex
925 cexp (double complex z)
927 double a, b;
928 double complex v;
930 a = REALPART (z);
931 b = IMAGPART (z);
932 COMPLEX_ASSIGN (v, cos (b), sin (b));
933 return exp (a) * v;
935 #endif
937 #if !defined(HAVE_CEXPL) && defined(HAVE_COSL) && defined(HAVE_SINL) && defined(EXPL)
938 #define HAVE_CEXPL 1
939 long double complex cexpl (long double complex z);
941 long double complex
942 cexpl (long double complex z)
944 long double a, b;
945 long double complex v;
947 a = REALPART (z);
948 b = IMAGPART (z);
949 COMPLEX_ASSIGN (v, cosl (b), sinl (b));
950 return expl (a) * v;
952 #endif
955 /* log(z) = log (cabs(z)) + i*carg(z) */
956 #if !defined(HAVE_CLOGF)
957 #define HAVE_CLOGF 1
958 float complex clogf (float complex z);
960 float complex
961 clogf (float complex z)
963 float complex v;
965 COMPLEX_ASSIGN (v, logf (cabsf (z)), cargf (z));
966 return v;
968 #endif
970 #if !defined(HAVE_CLOG)
971 #define HAVE_CLOG 1
972 double complex clog (double complex z);
974 double complex
975 clog (double complex z)
977 double complex v;
979 COMPLEX_ASSIGN (v, log (cabs (z)), carg (z));
980 return v;
982 #endif
984 #if !defined(HAVE_CLOGL) && defined(HAVE_LOGL) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
985 #define HAVE_CLOGL 1
986 long double complex clogl (long double complex z);
988 long double complex
989 clogl (long double complex z)
991 long double complex v;
993 COMPLEX_ASSIGN (v, logl (cabsl (z)), cargl (z));
994 return v;
996 #endif
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);
1004 float complex
1005 clog10f (float complex z)
1007 float complex v;
1009 COMPLEX_ASSIGN (v, log10f (cabsf (z)), cargf (z));
1010 return v;
1012 #endif
1014 #if !defined(HAVE_CLOG10)
1015 #define HAVE_CLOG10 1
1016 double complex clog10 (double complex z);
1018 double complex
1019 clog10 (double complex z)
1021 double complex v;
1023 COMPLEX_ASSIGN (v, log10 (cabs (z)), carg (z));
1024 return v;
1026 #endif
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);
1032 long double complex
1033 clog10l (long double complex z)
1035 long double complex v;
1037 COMPLEX_ASSIGN (v, log10l (cabsl (z)), cargl (z));
1038 return v;
1040 #endif
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);
1048 float complex
1049 cpowf (float complex base, float complex power)
1051 return cexpf (power * clogf (base));
1053 #endif
1055 #if !defined(HAVE_CPOW)
1056 #define HAVE_CPOW 1
1057 double complex cpow (double complex base, double complex power);
1059 double complex
1060 cpow (double complex base, double complex power)
1062 return cexp (power * clog (base));
1064 #endif
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);
1070 long double complex
1071 cpowl (long double complex base, long double complex power)
1073 return cexpl (power * clogl (base));
1075 #endif
1078 /* sqrt(z). Algorithm pulled from glibc. */
1079 #if !defined(HAVE_CSQRTF)
1080 #define HAVE_CSQRTF 1
1081 float complex csqrtf (float complex z);
1083 float complex
1084 csqrtf (float complex z)
1086 float re, im;
1087 float complex v;
1089 re = REALPART (z);
1090 im = IMAGPART (z);
1091 if (im == 0)
1093 if (re < 0)
1095 COMPLEX_ASSIGN (v, 0, copysignf (sqrtf (-re), im));
1097 else
1099 COMPLEX_ASSIGN (v, fabsf (sqrtf (re)), copysignf (0, im));
1102 else if (re == 0)
1104 float r;
1106 r = sqrtf (0.5 * fabsf (im));
1108 COMPLEX_ASSIGN (v, r, copysignf (r, im));
1110 else
1112 float d, r, s;
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. */
1117 if (re > 0)
1119 r = sqrtf (0.5 * d + 0.5 * re);
1120 s = (0.5 * im) / r;
1122 else
1124 s = sqrtf (0.5 * d - 0.5 * re);
1125 r = fabsf ((0.5 * im) / s);
1128 COMPLEX_ASSIGN (v, r, copysignf (s, im));
1130 return v;
1132 #endif
1134 #if !defined(HAVE_CSQRT)
1135 #define HAVE_CSQRT 1
1136 double complex csqrt (double complex z);
1138 double complex
1139 csqrt (double complex z)
1141 double re, im;
1142 double complex v;
1144 re = REALPART (z);
1145 im = IMAGPART (z);
1146 if (im == 0)
1148 if (re < 0)
1150 COMPLEX_ASSIGN (v, 0, copysign (sqrt (-re), im));
1152 else
1154 COMPLEX_ASSIGN (v, fabs (sqrt (re)), copysign (0, im));
1157 else if (re == 0)
1159 double r;
1161 r = sqrt (0.5 * fabs (im));
1163 COMPLEX_ASSIGN (v, r, copysign (r, im));
1165 else
1167 double d, r, s;
1169 d = hypot (re, im);
1170 /* Use the identity 2 Re res Im res = Im x
1171 to avoid cancellation error in d +/- Re x. */
1172 if (re > 0)
1174 r = sqrt (0.5 * d + 0.5 * re);
1175 s = (0.5 * im) / r;
1177 else
1179 s = sqrt (0.5 * d - 0.5 * re);
1180 r = fabs ((0.5 * im) / s);
1183 COMPLEX_ASSIGN (v, r, copysign (s, im));
1185 return v;
1187 #endif
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);
1193 long double complex
1194 csqrtl (long double complex z)
1196 long double re, im;
1197 long double complex v;
1199 re = REALPART (z);
1200 im = IMAGPART (z);
1201 if (im == 0)
1203 if (re < 0)
1205 COMPLEX_ASSIGN (v, 0, copysignl (sqrtl (-re), im));
1207 else
1209 COMPLEX_ASSIGN (v, fabsl (sqrtl (re)), copysignl (0, im));
1212 else if (re == 0)
1214 long double r;
1216 r = sqrtl (0.5 * fabsl (im));
1218 COMPLEX_ASSIGN (v, copysignl (r, im), r);
1220 else
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. */
1227 if (re > 0)
1229 r = sqrtl (0.5 * d + 0.5 * re);
1230 s = (0.5 * im) / r;
1232 else
1234 s = sqrtl (0.5 * d - 0.5 * re);
1235 r = fabsl ((0.5 * im) / s);
1238 COMPLEX_ASSIGN (v, r, copysignl (s, im));
1240 return v;
1242 #endif
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);
1250 float complex
1251 csinhf (float complex a)
1253 float r, i;
1254 float complex v;
1256 r = REALPART (a);
1257 i = IMAGPART (a);
1258 COMPLEX_ASSIGN (v, sinhf (r) * cosf (i), coshf (r) * sinf (i));
1259 return v;
1261 #endif
1263 #if !defined(HAVE_CSINH)
1264 #define HAVE_CSINH 1
1265 double complex csinh (double complex a);
1267 double complex
1268 csinh (double complex a)
1270 double r, i;
1271 double complex v;
1273 r = REALPART (a);
1274 i = IMAGPART (a);
1275 COMPLEX_ASSIGN (v, sinh (r) * cos (i), cosh (r) * sin (i));
1276 return v;
1278 #endif
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);
1284 long double complex
1285 csinhl (long double complex a)
1287 long double r, i;
1288 long double complex v;
1290 r = REALPART (a);
1291 i = IMAGPART (a);
1292 COMPLEX_ASSIGN (v, sinhl (r) * cosl (i), coshl (r) * sinl (i));
1293 return v;
1295 #endif
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);
1303 float complex
1304 ccoshf (float complex a)
1306 float r, i;
1307 float complex v;
1309 r = REALPART (a);
1310 i = IMAGPART (a);
1311 COMPLEX_ASSIGN (v, coshf (r) * cosf (i), sinhf (r) * sinf (i));
1312 return v;
1314 #endif
1316 #if !defined(HAVE_CCOSH)
1317 #define HAVE_CCOSH 1
1318 double complex ccosh (double complex a);
1320 double complex
1321 ccosh (double complex a)
1323 double r, i;
1324 double complex v;
1326 r = REALPART (a);
1327 i = IMAGPART (a);
1328 COMPLEX_ASSIGN (v, cosh (r) * cos (i), sinh (r) * sin (i));
1329 return v;
1331 #endif
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);
1337 long double complex
1338 ccoshl (long double complex a)
1340 long double r, i;
1341 long double complex v;
1343 r = REALPART (a);
1344 i = IMAGPART (a);
1345 COMPLEX_ASSIGN (v, coshl (r) * cosl (i), sinhl (r) * sinl (i));
1346 return v;
1348 #endif
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);
1356 float complex
1357 ctanhf (float complex a)
1359 float rt, it;
1360 float complex n, d;
1362 rt = tanhf (REALPART (a));
1363 it = tanf (IMAGPART (a));
1364 COMPLEX_ASSIGN (n, rt, it);
1365 COMPLEX_ASSIGN (d, 1, rt * it);
1367 return n / d;
1369 #endif
1371 #if !defined(HAVE_CTANH)
1372 #define HAVE_CTANH 1
1373 double complex ctanh (double complex a);
1374 double complex
1375 ctanh (double complex a)
1377 double rt, it;
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);
1385 return n / d;
1387 #endif
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);
1393 long double complex
1394 ctanhl (long double complex a)
1396 long double rt, it;
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);
1404 return n / d;
1406 #endif
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);
1414 float complex
1415 csinf (float complex a)
1417 float r, i;
1418 float complex v;
1420 r = REALPART (a);
1421 i = IMAGPART (a);
1422 COMPLEX_ASSIGN (v, sinf (r) * coshf (i), cosf (r) * sinhf (i));
1423 return v;
1425 #endif
1427 #if !defined(HAVE_CSIN)
1428 #define HAVE_CSIN 1
1429 double complex csin (double complex a);
1431 double complex
1432 csin (double complex a)
1434 double r, i;
1435 double complex v;
1437 r = REALPART (a);
1438 i = IMAGPART (a);
1439 COMPLEX_ASSIGN (v, sin (r) * cosh (i), cos (r) * sinh (i));
1440 return v;
1442 #endif
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);
1448 long double complex
1449 csinl (long double complex a)
1451 long double r, i;
1452 long double complex v;
1454 r = REALPART (a);
1455 i = IMAGPART (a);
1456 COMPLEX_ASSIGN (v, sinl (r) * coshl (i), cosl (r) * sinhl (i));
1457 return v;
1459 #endif
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);
1467 float complex
1468 ccosf (float complex a)
1470 float r, i;
1471 float complex v;
1473 r = REALPART (a);
1474 i = IMAGPART (a);
1475 COMPLEX_ASSIGN (v, cosf (r) * coshf (i), - (sinf (r) * sinhf (i)));
1476 return v;
1478 #endif
1480 #if !defined(HAVE_CCOS)
1481 #define HAVE_CCOS 1
1482 double complex ccos (double complex a);
1484 double complex
1485 ccos (double complex a)
1487 double r, i;
1488 double complex v;
1490 r = REALPART (a);
1491 i = IMAGPART (a);
1492 COMPLEX_ASSIGN (v, cos (r) * cosh (i), - (sin (r) * sinh (i)));
1493 return v;
1495 #endif
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);
1501 long double complex
1502 ccosl (long double complex a)
1504 long double r, i;
1505 long double complex v;
1507 r = REALPART (a);
1508 i = IMAGPART (a);
1509 COMPLEX_ASSIGN (v, cosl (r) * coshl (i), - (sinl (r) * sinhl (i)));
1510 return v;
1512 #endif
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);
1520 float complex
1521 ctanf (float complex a)
1523 float rt, it;
1524 float complex n, d;
1526 rt = tanf (REALPART (a));
1527 it = tanhf (IMAGPART (a));
1528 COMPLEX_ASSIGN (n, rt, it);
1529 COMPLEX_ASSIGN (d, 1, - (rt * it));
1531 return n / d;
1533 #endif
1535 #if !defined(HAVE_CTAN)
1536 #define HAVE_CTAN 1
1537 double complex ctan (double complex a);
1539 double complex
1540 ctan (double complex a)
1542 double rt, it;
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));
1550 return n / d;
1552 #endif
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);
1558 long double complex
1559 ctanl (long double complex a)
1561 long double rt, it;
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));
1569 return n / d;
1571 #endif
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);
1581 complex float
1582 casinf (complex float z)
1584 return -I*clogf (I*z + csqrtf (1.0f-z*z));
1586 #endif
1589 #if !defined(HAVE_CASIN) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1590 #define HAVE_CASIN 1
1591 complex double casin (complex double z);
1593 complex double
1594 casin (complex double z)
1596 return -I*clog (I*z + csqrt (1.0-z*z));
1598 #endif
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);
1605 complex long double
1606 casinl (complex long double z)
1608 return -I*clogl (I*z + csqrtl (1.0L-z*z));
1610 #endif
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);
1620 complex float
1621 cacosf (complex float z)
1623 return -I*clogf (z + I*csqrtf (1.0f-z*z));
1625 #endif
1628 #if !defined(HAVE_CACOS) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1629 #define HAVE_CACOS 1
1630 complex double cacos (complex double z);
1632 complex double
1633 cacos (complex double z)
1635 return -I*clog (z + I*csqrt (1.0-z*z));
1637 #endif
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);
1644 complex long double
1645 cacosl (complex long double z)
1647 return -I*clogl (z + I*csqrtl (1.0L-z*z));
1649 #endif
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);
1659 complex float
1660 catanf (complex float z)
1662 return I*clogf ((I+z)/(I-z))/2.0f;
1664 #endif
1667 #if !defined(HAVE_CATAN) && defined(HAVE_CLOG)
1668 #define HAVE_CACOS 1
1669 complex double catan (complex double z);
1671 complex double
1672 catan (complex double z)
1674 return I*clog ((I+z)/(I-z))/2.0;
1676 #endif
1679 #if !defined(HAVE_CATANL) && defined(HAVE_CLOGL)
1680 #define HAVE_CACOSL 1
1681 complex long double catanl (complex long double z);
1683 complex long double
1684 catanl (complex long double z)
1686 return I*clogl ((I+z)/(I-z))/2.0L;
1688 #endif
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);
1698 complex float
1699 casinhf (complex float z)
1701 return clogf (z + csqrtf (z*z+1.0f));
1703 #endif
1706 #if !defined(HAVE_CASINH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1707 #define HAVE_CASINH 1
1708 complex double casinh (complex double z);
1710 complex double
1711 casinh (complex double z)
1713 return clog (z + csqrt (z*z+1.0));
1715 #endif
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);
1722 complex long double
1723 casinhl (complex long double z)
1725 return clogl (z + csqrtl (z*z+1.0L));
1727 #endif
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);
1737 complex float
1738 cacoshf (complex float z)
1740 return clogf (z + csqrtf (z-1.0f) * csqrtf (z+1.0f));
1742 #endif
1745 #if !defined(HAVE_CACOSH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1746 #define HAVE_CACOSH 1
1747 complex double cacosh (complex double z);
1749 complex double
1750 cacosh (complex double z)
1752 return clog (z + csqrt (z-1.0) * csqrt (z+1.0));
1754 #endif
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);
1761 complex long double
1762 cacoshl (complex long double z)
1764 return clogl (z + csqrtl (z-1.0L) * csqrtl (z+1.0L));
1766 #endif
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);
1776 complex float
1777 catanhf (complex float z)
1779 return clogf ((1.0f+z)/(1.0f-z))/2.0f;
1781 #endif
1784 #if !defined(HAVE_CATANH) && defined(HAVE_CLOG)
1785 #define HAVE_CATANH 1
1786 complex double catanh (complex double z);
1788 complex double
1789 catanh (complex double z)
1791 return clog ((1.0+z)/(1.0-z))/2.0;
1793 #endif
1795 #if !defined(HAVE_CATANHL) && defined(HAVE_CLOGL)
1796 #define HAVE_CATANHL 1
1797 complex long double catanhl (complex long double z);
1799 complex long double
1800 catanhl (complex long double z)
1802 return clogl ((1.0L+z)/(1.0L-z))/2.0L;
1804 #endif
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. */
1814 #undef SQRTPI
1815 #define SQRTPI 0.9189385332046727417803297
1817 #undef PI
1818 #define PI 3.1415926535897932384626434
1820 double
1821 tgamma (double x)
1823 int i, n, parity;
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;
1848 if (eps == 0)
1849 eps = nextafter (1., 2.) - 1.;
1851 parity = 0;
1852 fact = 1;
1853 n = 0;
1854 y = x;
1856 if (isnan (x))
1857 return x;
1859 if (y <= 0)
1861 y = -x;
1862 y1 = trunc (y);
1863 res = y - y1;
1865 if (res != 0)
1867 if (y1 != trunc (y1*0.5l)*2)
1868 parity = 1;
1869 fact = -PI / sin (PI*res);
1870 y = y + 1;
1872 else
1873 return x == 0 ? copysign (xinf, x) : xnan;
1876 if (y < eps)
1878 if (y >= xminin)
1879 res = 1 / y;
1880 else
1881 return xinf;
1883 else if (y < 13)
1885 y1 = y;
1886 if (y < 1)
1888 z = y;
1889 y = y + 1;
1891 else
1893 n = (int)y - 1;
1894 y = y - n;
1895 z = y - 1;
1898 xnum = 0;
1899 xden = 1;
1900 for (i = 0; i < 8; i++)
1902 xnum = (xnum + p[i]) * z;
1903 xden = xden * z + q[i];
1906 res = xnum / xden + 1;
1908 if (y1 < y)
1909 res = res / y1;
1910 else if (y1 > y)
1911 for (i = 1; i <= n; i++)
1913 res = res * y;
1914 y = y + 1;
1917 else
1919 if (y < xbig)
1921 ysq = y * y;
1922 sum = c[6];
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);
1928 res = exp (sum);
1930 else
1931 return x < 0 ? xnan : xinf;
1934 if (parity)
1935 res = -res;
1936 if (fact != 1)
1937 res = fact / res;
1939 return res;
1941 #endif
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)
1956 double
1957 lgamma (double y)
1960 #undef SQRTPI
1961 #define SQRTPI 0.9189385332046727417803297
1963 #undef PI
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,
2005 5.7083835261e-03 };
2007 static double xbig = 2.55e305, xinf = __builtin_inf (), eps = 0,
2008 frtbig = 2.25e76;
2010 int i;
2011 double corr, res, xden, xm1, xm2, xm4, xnum, ysq;
2013 if (eps == 0)
2014 eps = __builtin_nextafter (1., 2.) - 1.;
2016 if ((y > 0) && (y <= xbig))
2018 if (y <= eps)
2019 res = -log (y);
2020 else if (y <= 1.5)
2022 if (y < PNT68)
2024 corr = -log (y);
2025 xm1 = y;
2027 else
2029 corr = 0;
2030 xm1 = (y - 0.5) - 0.5;
2033 if ((y <= 0.5) || (y >= PNT68))
2035 xden = 1;
2036 xnum = 0;
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)));
2044 else
2046 xm2 = (y - 0.5) - 0.5;
2047 xden = 1;
2048 xnum = 0;
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));
2057 else if (y <= 4)
2059 xm2 = y - 2;
2060 xden = 1;
2061 xnum = 0;
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));
2069 else if (y <= 12)
2071 xm4 = y - 4;
2072 xden = -1;
2073 xnum = 0;
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);
2081 else
2083 res = 0;
2084 if (y <= frtbig)
2086 res = c[6];
2087 ysq = y * y;
2088 for (i = 0; i < 6; i++)
2089 res = res / ysq + c[i];
2091 res = res/y;
2092 corr = log (y);
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. */
2102 if (y > -1.e-100)
2103 res = -2 * log (fabs (y)) - lgamma (-y);
2104 else
2105 res = log (PI / fabs (y * sin (PI * y))) - lgamma (-y);
2107 else
2108 res = xinf;
2110 return res;
2112 #endif
2115 #if defined(HAVE_TGAMMA) && !defined(HAVE_TGAMMAF)
2116 #define HAVE_TGAMMAF 1
2117 float tgammaf (float);
2119 float
2120 tgammaf (float x)
2122 return (float) tgamma ((double) x);
2124 #endif
2126 #if defined(HAVE_LGAMMA) && !defined(HAVE_LGAMMAF)
2127 #define HAVE_LGAMMAF 1
2128 float lgammaf (float);
2130 float
2131 lgammaf (float x)
2133 return (float) lgamma ((double) x);
2135 #endif