Daily bump.
[official-gcc.git] / libgfortran / intrinsics / c99_functions.c
blobee74b115ea7e33a8d87d5bfc1be28ed3712cb926
1 /* Implementation of various C99 functions
2 Copyright (C) 2004-2013 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 /* On a C99 system "I" (with I*I = -1) should be defined in complex.h;
31 if not, we define a fallback version here. */
32 #ifndef I
33 # if defined(_Imaginary_I)
34 # define I _Imaginary_I
35 # elif defined(_Complex_I)
36 # define I _Complex_I
37 # else
38 # define I (1.0fi)
39 # endif
40 #endif
42 /* Prototypes are included to silence -Wstrict-prototypes
43 -Wmissing-prototypes. */
46 /* Wrappers for systems without the various C99 single precision Bessel
47 functions. */
49 #if defined(HAVE_J0) && ! defined(HAVE_J0F)
50 #define HAVE_J0F 1
51 float j0f (float);
53 float
54 j0f (float x)
56 return (float) j0 ((double) x);
58 #endif
60 #if defined(HAVE_J1) && !defined(HAVE_J1F)
61 #define HAVE_J1F 1
62 float j1f (float);
64 float j1f (float x)
66 return (float) j1 ((double) x);
68 #endif
70 #if defined(HAVE_JN) && !defined(HAVE_JNF)
71 #define HAVE_JNF 1
72 float jnf (int, float);
74 float
75 jnf (int n, float x)
77 return (float) jn (n, (double) x);
79 #endif
81 #if defined(HAVE_Y0) && !defined(HAVE_Y0F)
82 #define HAVE_Y0F 1
83 float y0f (float);
85 float
86 y0f (float x)
88 return (float) y0 ((double) x);
90 #endif
92 #if defined(HAVE_Y1) && !defined(HAVE_Y1F)
93 #define HAVE_Y1F 1
94 float y1f (float);
96 float
97 y1f (float x)
99 return (float) y1 ((double) x);
101 #endif
103 #if defined(HAVE_YN) && !defined(HAVE_YNF)
104 #define HAVE_YNF 1
105 float ynf (int, float);
107 float
108 ynf (int n, float x)
110 return (float) yn (n, (double) x);
112 #endif
115 /* Wrappers for systems without the C99 erff() and erfcf() functions. */
117 #if defined(HAVE_ERF) && !defined(HAVE_ERFF)
118 #define HAVE_ERFF 1
119 float erff (float);
121 float
122 erff (float x)
124 return (float) erf ((double) x);
126 #endif
128 #if defined(HAVE_ERFC) && !defined(HAVE_ERFCF)
129 #define HAVE_ERFCF 1
130 float erfcf (float);
132 float
133 erfcf (float x)
135 return (float) erfc ((double) x);
137 #endif
140 #ifndef HAVE_ACOSF
141 #define HAVE_ACOSF 1
142 float acosf (float x);
144 float
145 acosf (float x)
147 return (float) acos (x);
149 #endif
151 #if HAVE_ACOSH && !HAVE_ACOSHF
152 float acoshf (float x);
154 float
155 acoshf (float x)
157 return (float) acosh ((double) x);
159 #endif
161 #ifndef HAVE_ASINF
162 #define HAVE_ASINF 1
163 float asinf (float x);
165 float
166 asinf (float x)
168 return (float) asin (x);
170 #endif
172 #if HAVE_ASINH && !HAVE_ASINHF
173 float asinhf (float x);
175 float
176 asinhf (float x)
178 return (float) asinh ((double) x);
180 #endif
182 #ifndef HAVE_ATAN2F
183 #define HAVE_ATAN2F 1
184 float atan2f (float y, float x);
186 float
187 atan2f (float y, float x)
189 return (float) atan2 (y, x);
191 #endif
193 #ifndef HAVE_ATANF
194 #define HAVE_ATANF 1
195 float atanf (float x);
197 float
198 atanf (float x)
200 return (float) atan (x);
202 #endif
204 #if HAVE_ATANH && !HAVE_ATANHF
205 float atanhf (float x);
207 float
208 atanhf (float x)
210 return (float) atanh ((double) x);
212 #endif
214 #ifndef HAVE_CEILF
215 #define HAVE_CEILF 1
216 float ceilf (float x);
218 float
219 ceilf (float x)
221 return (float) ceil (x);
223 #endif
225 #ifndef HAVE_COPYSIGNF
226 #define HAVE_COPYSIGNF 1
227 float copysignf (float x, float y);
229 float
230 copysignf (float x, float y)
232 return (float) copysign (x, y);
234 #endif
236 #ifndef HAVE_COSF
237 #define HAVE_COSF 1
238 float cosf (float x);
240 float
241 cosf (float x)
243 return (float) cos (x);
245 #endif
247 #ifndef HAVE_COSHF
248 #define HAVE_COSHF 1
249 float coshf (float x);
251 float
252 coshf (float x)
254 return (float) cosh (x);
256 #endif
258 #ifndef HAVE_EXPF
259 #define HAVE_EXPF 1
260 float expf (float x);
262 float
263 expf (float x)
265 return (float) exp (x);
267 #endif
269 #ifndef HAVE_FABSF
270 #define HAVE_FABSF 1
271 float fabsf (float x);
273 float
274 fabsf (float x)
276 return (float) fabs (x);
278 #endif
280 #ifndef HAVE_FLOORF
281 #define HAVE_FLOORF 1
282 float floorf (float x);
284 float
285 floorf (float x)
287 return (float) floor (x);
289 #endif
291 #ifndef HAVE_FMODF
292 #define HAVE_FMODF 1
293 float fmodf (float x, float y);
295 float
296 fmodf (float x, float y)
298 return (float) fmod (x, y);
300 #endif
302 #ifndef HAVE_FREXPF
303 #define HAVE_FREXPF 1
304 float frexpf (float x, int *exp);
306 float
307 frexpf (float x, int *exp)
309 return (float) frexp (x, exp);
311 #endif
313 #ifndef HAVE_HYPOTF
314 #define HAVE_HYPOTF 1
315 float hypotf (float x, float y);
317 float
318 hypotf (float x, float y)
320 return (float) hypot (x, y);
322 #endif
324 #ifndef HAVE_LOGF
325 #define HAVE_LOGF 1
326 float logf (float x);
328 float
329 logf (float x)
331 return (float) log (x);
333 #endif
335 #ifndef HAVE_LOG10F
336 #define HAVE_LOG10F 1
337 float log10f (float x);
339 float
340 log10f (float x)
342 return (float) log10 (x);
344 #endif
346 #ifndef HAVE_SCALBN
347 #define HAVE_SCALBN 1
348 double scalbn (double x, int y);
350 double
351 scalbn (double x, int y)
353 #if (FLT_RADIX == 2) && defined(HAVE_LDEXP)
354 return ldexp (x, y);
355 #else
356 return x * pow (FLT_RADIX, y);
357 #endif
359 #endif
361 #ifndef HAVE_SCALBNF
362 #define HAVE_SCALBNF 1
363 float scalbnf (float x, int y);
365 float
366 scalbnf (float x, int y)
368 return (float) scalbn (x, y);
370 #endif
372 #ifndef HAVE_SINF
373 #define HAVE_SINF 1
374 float sinf (float x);
376 float
377 sinf (float x)
379 return (float) sin (x);
381 #endif
383 #ifndef HAVE_SINHF
384 #define HAVE_SINHF 1
385 float sinhf (float x);
387 float
388 sinhf (float x)
390 return (float) sinh (x);
392 #endif
394 #ifndef HAVE_SQRTF
395 #define HAVE_SQRTF 1
396 float sqrtf (float x);
398 float
399 sqrtf (float x)
401 return (float) sqrt (x);
403 #endif
405 #ifndef HAVE_TANF
406 #define HAVE_TANF 1
407 float tanf (float x);
409 float
410 tanf (float x)
412 return (float) tan (x);
414 #endif
416 #ifndef HAVE_TANHF
417 #define HAVE_TANHF 1
418 float tanhf (float x);
420 float
421 tanhf (float x)
423 return (float) tanh (x);
425 #endif
427 #ifndef HAVE_TRUNC
428 #define HAVE_TRUNC 1
429 double trunc (double x);
431 double
432 trunc (double x)
434 if (!isfinite (x))
435 return x;
437 if (x < 0.0)
438 return - floor (-x);
439 else
440 return floor (x);
442 #endif
444 #ifndef HAVE_TRUNCF
445 #define HAVE_TRUNCF 1
446 float truncf (float x);
448 float
449 truncf (float x)
451 return (float) trunc (x);
453 #endif
455 #ifndef HAVE_NEXTAFTERF
456 #define HAVE_NEXTAFTERF 1
457 /* This is a portable implementation of nextafterf that is intended to be
458 independent of the floating point format or its in memory representation.
459 This implementation works correctly with denormalized values. */
460 float nextafterf (float x, float y);
462 float
463 nextafterf (float x, float y)
465 /* This variable is marked volatile to avoid excess precision problems
466 on some platforms, including IA-32. */
467 volatile float delta;
468 float absx, denorm_min;
470 if (isnan (x) || isnan (y))
471 return x + y;
472 if (x == y)
473 return x;
474 if (!isfinite (x))
475 return x > 0 ? __FLT_MAX__ : - __FLT_MAX__;
477 /* absx = fabsf (x); */
478 absx = (x < 0.0) ? -x : x;
480 /* __FLT_DENORM_MIN__ is non-zero iff the target supports denormals. */
481 if (__FLT_DENORM_MIN__ == 0.0f)
482 denorm_min = __FLT_MIN__;
483 else
484 denorm_min = __FLT_DENORM_MIN__;
486 if (absx < __FLT_MIN__)
487 delta = denorm_min;
488 else
490 float frac;
491 int exp;
493 /* Discard the fraction from x. */
494 frac = frexpf (absx, &exp);
495 delta = scalbnf (0.5f, exp);
497 /* Scale x by the epsilon of the representation. By rights we should
498 have been able to combine this with scalbnf, but some targets don't
499 get that correct with denormals. */
500 delta *= __FLT_EPSILON__;
502 /* If we're going to be reducing the absolute value of X, and doing so
503 would reduce the exponent of X, then the delta to be applied is
504 one exponent smaller. */
505 if (frac == 0.5f && (y < x) == (x > 0))
506 delta *= 0.5f;
508 /* If that underflows to zero, then we're back to the minimum. */
509 if (delta == 0.0f)
510 delta = denorm_min;
513 if (y < x)
514 delta = -delta;
516 return x + delta;
518 #endif
521 #ifndef HAVE_POWF
522 #define HAVE_POWF 1
523 float powf (float x, float y);
525 float
526 powf (float x, float y)
528 return (float) pow (x, y);
530 #endif
533 #ifndef HAVE_ROUND
534 #define HAVE_ROUND 1
535 /* Round to nearest integral value. If the argument is halfway between two
536 integral values then round away from zero. */
537 double round (double x);
539 double
540 round (double x)
542 double t;
543 if (!isfinite (x))
544 return (x);
546 if (x >= 0.0)
548 t = floor (x);
549 if (t - x <= -0.5)
550 t += 1.0;
551 return (t);
553 else
555 t = floor (-x);
556 if (t + x <= -0.5)
557 t += 1.0;
558 return (-t);
561 #endif
564 /* Algorithm by Steven G. Kargl. */
566 #if !defined(HAVE_ROUNDL)
567 #define HAVE_ROUNDL 1
568 long double roundl (long double x);
570 #if defined(HAVE_CEILL)
571 /* Round to nearest integral value. If the argument is halfway between two
572 integral values then round away from zero. */
574 long double
575 roundl (long double x)
577 long double t;
578 if (!isfinite (x))
579 return (x);
581 if (x >= 0.0)
583 t = ceill (x);
584 if (t - x > 0.5)
585 t -= 1.0;
586 return (t);
588 else
590 t = ceill (-x);
591 if (t + x > 0.5)
592 t -= 1.0;
593 return (-t);
596 #else
598 /* Poor version of roundl for system that don't have ceill. */
599 long double
600 roundl (long double x)
602 if (x > DBL_MAX || x < -DBL_MAX)
604 #ifdef HAVE_NEXTAFTERL
605 long double prechalf = nextafterl (0.5L, LDBL_MAX);
606 #else
607 static long double prechalf = 0.5L;
608 #endif
609 return (GFC_INTEGER_LARGEST) (x + (x > 0 ? prechalf : -prechalf));
611 else
612 /* Use round(). */
613 return round ((double) x);
616 #endif
617 #endif
619 #ifndef HAVE_ROUNDF
620 #define HAVE_ROUNDF 1
621 /* Round to nearest integral value. If the argument is halfway between two
622 integral values then round away from zero. */
623 float roundf (float x);
625 float
626 roundf (float x)
628 float t;
629 if (!isfinite (x))
630 return (x);
632 if (x >= 0.0)
634 t = floorf (x);
635 if (t - x <= -0.5)
636 t += 1.0;
637 return (t);
639 else
641 t = floorf (-x);
642 if (t + x <= -0.5)
643 t += 1.0;
644 return (-t);
647 #endif
650 /* lround{f,,l} and llround{f,,l} functions. */
652 #if !defined(HAVE_LROUNDF) && defined(HAVE_ROUNDF)
653 #define HAVE_LROUNDF 1
654 long int lroundf (float x);
656 long int
657 lroundf (float x)
659 return (long int) roundf (x);
661 #endif
663 #if !defined(HAVE_LROUND) && defined(HAVE_ROUND)
664 #define HAVE_LROUND 1
665 long int lround (double x);
667 long int
668 lround (double x)
670 return (long int) round (x);
672 #endif
674 #if !defined(HAVE_LROUNDL) && defined(HAVE_ROUNDL)
675 #define HAVE_LROUNDL 1
676 long int lroundl (long double x);
678 long int
679 lroundl (long double x)
681 return (long long int) roundl (x);
683 #endif
685 #if !defined(HAVE_LLROUNDF) && defined(HAVE_ROUNDF)
686 #define HAVE_LLROUNDF 1
687 long long int llroundf (float x);
689 long long int
690 llroundf (float x)
692 return (long long int) roundf (x);
694 #endif
696 #if !defined(HAVE_LLROUND) && defined(HAVE_ROUND)
697 #define HAVE_LLROUND 1
698 long long int llround (double x);
700 long long int
701 llround (double x)
703 return (long long int) round (x);
705 #endif
707 #if !defined(HAVE_LLROUNDL) && defined(HAVE_ROUNDL)
708 #define HAVE_LLROUNDL 1
709 long long int llroundl (long double x);
711 long long int
712 llroundl (long double x)
714 return (long long int) roundl (x);
716 #endif
719 #ifndef HAVE_LOG10L
720 #define HAVE_LOG10L 1
721 /* log10 function for long double variables. The version provided here
722 reduces the argument until it fits into a double, then use log10. */
723 long double log10l (long double x);
725 long double
726 log10l (long double x)
728 #if LDBL_MAX_EXP > DBL_MAX_EXP
729 if (x > DBL_MAX)
731 double val;
732 int p2_result = 0;
733 if (x > 0x1p16383L) { p2_result += 16383; x /= 0x1p16383L; }
734 if (x > 0x1p8191L) { p2_result += 8191; x /= 0x1p8191L; }
735 if (x > 0x1p4095L) { p2_result += 4095; x /= 0x1p4095L; }
736 if (x > 0x1p2047L) { p2_result += 2047; x /= 0x1p2047L; }
737 if (x > 0x1p1023L) { p2_result += 1023; x /= 0x1p1023L; }
738 val = log10 ((double) x);
739 return (val + p2_result * .30102999566398119521373889472449302L);
741 #endif
742 #if LDBL_MIN_EXP < DBL_MIN_EXP
743 if (x < DBL_MIN)
745 double val;
746 int p2_result = 0;
747 if (x < 0x1p-16380L) { p2_result += 16380; x /= 0x1p-16380L; }
748 if (x < 0x1p-8189L) { p2_result += 8189; x /= 0x1p-8189L; }
749 if (x < 0x1p-4093L) { p2_result += 4093; x /= 0x1p-4093L; }
750 if (x < 0x1p-2045L) { p2_result += 2045; x /= 0x1p-2045L; }
751 if (x < 0x1p-1021L) { p2_result += 1021; x /= 0x1p-1021L; }
752 val = fabs (log10 ((double) x));
753 return (- val - p2_result * .30102999566398119521373889472449302L);
755 #endif
756 return log10 (x);
758 #endif
761 #ifndef HAVE_FLOORL
762 #define HAVE_FLOORL 1
763 long double floorl (long double x);
765 long double
766 floorl (long double x)
768 /* Zero, possibly signed. */
769 if (x == 0)
770 return x;
772 /* Large magnitude. */
773 if (x > DBL_MAX || x < (-DBL_MAX))
774 return x;
776 /* Small positive values. */
777 if (x >= 0 && x < DBL_MIN)
778 return 0;
780 /* Small negative values. */
781 if (x < 0 && x > (-DBL_MIN))
782 return -1;
784 return floor (x);
786 #endif
789 #ifndef HAVE_FMODL
790 #define HAVE_FMODL 1
791 long double fmodl (long double x, long double y);
793 long double
794 fmodl (long double x, long double y)
796 if (y == 0.0L)
797 return 0.0L;
799 /* Need to check that the result has the same sign as x and magnitude
800 less than the magnitude of y. */
801 return x - floorl (x / y) * y;
803 #endif
806 #if !defined(HAVE_CABSF)
807 #define HAVE_CABSF 1
808 float cabsf (float complex z);
810 float
811 cabsf (float complex z)
813 return hypotf (REALPART (z), IMAGPART (z));
815 #endif
817 #if !defined(HAVE_CABS)
818 #define HAVE_CABS 1
819 double cabs (double complex z);
821 double
822 cabs (double complex z)
824 return hypot (REALPART (z), IMAGPART (z));
826 #endif
828 #if !defined(HAVE_CABSL) && defined(HAVE_HYPOTL)
829 #define HAVE_CABSL 1
830 long double cabsl (long double complex z);
832 long double
833 cabsl (long double complex z)
835 return hypotl (REALPART (z), IMAGPART (z));
837 #endif
840 #if !defined(HAVE_CARGF)
841 #define HAVE_CARGF 1
842 float cargf (float complex z);
844 float
845 cargf (float complex z)
847 return atan2f (IMAGPART (z), REALPART (z));
849 #endif
851 #if !defined(HAVE_CARG)
852 #define HAVE_CARG 1
853 double carg (double complex z);
855 double
856 carg (double complex z)
858 return atan2 (IMAGPART (z), REALPART (z));
860 #endif
862 #if !defined(HAVE_CARGL) && defined(HAVE_ATAN2L)
863 #define HAVE_CARGL 1
864 long double cargl (long double complex z);
866 long double
867 cargl (long double complex z)
869 return atan2l (IMAGPART (z), REALPART (z));
871 #endif
874 /* exp(z) = exp(a)*(cos(b) + i sin(b)) */
875 #if !defined(HAVE_CEXPF)
876 #define HAVE_CEXPF 1
877 float complex cexpf (float complex z);
879 float complex
880 cexpf (float complex z)
882 float a, b;
883 float complex v;
885 a = REALPART (z);
886 b = IMAGPART (z);
887 COMPLEX_ASSIGN (v, cosf (b), sinf (b));
888 return expf (a) * v;
890 #endif
892 #if !defined(HAVE_CEXP)
893 #define HAVE_CEXP 1
894 double complex cexp (double complex z);
896 double complex
897 cexp (double complex z)
899 double a, b;
900 double complex v;
902 a = REALPART (z);
903 b = IMAGPART (z);
904 COMPLEX_ASSIGN (v, cos (b), sin (b));
905 return exp (a) * v;
907 #endif
909 #if !defined(HAVE_CEXPL) && defined(HAVE_COSL) && defined(HAVE_SINL) && defined(EXPL)
910 #define HAVE_CEXPL 1
911 long double complex cexpl (long double complex z);
913 long double complex
914 cexpl (long double complex z)
916 long double a, b;
917 long double complex v;
919 a = REALPART (z);
920 b = IMAGPART (z);
921 COMPLEX_ASSIGN (v, cosl (b), sinl (b));
922 return expl (a) * v;
924 #endif
927 /* log(z) = log (cabs(z)) + i*carg(z) */
928 #if !defined(HAVE_CLOGF)
929 #define HAVE_CLOGF 1
930 float complex clogf (float complex z);
932 float complex
933 clogf (float complex z)
935 float complex v;
937 COMPLEX_ASSIGN (v, logf (cabsf (z)), cargf (z));
938 return v;
940 #endif
942 #if !defined(HAVE_CLOG)
943 #define HAVE_CLOG 1
944 double complex clog (double complex z);
946 double complex
947 clog (double complex z)
949 double complex v;
951 COMPLEX_ASSIGN (v, log (cabs (z)), carg (z));
952 return v;
954 #endif
956 #if !defined(HAVE_CLOGL) && defined(HAVE_LOGL) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
957 #define HAVE_CLOGL 1
958 long double complex clogl (long double complex z);
960 long double complex
961 clogl (long double complex z)
963 long double complex v;
965 COMPLEX_ASSIGN (v, logl (cabsl (z)), cargl (z));
966 return v;
968 #endif
971 /* log10(z) = log10 (cabs(z)) + i*carg(z) */
972 #if !defined(HAVE_CLOG10F)
973 #define HAVE_CLOG10F 1
974 float complex clog10f (float complex z);
976 float complex
977 clog10f (float complex z)
979 float complex v;
981 COMPLEX_ASSIGN (v, log10f (cabsf (z)), cargf (z));
982 return v;
984 #endif
986 #if !defined(HAVE_CLOG10)
987 #define HAVE_CLOG10 1
988 double complex clog10 (double complex z);
990 double complex
991 clog10 (double complex z)
993 double complex v;
995 COMPLEX_ASSIGN (v, log10 (cabs (z)), carg (z));
996 return v;
998 #endif
1000 #if !defined(HAVE_CLOG10L) && defined(HAVE_LOG10L) && defined(HAVE_CABSL) && defined(HAVE_CARGL)
1001 #define HAVE_CLOG10L 1
1002 long double complex clog10l (long double complex z);
1004 long double complex
1005 clog10l (long double complex z)
1007 long double complex v;
1009 COMPLEX_ASSIGN (v, log10l (cabsl (z)), cargl (z));
1010 return v;
1012 #endif
1015 /* pow(base, power) = cexp (power * clog (base)) */
1016 #if !defined(HAVE_CPOWF)
1017 #define HAVE_CPOWF 1
1018 float complex cpowf (float complex base, float complex power);
1020 float complex
1021 cpowf (float complex base, float complex power)
1023 return cexpf (power * clogf (base));
1025 #endif
1027 #if !defined(HAVE_CPOW)
1028 #define HAVE_CPOW 1
1029 double complex cpow (double complex base, double complex power);
1031 double complex
1032 cpow (double complex base, double complex power)
1034 return cexp (power * clog (base));
1036 #endif
1038 #if !defined(HAVE_CPOWL) && defined(HAVE_CEXPL) && defined(HAVE_CLOGL)
1039 #define HAVE_CPOWL 1
1040 long double complex cpowl (long double complex base, long double complex power);
1042 long double complex
1043 cpowl (long double complex base, long double complex power)
1045 return cexpl (power * clogl (base));
1047 #endif
1050 /* sqrt(z). Algorithm pulled from glibc. */
1051 #if !defined(HAVE_CSQRTF)
1052 #define HAVE_CSQRTF 1
1053 float complex csqrtf (float complex z);
1055 float complex
1056 csqrtf (float complex z)
1058 float re, im;
1059 float complex v;
1061 re = REALPART (z);
1062 im = IMAGPART (z);
1063 if (im == 0)
1065 if (re < 0)
1067 COMPLEX_ASSIGN (v, 0, copysignf (sqrtf (-re), im));
1069 else
1071 COMPLEX_ASSIGN (v, fabsf (sqrtf (re)), copysignf (0, im));
1074 else if (re == 0)
1076 float r;
1078 r = sqrtf (0.5 * fabsf (im));
1080 COMPLEX_ASSIGN (v, r, copysignf (r, im));
1082 else
1084 float d, r, s;
1086 d = hypotf (re, im);
1087 /* Use the identity 2 Re res Im res = Im x
1088 to avoid cancellation error in d +/- Re x. */
1089 if (re > 0)
1091 r = sqrtf (0.5 * d + 0.5 * re);
1092 s = (0.5 * im) / r;
1094 else
1096 s = sqrtf (0.5 * d - 0.5 * re);
1097 r = fabsf ((0.5 * im) / s);
1100 COMPLEX_ASSIGN (v, r, copysignf (s, im));
1102 return v;
1104 #endif
1106 #if !defined(HAVE_CSQRT)
1107 #define HAVE_CSQRT 1
1108 double complex csqrt (double complex z);
1110 double complex
1111 csqrt (double complex z)
1113 double re, im;
1114 double complex v;
1116 re = REALPART (z);
1117 im = IMAGPART (z);
1118 if (im == 0)
1120 if (re < 0)
1122 COMPLEX_ASSIGN (v, 0, copysign (sqrt (-re), im));
1124 else
1126 COMPLEX_ASSIGN (v, fabs (sqrt (re)), copysign (0, im));
1129 else if (re == 0)
1131 double r;
1133 r = sqrt (0.5 * fabs (im));
1135 COMPLEX_ASSIGN (v, r, copysign (r, im));
1137 else
1139 double d, r, s;
1141 d = hypot (re, im);
1142 /* Use the identity 2 Re res Im res = Im x
1143 to avoid cancellation error in d +/- Re x. */
1144 if (re > 0)
1146 r = sqrt (0.5 * d + 0.5 * re);
1147 s = (0.5 * im) / r;
1149 else
1151 s = sqrt (0.5 * d - 0.5 * re);
1152 r = fabs ((0.5 * im) / s);
1155 COMPLEX_ASSIGN (v, r, copysign (s, im));
1157 return v;
1159 #endif
1161 #if !defined(HAVE_CSQRTL) && defined(HAVE_COPYSIGNL) && defined(HAVE_SQRTL) && defined(HAVE_FABSL) && defined(HAVE_HYPOTL)
1162 #define HAVE_CSQRTL 1
1163 long double complex csqrtl (long double complex z);
1165 long double complex
1166 csqrtl (long double complex z)
1168 long double re, im;
1169 long double complex v;
1171 re = REALPART (z);
1172 im = IMAGPART (z);
1173 if (im == 0)
1175 if (re < 0)
1177 COMPLEX_ASSIGN (v, 0, copysignl (sqrtl (-re), im));
1179 else
1181 COMPLEX_ASSIGN (v, fabsl (sqrtl (re)), copysignl (0, im));
1184 else if (re == 0)
1186 long double r;
1188 r = sqrtl (0.5 * fabsl (im));
1190 COMPLEX_ASSIGN (v, copysignl (r, im), r);
1192 else
1194 long double d, r, s;
1196 d = hypotl (re, im);
1197 /* Use the identity 2 Re res Im res = Im x
1198 to avoid cancellation error in d +/- Re x. */
1199 if (re > 0)
1201 r = sqrtl (0.5 * d + 0.5 * re);
1202 s = (0.5 * im) / r;
1204 else
1206 s = sqrtl (0.5 * d - 0.5 * re);
1207 r = fabsl ((0.5 * im) / s);
1210 COMPLEX_ASSIGN (v, r, copysignl (s, im));
1212 return v;
1214 #endif
1217 /* sinh(a + i b) = sinh(a) cos(b) + i cosh(a) sin(b) */
1218 #if !defined(HAVE_CSINHF)
1219 #define HAVE_CSINHF 1
1220 float complex csinhf (float complex a);
1222 float complex
1223 csinhf (float complex a)
1225 float r, i;
1226 float complex v;
1228 r = REALPART (a);
1229 i = IMAGPART (a);
1230 COMPLEX_ASSIGN (v, sinhf (r) * cosf (i), coshf (r) * sinf (i));
1231 return v;
1233 #endif
1235 #if !defined(HAVE_CSINH)
1236 #define HAVE_CSINH 1
1237 double complex csinh (double complex a);
1239 double complex
1240 csinh (double complex a)
1242 double r, i;
1243 double complex v;
1245 r = REALPART (a);
1246 i = IMAGPART (a);
1247 COMPLEX_ASSIGN (v, sinh (r) * cos (i), cosh (r) * sin (i));
1248 return v;
1250 #endif
1252 #if !defined(HAVE_CSINHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1253 #define HAVE_CSINHL 1
1254 long double complex csinhl (long double complex a);
1256 long double complex
1257 csinhl (long double complex a)
1259 long double r, i;
1260 long double complex v;
1262 r = REALPART (a);
1263 i = IMAGPART (a);
1264 COMPLEX_ASSIGN (v, sinhl (r) * cosl (i), coshl (r) * sinl (i));
1265 return v;
1267 #endif
1270 /* cosh(a + i b) = cosh(a) cos(b) + i sinh(a) sin(b) */
1271 #if !defined(HAVE_CCOSHF)
1272 #define HAVE_CCOSHF 1
1273 float complex ccoshf (float complex a);
1275 float complex
1276 ccoshf (float complex a)
1278 float r, i;
1279 float complex v;
1281 r = REALPART (a);
1282 i = IMAGPART (a);
1283 COMPLEX_ASSIGN (v, coshf (r) * cosf (i), sinhf (r) * sinf (i));
1284 return v;
1286 #endif
1288 #if !defined(HAVE_CCOSH)
1289 #define HAVE_CCOSH 1
1290 double complex ccosh (double complex a);
1292 double complex
1293 ccosh (double complex a)
1295 double r, i;
1296 double complex v;
1298 r = REALPART (a);
1299 i = IMAGPART (a);
1300 COMPLEX_ASSIGN (v, cosh (r) * cos (i), sinh (r) * sin (i));
1301 return v;
1303 #endif
1305 #if !defined(HAVE_CCOSHL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1306 #define HAVE_CCOSHL 1
1307 long double complex ccoshl (long double complex a);
1309 long double complex
1310 ccoshl (long double complex a)
1312 long double r, i;
1313 long double complex v;
1315 r = REALPART (a);
1316 i = IMAGPART (a);
1317 COMPLEX_ASSIGN (v, coshl (r) * cosl (i), sinhl (r) * sinl (i));
1318 return v;
1320 #endif
1323 /* tanh(a + i b) = (tanh(a) + i tan(b)) / (1 + i tanh(a) tan(b)) */
1324 #if !defined(HAVE_CTANHF)
1325 #define HAVE_CTANHF 1
1326 float complex ctanhf (float complex a);
1328 float complex
1329 ctanhf (float complex a)
1331 float rt, it;
1332 float complex n, d;
1334 rt = tanhf (REALPART (a));
1335 it = tanf (IMAGPART (a));
1336 COMPLEX_ASSIGN (n, rt, it);
1337 COMPLEX_ASSIGN (d, 1, rt * it);
1339 return n / d;
1341 #endif
1343 #if !defined(HAVE_CTANH)
1344 #define HAVE_CTANH 1
1345 double complex ctanh (double complex a);
1346 double complex
1347 ctanh (double complex a)
1349 double rt, it;
1350 double complex n, d;
1352 rt = tanh (REALPART (a));
1353 it = tan (IMAGPART (a));
1354 COMPLEX_ASSIGN (n, rt, it);
1355 COMPLEX_ASSIGN (d, 1, rt * it);
1357 return n / d;
1359 #endif
1361 #if !defined(HAVE_CTANHL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1362 #define HAVE_CTANHL 1
1363 long double complex ctanhl (long double complex a);
1365 long double complex
1366 ctanhl (long double complex a)
1368 long double rt, it;
1369 long double complex n, d;
1371 rt = tanhl (REALPART (a));
1372 it = tanl (IMAGPART (a));
1373 COMPLEX_ASSIGN (n, rt, it);
1374 COMPLEX_ASSIGN (d, 1, rt * it);
1376 return n / d;
1378 #endif
1381 /* sin(a + i b) = sin(a) cosh(b) + i cos(a) sinh(b) */
1382 #if !defined(HAVE_CSINF)
1383 #define HAVE_CSINF 1
1384 float complex csinf (float complex a);
1386 float complex
1387 csinf (float complex a)
1389 float r, i;
1390 float complex v;
1392 r = REALPART (a);
1393 i = IMAGPART (a);
1394 COMPLEX_ASSIGN (v, sinf (r) * coshf (i), cosf (r) * sinhf (i));
1395 return v;
1397 #endif
1399 #if !defined(HAVE_CSIN)
1400 #define HAVE_CSIN 1
1401 double complex csin (double complex a);
1403 double complex
1404 csin (double complex a)
1406 double r, i;
1407 double complex v;
1409 r = REALPART (a);
1410 i = IMAGPART (a);
1411 COMPLEX_ASSIGN (v, sin (r) * cosh (i), cos (r) * sinh (i));
1412 return v;
1414 #endif
1416 #if !defined(HAVE_CSINL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1417 #define HAVE_CSINL 1
1418 long double complex csinl (long double complex a);
1420 long double complex
1421 csinl (long double complex a)
1423 long double r, i;
1424 long double complex v;
1426 r = REALPART (a);
1427 i = IMAGPART (a);
1428 COMPLEX_ASSIGN (v, sinl (r) * coshl (i), cosl (r) * sinhl (i));
1429 return v;
1431 #endif
1434 /* cos(a + i b) = cos(a) cosh(b) - i sin(a) sinh(b) */
1435 #if !defined(HAVE_CCOSF)
1436 #define HAVE_CCOSF 1
1437 float complex ccosf (float complex a);
1439 float complex
1440 ccosf (float complex a)
1442 float r, i;
1443 float complex v;
1445 r = REALPART (a);
1446 i = IMAGPART (a);
1447 COMPLEX_ASSIGN (v, cosf (r) * coshf (i), - (sinf (r) * sinhf (i)));
1448 return v;
1450 #endif
1452 #if !defined(HAVE_CCOS)
1453 #define HAVE_CCOS 1
1454 double complex ccos (double complex a);
1456 double complex
1457 ccos (double complex a)
1459 double r, i;
1460 double complex v;
1462 r = REALPART (a);
1463 i = IMAGPART (a);
1464 COMPLEX_ASSIGN (v, cos (r) * cosh (i), - (sin (r) * sinh (i)));
1465 return v;
1467 #endif
1469 #if !defined(HAVE_CCOSL) && defined(HAVE_COSL) && defined(HAVE_COSHL) && defined(HAVE_SINL) && defined(HAVE_SINHL)
1470 #define HAVE_CCOSL 1
1471 long double complex ccosl (long double complex a);
1473 long double complex
1474 ccosl (long double complex a)
1476 long double r, i;
1477 long double complex v;
1479 r = REALPART (a);
1480 i = IMAGPART (a);
1481 COMPLEX_ASSIGN (v, cosl (r) * coshl (i), - (sinl (r) * sinhl (i)));
1482 return v;
1484 #endif
1487 /* tan(a + i b) = (tan(a) + i tanh(b)) / (1 - i tan(a) tanh(b)) */
1488 #if !defined(HAVE_CTANF)
1489 #define HAVE_CTANF 1
1490 float complex ctanf (float complex a);
1492 float complex
1493 ctanf (float complex a)
1495 float rt, it;
1496 float complex n, d;
1498 rt = tanf (REALPART (a));
1499 it = tanhf (IMAGPART (a));
1500 COMPLEX_ASSIGN (n, rt, it);
1501 COMPLEX_ASSIGN (d, 1, - (rt * it));
1503 return n / d;
1505 #endif
1507 #if !defined(HAVE_CTAN)
1508 #define HAVE_CTAN 1
1509 double complex ctan (double complex a);
1511 double complex
1512 ctan (double complex a)
1514 double rt, it;
1515 double complex n, d;
1517 rt = tan (REALPART (a));
1518 it = tanh (IMAGPART (a));
1519 COMPLEX_ASSIGN (n, rt, it);
1520 COMPLEX_ASSIGN (d, 1, - (rt * it));
1522 return n / d;
1524 #endif
1526 #if !defined(HAVE_CTANL) && defined(HAVE_TANL) && defined(HAVE_TANHL)
1527 #define HAVE_CTANL 1
1528 long double complex ctanl (long double complex a);
1530 long double complex
1531 ctanl (long double complex a)
1533 long double rt, it;
1534 long double complex n, d;
1536 rt = tanl (REALPART (a));
1537 it = tanhl (IMAGPART (a));
1538 COMPLEX_ASSIGN (n, rt, it);
1539 COMPLEX_ASSIGN (d, 1, - (rt * it));
1541 return n / d;
1543 #endif
1546 /* Complex ASIN. Returns wrongly NaN for infinite arguments.
1547 Algorithm taken from Abramowitz & Stegun. */
1549 #if !defined(HAVE_CASINF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1550 #define HAVE_CASINF 1
1551 complex float casinf (complex float z);
1553 complex float
1554 casinf (complex float z)
1556 return -I*clogf (I*z + csqrtf (1.0f-z*z));
1558 #endif
1561 #if !defined(HAVE_CASIN) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1562 #define HAVE_CASIN 1
1563 complex double casin (complex double z);
1565 complex double
1566 casin (complex double z)
1568 return -I*clog (I*z + csqrt (1.0-z*z));
1570 #endif
1573 #if !defined(HAVE_CASINL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1574 #define HAVE_CASINL 1
1575 complex long double casinl (complex long double z);
1577 complex long double
1578 casinl (complex long double z)
1580 return -I*clogl (I*z + csqrtl (1.0L-z*z));
1582 #endif
1585 /* Complex ACOS. Returns wrongly NaN for infinite arguments.
1586 Algorithm taken from Abramowitz & Stegun. */
1588 #if !defined(HAVE_CACOSF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1589 #define HAVE_CACOSF 1
1590 complex float cacosf (complex float z);
1592 complex float
1593 cacosf (complex float z)
1595 return -I*clogf (z + I*csqrtf (1.0f-z*z));
1597 #endif
1600 #if !defined(HAVE_CACOS) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1601 #define HAVE_CACOS 1
1602 complex double cacos (complex double z);
1604 complex double
1605 cacos (complex double z)
1607 return -I*clog (z + I*csqrt (1.0-z*z));
1609 #endif
1612 #if !defined(HAVE_CACOSL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1613 #define HAVE_CACOSL 1
1614 complex long double cacosl (complex long double z);
1616 complex long double
1617 cacosl (complex long double z)
1619 return -I*clogl (z + I*csqrtl (1.0L-z*z));
1621 #endif
1624 /* Complex ATAN. Returns wrongly NaN for infinite arguments.
1625 Algorithm taken from Abramowitz & Stegun. */
1627 #if !defined(HAVE_CATANF) && defined(HAVE_CLOGF)
1628 #define HAVE_CACOSF 1
1629 complex float catanf (complex float z);
1631 complex float
1632 catanf (complex float z)
1634 return I*clogf ((I+z)/(I-z))/2.0f;
1636 #endif
1639 #if !defined(HAVE_CATAN) && defined(HAVE_CLOG)
1640 #define HAVE_CACOS 1
1641 complex double catan (complex double z);
1643 complex double
1644 catan (complex double z)
1646 return I*clog ((I+z)/(I-z))/2.0;
1648 #endif
1651 #if !defined(HAVE_CATANL) && defined(HAVE_CLOGL)
1652 #define HAVE_CACOSL 1
1653 complex long double catanl (complex long double z);
1655 complex long double
1656 catanl (complex long double z)
1658 return I*clogl ((I+z)/(I-z))/2.0L;
1660 #endif
1663 /* Complex ASINH. Returns wrongly NaN for infinite arguments.
1664 Algorithm taken from Abramowitz & Stegun. */
1666 #if !defined(HAVE_CASINHF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1667 #define HAVE_CASINHF 1
1668 complex float casinhf (complex float z);
1670 complex float
1671 casinhf (complex float z)
1673 return clogf (z + csqrtf (z*z+1.0f));
1675 #endif
1678 #if !defined(HAVE_CASINH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1679 #define HAVE_CASINH 1
1680 complex double casinh (complex double z);
1682 complex double
1683 casinh (complex double z)
1685 return clog (z + csqrt (z*z+1.0));
1687 #endif
1690 #if !defined(HAVE_CASINHL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1691 #define HAVE_CASINHL 1
1692 complex long double casinhl (complex long double z);
1694 complex long double
1695 casinhl (complex long double z)
1697 return clogl (z + csqrtl (z*z+1.0L));
1699 #endif
1702 /* Complex ACOSH. Returns wrongly NaN for infinite arguments.
1703 Algorithm taken from Abramowitz & Stegun. */
1705 #if !defined(HAVE_CACOSHF) && defined(HAVE_CLOGF) && defined(HAVE_CSQRTF)
1706 #define HAVE_CACOSHF 1
1707 complex float cacoshf (complex float z);
1709 complex float
1710 cacoshf (complex float z)
1712 return clogf (z + csqrtf (z-1.0f) * csqrtf (z+1.0f));
1714 #endif
1717 #if !defined(HAVE_CACOSH) && defined(HAVE_CLOG) && defined(HAVE_CSQRT)
1718 #define HAVE_CACOSH 1
1719 complex double cacosh (complex double z);
1721 complex double
1722 cacosh (complex double z)
1724 return clog (z + csqrt (z-1.0) * csqrt (z+1.0));
1726 #endif
1729 #if !defined(HAVE_CACOSHL) && defined(HAVE_CLOGL) && defined(HAVE_CSQRTL)
1730 #define HAVE_CACOSHL 1
1731 complex long double cacoshl (complex long double z);
1733 complex long double
1734 cacoshl (complex long double z)
1736 return clogl (z + csqrtl (z-1.0L) * csqrtl (z+1.0L));
1738 #endif
1741 /* Complex ATANH. Returns wrongly NaN for infinite arguments.
1742 Algorithm taken from Abramowitz & Stegun. */
1744 #if !defined(HAVE_CATANHF) && defined(HAVE_CLOGF)
1745 #define HAVE_CATANHF 1
1746 complex float catanhf (complex float z);
1748 complex float
1749 catanhf (complex float z)
1751 return clogf ((1.0f+z)/(1.0f-z))/2.0f;
1753 #endif
1756 #if !defined(HAVE_CATANH) && defined(HAVE_CLOG)
1757 #define HAVE_CATANH 1
1758 complex double catanh (complex double z);
1760 complex double
1761 catanh (complex double z)
1763 return clog ((1.0+z)/(1.0-z))/2.0;
1765 #endif
1767 #if !defined(HAVE_CATANHL) && defined(HAVE_CLOGL)
1768 #define HAVE_CATANHL 1
1769 complex long double catanhl (complex long double z);
1771 complex long double
1772 catanhl (complex long double z)
1774 return clogl ((1.0L+z)/(1.0L-z))/2.0L;
1776 #endif
1779 #if !defined(HAVE_TGAMMA)
1780 #define HAVE_TGAMMA 1
1781 double tgamma (double);
1783 /* Fallback tgamma() function. Uses the algorithm from
1784 http://www.netlib.org/specfun/gamma and references therein. */
1786 #undef SQRTPI
1787 #define SQRTPI 0.9189385332046727417803297
1789 #undef PI
1790 #define PI 3.1415926535897932384626434
1792 double
1793 tgamma (double x)
1795 int i, n, parity;
1796 double fact, res, sum, xden, xnum, y, y1, ysq, z;
1798 static double p[8] = {
1799 -1.71618513886549492533811e0, 2.47656508055759199108314e1,
1800 -3.79804256470945635097577e2, 6.29331155312818442661052e2,
1801 8.66966202790413211295064e2, -3.14512729688483675254357e4,
1802 -3.61444134186911729807069e4, 6.64561438202405440627855e4 };
1804 static double q[8] = {
1805 -3.08402300119738975254353e1, 3.15350626979604161529144e2,
1806 -1.01515636749021914166146e3, -3.10777167157231109440444e3,
1807 2.25381184209801510330112e4, 4.75584627752788110767815e3,
1808 -1.34659959864969306392456e5, -1.15132259675553483497211e5 };
1810 static double c[7] = { -1.910444077728e-03,
1811 8.4171387781295e-04, -5.952379913043012e-04,
1812 7.93650793500350248e-04, -2.777777777777681622553e-03,
1813 8.333333333333333331554247e-02, 5.7083835261e-03 };
1815 static const double xminin = 2.23e-308;
1816 static const double xbig = 171.624;
1817 static const double xnan = __builtin_nan ("0x0"), xinf = __builtin_inf ();
1818 static double eps = 0;
1820 if (eps == 0)
1821 eps = nextafter (1., 2.) - 1.;
1823 parity = 0;
1824 fact = 1;
1825 n = 0;
1826 y = x;
1828 if (isnan (x))
1829 return x;
1831 if (y <= 0)
1833 y = -x;
1834 y1 = trunc (y);
1835 res = y - y1;
1837 if (res != 0)
1839 if (y1 != trunc (y1*0.5l)*2)
1840 parity = 1;
1841 fact = -PI / sin (PI*res);
1842 y = y + 1;
1844 else
1845 return x == 0 ? copysign (xinf, x) : xnan;
1848 if (y < eps)
1850 if (y >= xminin)
1851 res = 1 / y;
1852 else
1853 return xinf;
1855 else if (y < 13)
1857 y1 = y;
1858 if (y < 1)
1860 z = y;
1861 y = y + 1;
1863 else
1865 n = (int)y - 1;
1866 y = y - n;
1867 z = y - 1;
1870 xnum = 0;
1871 xden = 1;
1872 for (i = 0; i < 8; i++)
1874 xnum = (xnum + p[i]) * z;
1875 xden = xden * z + q[i];
1878 res = xnum / xden + 1;
1880 if (y1 < y)
1881 res = res / y1;
1882 else if (y1 > y)
1883 for (i = 1; i <= n; i++)
1885 res = res * y;
1886 y = y + 1;
1889 else
1891 if (y < xbig)
1893 ysq = y * y;
1894 sum = c[6];
1895 for (i = 0; i < 6; i++)
1896 sum = sum / ysq + c[i];
1898 sum = sum/y - y + SQRTPI;
1899 sum = sum + (y - 0.5) * log (y);
1900 res = exp (sum);
1902 else
1903 return x < 0 ? xnan : xinf;
1906 if (parity)
1907 res = -res;
1908 if (fact != 1)
1909 res = fact / res;
1911 return res;
1913 #endif
1917 #if !defined(HAVE_LGAMMA)
1918 #define HAVE_LGAMMA 1
1919 double lgamma (double);
1921 /* Fallback lgamma() function. Uses the algorithm from
1922 http://www.netlib.org/specfun/algama and references therein,
1923 except for negative arguments (where netlib would return +Inf)
1924 where we use the following identity:
1925 lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y)
1928 double
1929 lgamma (double y)
1932 #undef SQRTPI
1933 #define SQRTPI 0.9189385332046727417803297
1935 #undef PI
1936 #define PI 3.1415926535897932384626434
1938 #define PNT68 0.6796875
1939 #define D1 -0.5772156649015328605195174
1940 #define D2 0.4227843350984671393993777
1941 #define D4 1.791759469228055000094023
1943 static double p1[8] = {
1944 4.945235359296727046734888e0, 2.018112620856775083915565e2,
1945 2.290838373831346393026739e3, 1.131967205903380828685045e4,
1946 2.855724635671635335736389e4, 3.848496228443793359990269e4,
1947 2.637748787624195437963534e4, 7.225813979700288197698961e3 };
1948 static double q1[8] = {
1949 6.748212550303777196073036e1, 1.113332393857199323513008e3,
1950 7.738757056935398733233834e3, 2.763987074403340708898585e4,
1951 5.499310206226157329794414e4, 6.161122180066002127833352e4,
1952 3.635127591501940507276287e4, 8.785536302431013170870835e3 };
1953 static double p2[8] = {
1954 4.974607845568932035012064e0, 5.424138599891070494101986e2,
1955 1.550693864978364947665077e4, 1.847932904445632425417223e5,
1956 1.088204769468828767498470e6, 3.338152967987029735917223e6,
1957 5.106661678927352456275255e6, 3.074109054850539556250927e6 };
1958 static double q2[8] = {
1959 1.830328399370592604055942e2, 7.765049321445005871323047e3,
1960 1.331903827966074194402448e5, 1.136705821321969608938755e6,
1961 5.267964117437946917577538e6, 1.346701454311101692290052e7,
1962 1.782736530353274213975932e7, 9.533095591844353613395747e6 };
1963 static double p4[8] = {
1964 1.474502166059939948905062e4, 2.426813369486704502836312e6,
1965 1.214755574045093227939592e8, 2.663432449630976949898078e9,
1966 2.940378956634553899906876e10, 1.702665737765398868392998e11,
1967 4.926125793377430887588120e11, 5.606251856223951465078242e11 };
1968 static double q4[8] = {
1969 2.690530175870899333379843e3, 6.393885654300092398984238e5,
1970 4.135599930241388052042842e7, 1.120872109616147941376570e9,
1971 1.488613728678813811542398e10, 1.016803586272438228077304e11,
1972 3.417476345507377132798597e11, 4.463158187419713286462081e11 };
1973 static double c[7] = {
1974 -1.910444077728e-03, 8.4171387781295e-04,
1975 -5.952379913043012e-04, 7.93650793500350248e-04,
1976 -2.777777777777681622553e-03, 8.333333333333333331554247e-02,
1977 5.7083835261e-03 };
1979 static double xbig = 2.55e305, xinf = __builtin_inf (), eps = 0,
1980 frtbig = 2.25e76;
1982 int i;
1983 double corr, res, xden, xm1, xm2, xm4, xnum, ysq;
1985 if (eps == 0)
1986 eps = __builtin_nextafter (1., 2.) - 1.;
1988 if ((y > 0) && (y <= xbig))
1990 if (y <= eps)
1991 res = -log (y);
1992 else if (y <= 1.5)
1994 if (y < PNT68)
1996 corr = -log (y);
1997 xm1 = y;
1999 else
2001 corr = 0;
2002 xm1 = (y - 0.5) - 0.5;
2005 if ((y <= 0.5) || (y >= PNT68))
2007 xden = 1;
2008 xnum = 0;
2009 for (i = 0; i < 8; i++)
2011 xnum = xnum*xm1 + p1[i];
2012 xden = xden*xm1 + q1[i];
2014 res = corr + (xm1 * (D1 + xm1*(xnum/xden)));
2016 else
2018 xm2 = (y - 0.5) - 0.5;
2019 xden = 1;
2020 xnum = 0;
2021 for (i = 0; i < 8; i++)
2023 xnum = xnum*xm2 + p2[i];
2024 xden = xden*xm2 + q2[i];
2026 res = corr + xm2 * (D2 + xm2*(xnum/xden));
2029 else if (y <= 4)
2031 xm2 = y - 2;
2032 xden = 1;
2033 xnum = 0;
2034 for (i = 0; i < 8; i++)
2036 xnum = xnum*xm2 + p2[i];
2037 xden = xden*xm2 + q2[i];
2039 res = xm2 * (D2 + xm2*(xnum/xden));
2041 else if (y <= 12)
2043 xm4 = y - 4;
2044 xden = -1;
2045 xnum = 0;
2046 for (i = 0; i < 8; i++)
2048 xnum = xnum*xm4 + p4[i];
2049 xden = xden*xm4 + q4[i];
2051 res = D4 + xm4*(xnum/xden);
2053 else
2055 res = 0;
2056 if (y <= frtbig)
2058 res = c[6];
2059 ysq = y * y;
2060 for (i = 0; i < 6; i++)
2061 res = res / ysq + c[i];
2063 res = res/y;
2064 corr = log (y);
2065 res = res + SQRTPI - 0.5*corr;
2066 res = res + y*(corr-1);
2069 else if (y < 0 && __builtin_floor (y) != y)
2071 /* lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y)
2072 For abs(y) very close to zero, we use a series expansion to
2073 the first order in y to avoid overflow. */
2074 if (y > -1.e-100)
2075 res = -2 * log (fabs (y)) - lgamma (-y);
2076 else
2077 res = log (PI / fabs (y * sin (PI * y))) - lgamma (-y);
2079 else
2080 res = xinf;
2082 return res;
2084 #endif
2087 #if defined(HAVE_TGAMMA) && !defined(HAVE_TGAMMAF)
2088 #define HAVE_TGAMMAF 1
2089 float tgammaf (float);
2091 float
2092 tgammaf (float x)
2094 return (float) tgamma ((double) x);
2096 #endif
2098 #if defined(HAVE_LGAMMA) && !defined(HAVE_LGAMMAF)
2099 #define HAVE_LGAMMAF 1
2100 float lgammaf (float);
2102 float
2103 lgammaf (float x)
2105 return (float) lgamma ((double) x);
2107 #endif