diff options
-rw-r--r-- | libgfortran/ChangeLog | 12 | ||||
-rw-r--r-- | libgfortran/c99_protos.h | 22 | ||||
-rw-r--r-- | libgfortran/config.h.in | 18 | ||||
-rwxr-xr-x | libgfortran/configure | 462 | ||||
-rw-r--r-- | libgfortran/configure.ac | 6 | ||||
-rw-r--r-- | libgfortran/intrinsics/c99_functions.c | 332 |
6 files changed, 852 insertions, 0 deletions
diff --git a/libgfortran/ChangeLog b/libgfortran/ChangeLog index 801f04d..0fecb6f 100644 --- a/libgfortran/ChangeLog +++ b/libgfortran/ChangeLog @@ -1,3 +1,15 @@ +2007-11-16 Francois-Xavier Coudert <fxcoudert@gcc.gnu.org> + + PR libfortran/33583 + PR libfortran/33698 + * intrinsics/c99_functions.c (tgamma, tgammaf, lgamma, lgammaf): + New fallback functions. + * c99_protos.h (tgamma, tgammaf, lgamma, lgammaf): New prototypes. + * configure.ac: Add checks for tgamma, tgammaf, tgammal, lgamma, + lgammaf and lgammal. + * config.h.in: Regenerate. + * configure: Regenerate. + 2007-11-08 Francois-Xavier Coudert <fxcoudert@gcc.gnu.org> * mk-kinds-h.sh: Change sed syntax. diff --git a/libgfortran/c99_protos.h b/libgfortran/c99_protos.h index eeeb6d8..59cbe4c 100644 --- a/libgfortran/c99_protos.h +++ b/libgfortran/c99_protos.h @@ -502,5 +502,27 @@ extern long double complex ctanl (long double complex); #endif +/* Gamma-related prototypes. */ +#if !defined(HAVE_TGAMMA) +#define HAVE_TGAMMA 1 +extern double tgamma (double); +#endif + +#if !defined(HAVE_LGAMMA) +#define HAVE_LGAMMA 1 +extern double lgamma (double); +#endif + +#if defined(HAVE_TGAMMA) && !defined(HAVE_TGAMMAF) +#define HAVE_TGAMMAF 1 +extern float tgammaf (float); +#endif + +#if defined(HAVE_LGAMMA) && !defined(HAVE_LGAMMAF) +#define HAVE_LGAMMAF 1 +extern float lgammaf (float); +#endif + + #endif /* C99_PROTOS_H */ diff --git a/libgfortran/config.h.in b/libgfortran/config.h.in index ffa419c..72d46d9 100644 --- a/libgfortran/config.h.in +++ b/libgfortran/config.h.in @@ -483,6 +483,15 @@ /* libm includes ldexpl */ #undef HAVE_LDEXPL +/* libm includes lgamma */ +#undef HAVE_LGAMMA + +/* libm includes lgammaf */ +#undef HAVE_LGAMMAF + +/* libm includes lgammal */ +#undef HAVE_LGAMMAL + /* Define to 1 if you have the `link' function. */ #undef HAVE_LINK @@ -705,6 +714,15 @@ /* libm includes tanl */ #undef HAVE_TANL +/* libm includes tgamma */ +#undef HAVE_TGAMMA + +/* libm includes tgammaf */ +#undef HAVE_TGAMMAF + +/* libm includes tgammal */ +#undef HAVE_TGAMMAL + /* Define to 1 if you have the `time' function. */ #undef HAVE_TIME diff --git a/libgfortran/configure b/libgfortran/configure index 5132c28..f33516d 100755 --- a/libgfortran/configure +++ b/libgfortran/configure @@ -31425,6 +31425,468 @@ _ACEOF fi +echo "$as_me:$LINENO: checking for tgamma in -lm" >&5 +echo $ECHO_N "checking for tgamma in -lm... $ECHO_C" >&6 +if test "${ac_cv_lib_m_tgamma+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + ac_check_lib_save_LIBS=$LIBS +LIBS="-lm $LIBS" +if test x$gcc_no_link = xyes; then + { { echo "$as_me:$LINENO: error: Link tests are not allowed after GCC_NO_EXECUTABLES." >&5 +echo "$as_me: error: Link tests are not allowed after GCC_NO_EXECUTABLES." >&2;} + { (exit 1); exit 1; }; } +fi +cat >conftest.$ac_ext <<_ACEOF +/* confdefs.h. */ +_ACEOF +cat confdefs.h >>conftest.$ac_ext +cat >>conftest.$ac_ext <<_ACEOF +/* end confdefs.h. */ + +/* Override any gcc2 internal prototype to avoid an error. */ +#ifdef __cplusplus +extern "C" +#endif +/* We use char because int might match the return type of a gcc2 + builtin and then its argument prototype would still apply. */ +char tgamma (); +int +main () +{ +tgamma (); + ; + return 0; +} +_ACEOF +rm -f conftest.$ac_objext conftest$ac_exeext +if { (eval echo "$as_me:$LINENO: \"$ac_link\"") >&5 + (eval $ac_link) 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && + { ac_try='test -z "$ac_c_werror_flag" + || test ! -s conftest.err' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; } && + { ac_try='test -s conftest$ac_exeext' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; }; then + ac_cv_lib_m_tgamma=yes +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + +ac_cv_lib_m_tgamma=no +fi +rm -f conftest.err conftest.$ac_objext \ + conftest$ac_exeext conftest.$ac_ext +LIBS=$ac_check_lib_save_LIBS +fi +echo "$as_me:$LINENO: result: $ac_cv_lib_m_tgamma" >&5 +echo "${ECHO_T}$ac_cv_lib_m_tgamma" >&6 +if test $ac_cv_lib_m_tgamma = yes; then + +cat >>confdefs.h <<\_ACEOF +#define HAVE_TGAMMA 1 +_ACEOF + +fi + +echo "$as_me:$LINENO: checking for tgammaf in -lm" >&5 +echo $ECHO_N "checking for tgammaf in -lm... $ECHO_C" >&6 +if test "${ac_cv_lib_m_tgammaf+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + ac_check_lib_save_LIBS=$LIBS +LIBS="-lm $LIBS" +if test x$gcc_no_link = xyes; then + { { echo "$as_me:$LINENO: error: Link tests are not allowed after GCC_NO_EXECUTABLES." >&5 +echo "$as_me: error: Link tests are not allowed after GCC_NO_EXECUTABLES." >&2;} + { (exit 1); exit 1; }; } +fi +cat >conftest.$ac_ext <<_ACEOF +/* confdefs.h. */ +_ACEOF +cat confdefs.h >>conftest.$ac_ext +cat >>conftest.$ac_ext <<_ACEOF +/* end confdefs.h. */ + +/* Override any gcc2 internal prototype to avoid an error. */ +#ifdef __cplusplus +extern "C" +#endif +/* We use char because int might match the return type of a gcc2 + builtin and then its argument prototype would still apply. */ +char tgammaf (); +int +main () +{ +tgammaf (); + ; + return 0; +} +_ACEOF +rm -f conftest.$ac_objext conftest$ac_exeext +if { (eval echo "$as_me:$LINENO: \"$ac_link\"") >&5 + (eval $ac_link) 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && + { ac_try='test -z "$ac_c_werror_flag" + || test ! -s conftest.err' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; } && + { ac_try='test -s conftest$ac_exeext' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; }; then + ac_cv_lib_m_tgammaf=yes +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + +ac_cv_lib_m_tgammaf=no +fi +rm -f conftest.err conftest.$ac_objext \ + conftest$ac_exeext conftest.$ac_ext +LIBS=$ac_check_lib_save_LIBS +fi +echo "$as_me:$LINENO: result: $ac_cv_lib_m_tgammaf" >&5 +echo "${ECHO_T}$ac_cv_lib_m_tgammaf" >&6 +if test $ac_cv_lib_m_tgammaf = yes; then + +cat >>confdefs.h <<\_ACEOF +#define HAVE_TGAMMAF 1 +_ACEOF + +fi + +echo "$as_me:$LINENO: checking for tgammal in -lm" >&5 +echo $ECHO_N "checking for tgammal in -lm... $ECHO_C" >&6 +if test "${ac_cv_lib_m_tgammal+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + ac_check_lib_save_LIBS=$LIBS +LIBS="-lm $LIBS" +if test x$gcc_no_link = xyes; then + { { echo "$as_me:$LINENO: error: Link tests are not allowed after GCC_NO_EXECUTABLES." >&5 +echo "$as_me: error: Link tests are not allowed after GCC_NO_EXECUTABLES." >&2;} + { (exit 1); exit 1; }; } +fi +cat >conftest.$ac_ext <<_ACEOF +/* confdefs.h. */ +_ACEOF +cat confdefs.h >>conftest.$ac_ext +cat >>conftest.$ac_ext <<_ACEOF +/* end confdefs.h. */ + +/* Override any gcc2 internal prototype to avoid an error. */ +#ifdef __cplusplus +extern "C" +#endif +/* We use char because int might match the return type of a gcc2 + builtin and then its argument prototype would still apply. */ +char tgammal (); +int +main () +{ +tgammal (); + ; + return 0; +} +_ACEOF +rm -f conftest.$ac_objext conftest$ac_exeext +if { (eval echo "$as_me:$LINENO: \"$ac_link\"") >&5 + (eval $ac_link) 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && + { ac_try='test -z "$ac_c_werror_flag" + || test ! -s conftest.err' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; } && + { ac_try='test -s conftest$ac_exeext' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; }; then + ac_cv_lib_m_tgammal=yes +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + +ac_cv_lib_m_tgammal=no +fi +rm -f conftest.err conftest.$ac_objext \ + conftest$ac_exeext conftest.$ac_ext +LIBS=$ac_check_lib_save_LIBS +fi +echo "$as_me:$LINENO: result: $ac_cv_lib_m_tgammal" >&5 +echo "${ECHO_T}$ac_cv_lib_m_tgammal" >&6 +if test $ac_cv_lib_m_tgammal = yes; then + +cat >>confdefs.h <<\_ACEOF +#define HAVE_TGAMMAL 1 +_ACEOF + +fi + +echo "$as_me:$LINENO: checking for lgamma in -lm" >&5 +echo $ECHO_N "checking for lgamma in -lm... $ECHO_C" >&6 +if test "${ac_cv_lib_m_lgamma+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + ac_check_lib_save_LIBS=$LIBS +LIBS="-lm $LIBS" +if test x$gcc_no_link = xyes; then + { { echo "$as_me:$LINENO: error: Link tests are not allowed after GCC_NO_EXECUTABLES." >&5 +echo "$as_me: error: Link tests are not allowed after GCC_NO_EXECUTABLES." >&2;} + { (exit 1); exit 1; }; } +fi +cat >conftest.$ac_ext <<_ACEOF +/* confdefs.h. */ +_ACEOF +cat confdefs.h >>conftest.$ac_ext +cat >>conftest.$ac_ext <<_ACEOF +/* end confdefs.h. */ + +/* Override any gcc2 internal prototype to avoid an error. */ +#ifdef __cplusplus +extern "C" +#endif +/* We use char because int might match the return type of a gcc2 + builtin and then its argument prototype would still apply. */ +char lgamma (); +int +main () +{ +lgamma (); + ; + return 0; +} +_ACEOF +rm -f conftest.$ac_objext conftest$ac_exeext +if { (eval echo "$as_me:$LINENO: \"$ac_link\"") >&5 + (eval $ac_link) 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && + { ac_try='test -z "$ac_c_werror_flag" + || test ! -s conftest.err' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; } && + { ac_try='test -s conftest$ac_exeext' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; }; then + ac_cv_lib_m_lgamma=yes +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + +ac_cv_lib_m_lgamma=no +fi +rm -f conftest.err conftest.$ac_objext \ + conftest$ac_exeext conftest.$ac_ext +LIBS=$ac_check_lib_save_LIBS +fi +echo "$as_me:$LINENO: result: $ac_cv_lib_m_lgamma" >&5 +echo "${ECHO_T}$ac_cv_lib_m_lgamma" >&6 +if test $ac_cv_lib_m_lgamma = yes; then + +cat >>confdefs.h <<\_ACEOF +#define HAVE_LGAMMA 1 +_ACEOF + +fi + +echo "$as_me:$LINENO: checking for lgammaf in -lm" >&5 +echo $ECHO_N "checking for lgammaf in -lm... $ECHO_C" >&6 +if test "${ac_cv_lib_m_lgammaf+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + ac_check_lib_save_LIBS=$LIBS +LIBS="-lm $LIBS" +if test x$gcc_no_link = xyes; then + { { echo "$as_me:$LINENO: error: Link tests are not allowed after GCC_NO_EXECUTABLES." >&5 +echo "$as_me: error: Link tests are not allowed after GCC_NO_EXECUTABLES." >&2;} + { (exit 1); exit 1; }; } +fi +cat >conftest.$ac_ext <<_ACEOF +/* confdefs.h. */ +_ACEOF +cat confdefs.h >>conftest.$ac_ext +cat >>conftest.$ac_ext <<_ACEOF +/* end confdefs.h. */ + +/* Override any gcc2 internal prototype to avoid an error. */ +#ifdef __cplusplus +extern "C" +#endif +/* We use char because int might match the return type of a gcc2 + builtin and then its argument prototype would still apply. */ +char lgammaf (); +int +main () +{ +lgammaf (); + ; + return 0; +} +_ACEOF +rm -f conftest.$ac_objext conftest$ac_exeext +if { (eval echo "$as_me:$LINENO: \"$ac_link\"") >&5 + (eval $ac_link) 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && + { ac_try='test -z "$ac_c_werror_flag" + || test ! -s conftest.err' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; } && + { ac_try='test -s conftest$ac_exeext' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; }; then + ac_cv_lib_m_lgammaf=yes +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + +ac_cv_lib_m_lgammaf=no +fi +rm -f conftest.err conftest.$ac_objext \ + conftest$ac_exeext conftest.$ac_ext +LIBS=$ac_check_lib_save_LIBS +fi +echo "$as_me:$LINENO: result: $ac_cv_lib_m_lgammaf" >&5 +echo "${ECHO_T}$ac_cv_lib_m_lgammaf" >&6 +if test $ac_cv_lib_m_lgammaf = yes; then + +cat >>confdefs.h <<\_ACEOF +#define HAVE_LGAMMAF 1 +_ACEOF + +fi + +echo "$as_me:$LINENO: checking for lgammal in -lm" >&5 +echo $ECHO_N "checking for lgammal in -lm... $ECHO_C" >&6 +if test "${ac_cv_lib_m_lgammal+set}" = set; then + echo $ECHO_N "(cached) $ECHO_C" >&6 +else + ac_check_lib_save_LIBS=$LIBS +LIBS="-lm $LIBS" +if test x$gcc_no_link = xyes; then + { { echo "$as_me:$LINENO: error: Link tests are not allowed after GCC_NO_EXECUTABLES." >&5 +echo "$as_me: error: Link tests are not allowed after GCC_NO_EXECUTABLES." >&2;} + { (exit 1); exit 1; }; } +fi +cat >conftest.$ac_ext <<_ACEOF +/* confdefs.h. */ +_ACEOF +cat confdefs.h >>conftest.$ac_ext +cat >>conftest.$ac_ext <<_ACEOF +/* end confdefs.h. */ + +/* Override any gcc2 internal prototype to avoid an error. */ +#ifdef __cplusplus +extern "C" +#endif +/* We use char because int might match the return type of a gcc2 + builtin and then its argument prototype would still apply. */ +char lgammal (); +int +main () +{ +lgammal (); + ; + return 0; +} +_ACEOF +rm -f conftest.$ac_objext conftest$ac_exeext +if { (eval echo "$as_me:$LINENO: \"$ac_link\"") >&5 + (eval $ac_link) 2>conftest.er1 + ac_status=$? + grep -v '^ *+' conftest.er1 >conftest.err + rm -f conftest.er1 + cat conftest.err >&5 + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); } && + { ac_try='test -z "$ac_c_werror_flag" + || test ! -s conftest.err' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; } && + { ac_try='test -s conftest$ac_exeext' + { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5 + (eval $ac_try) 2>&5 + ac_status=$? + echo "$as_me:$LINENO: \$? = $ac_status" >&5 + (exit $ac_status); }; }; then + ac_cv_lib_m_lgammal=yes +else + echo "$as_me: failed program was:" >&5 +sed 's/^/| /' conftest.$ac_ext >&5 + +ac_cv_lib_m_lgammal=no +fi +rm -f conftest.err conftest.$ac_objext \ + conftest$ac_exeext conftest.$ac_ext +LIBS=$ac_check_lib_save_LIBS +fi +echo "$as_me:$LINENO: result: $ac_cv_lib_m_lgammal" >&5 +echo "${ECHO_T}$ac_cv_lib_m_lgammal" >&6 +if test $ac_cv_lib_m_lgammal = yes; then + +cat >>confdefs.h <<\_ACEOF +#define HAVE_LGAMMAL 1 +_ACEOF + +fi + # On AIX, clog is present in libm as __clog echo "$as_me:$LINENO: checking for __clog in -lm" >&5 diff --git a/libgfortran/configure.ac b/libgfortran/configure.ac index 0d153d4..a1caf3b 100644 --- a/libgfortran/configure.ac +++ b/libgfortran/configure.ac @@ -379,6 +379,12 @@ AC_CHECK_LIB([m],[y1l],[AC_DEFINE([HAVE_Y1L],[1],[libm includes y1l])]) AC_CHECK_LIB([m],[ynf],[AC_DEFINE([HAVE_YNF],[1],[libm includes ynf])]) AC_CHECK_LIB([m],[yn],[AC_DEFINE([HAVE_YN],[1],[libm includes yn])]) AC_CHECK_LIB([m],[ynl],[AC_DEFINE([HAVE_YNL],[1],[libm includes ynl])]) +AC_CHECK_LIB([m],[tgamma],[AC_DEFINE([HAVE_TGAMMA],[1],[libm includes tgamma])]) +AC_CHECK_LIB([m],[tgammaf],[AC_DEFINE([HAVE_TGAMMAF],[1],[libm includes tgammaf])]) +AC_CHECK_LIB([m],[tgammal],[AC_DEFINE([HAVE_TGAMMAL],[1],[libm includes tgammal])]) +AC_CHECK_LIB([m],[lgamma],[AC_DEFINE([HAVE_LGAMMA],[1],[libm includes lgamma])]) +AC_CHECK_LIB([m],[lgammaf],[AC_DEFINE([HAVE_LGAMMAF],[1],[libm includes lgammaf])]) +AC_CHECK_LIB([m],[lgammal],[AC_DEFINE([HAVE_LGAMMAL],[1],[libm includes lgammal])]) # On AIX, clog is present in libm as __clog AC_CHECK_LIB([m],[__clog],[AC_DEFINE([HAVE_CLOG],[1],[libm includes clog])]) diff --git a/libgfortran/intrinsics/c99_functions.c b/libgfortran/intrinsics/c99_functions.c index c9c4796..13d5503 100644 --- a/libgfortran/intrinsics/c99_functions.c +++ b/libgfortran/intrinsics/c99_functions.c @@ -1414,3 +1414,335 @@ ctanl (long double complex a) } #endif + +#if !defined(HAVE_TGAMMA) +#define HAVE_TGAMMA 1 + +extern double tgamma (double); + +/* Fallback tgamma() function. Uses the algorithm from + http://www.netlib.org/specfun/gamma and references therein. */ + +#undef SQRTPI +#define SQRTPI 0.9189385332046727417803297 + +#undef PI +#define PI 3.1415926535897932384626434 + +double +tgamma (double x) +{ + int i, n, parity; + double fact, res, sum, xden, xnum, y, y1, ysq, z; + + static double p[8] = { + -1.71618513886549492533811e0, 2.47656508055759199108314e1, + -3.79804256470945635097577e2, 6.29331155312818442661052e2, + 8.66966202790413211295064e2, -3.14512729688483675254357e4, + -3.61444134186911729807069e4, 6.64561438202405440627855e4 }; + + static double q[8] = { + -3.08402300119738975254353e1, 3.15350626979604161529144e2, + -1.01515636749021914166146e3, -3.10777167157231109440444e3, + 2.25381184209801510330112e4, 4.75584627752788110767815e3, + -1.34659959864969306392456e5, -1.15132259675553483497211e5 }; + + static double c[7] = { -1.910444077728e-03, + 8.4171387781295e-04, -5.952379913043012e-04, + 7.93650793500350248e-04, -2.777777777777681622553e-03, + 8.333333333333333331554247e-02, 5.7083835261e-03 }; + + static const double xminin = 2.23e-308; + static const double xbig = 171.624; + static const double xnan = __builtin_nan ("0x0"), xinf = __builtin_inf (); + static double eps = 0; + + if (eps == 0) + eps = nextafter(1., 2.) - 1.; + + parity = 0; + fact = 1; + n = 0; + y = x; + + if (__builtin_isnan (x)) + return x; + + if (y <= 0) + { + y = -x; + y1 = trunc(y); + res = y - y1; + + if (res != 0) + { + if (y1 != trunc(y1*0.5l)*2) + parity = 1; + fact = -PI / sin(PI*res); + y = y + 1; + } + else + return x == 0 ? copysign (xinf, x) : xnan; + } + + if (y < eps) + { + if (y >= xminin) + res = 1 / y; + else + return xinf; + } + else if (y < 13) + { + y1 = y; + if (y < 1) + { + z = y; + y = y + 1; + } + else + { + n = (int)y - 1; + y = y - n; + z = y - 1; + } + + xnum = 0; + xden = 1; + for (i = 0; i < 8; i++) + { + xnum = (xnum + p[i]) * z; + xden = xden * z + q[i]; + } + + res = xnum / xden + 1; + + if (y1 < y) + res = res / y1; + else if (y1 > y) + for (i = 1; i <= n; i++) + { + res = res * y; + y = y + 1; + } + } + else + { + if (y < xbig) + { + ysq = y * y; + sum = c[6]; + for (i = 0; i < 6; i++) + sum = sum / ysq + c[i]; + + sum = sum/y - y + SQRTPI; + sum = sum + (y - 0.5) * log(y); + res = exp(sum); + } + else + return x < 0 ? xnan : xinf; + } + + if (parity) + res = -res; + if (fact != 1) + res = fact / res; + + return res; +} +#endif + + + +#if !defined(HAVE_LGAMMA) +#define HAVE_LGAMMA 1 + +extern double lgamma (double); + +/* Fallback lgamma() function. Uses the algorithm from + http://www.netlib.org/specfun/algama and references therein, + except for negative arguments (where netlib would return +Inf) + where we use the following identity: + lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y) + */ + +double +lgamma (double y) +{ + +#undef SQRTPI +#define SQRTPI 0.9189385332046727417803297 + +#undef PI +#define PI 3.1415926535897932384626434 + +#define PNT68 0.6796875 +#define D1 -0.5772156649015328605195174 +#define D2 0.4227843350984671393993777 +#define D4 1.791759469228055000094023 + + static double p1[8] = { + 4.945235359296727046734888e0, 2.018112620856775083915565e2, + 2.290838373831346393026739e3, 1.131967205903380828685045e4, + 2.855724635671635335736389e4, 3.848496228443793359990269e4, + 2.637748787624195437963534e4, 7.225813979700288197698961e3 }; + static double q1[8] = { + 6.748212550303777196073036e1, 1.113332393857199323513008e3, + 7.738757056935398733233834e3, 2.763987074403340708898585e4, + 5.499310206226157329794414e4, 6.161122180066002127833352e4, + 3.635127591501940507276287e4, 8.785536302431013170870835e3 }; + static double p2[8] = { + 4.974607845568932035012064e0, 5.424138599891070494101986e2, + 1.550693864978364947665077e4, 1.847932904445632425417223e5, + 1.088204769468828767498470e6, 3.338152967987029735917223e6, + 5.106661678927352456275255e6, 3.074109054850539556250927e6 }; + static double q2[8] = { + 1.830328399370592604055942e2, 7.765049321445005871323047e3, + 1.331903827966074194402448e5, 1.136705821321969608938755e6, + 5.267964117437946917577538e6, 1.346701454311101692290052e7, + 1.782736530353274213975932e7, 9.533095591844353613395747e6 }; + static double p4[8] = { + 1.474502166059939948905062e4, 2.426813369486704502836312e6, + 1.214755574045093227939592e8, 2.663432449630976949898078e9, + 2.940378956634553899906876e10, 1.702665737765398868392998e11, + 4.926125793377430887588120e11, 5.606251856223951465078242e11 }; + static double q4[8] = { + 2.690530175870899333379843e3, 6.393885654300092398984238e5, + 4.135599930241388052042842e7, 1.120872109616147941376570e9, + 1.488613728678813811542398e10, 1.016803586272438228077304e11, + 3.417476345507377132798597e11, 4.463158187419713286462081e11 }; + static double c[7] = { + -1.910444077728e-03, 8.4171387781295e-04, + -5.952379913043012e-04, 7.93650793500350248e-04, + -2.777777777777681622553e-03, 8.333333333333333331554247e-02, + 5.7083835261e-03 }; + + static double xbig = 2.55e305, xinf = __builtin_inf (), eps = 0, + frtbig = 2.25e76; + + int i; + double corr, res, xden, xm1, xm2, xm4, xnum, ysq; + + if (eps == 0) + eps = __builtin_nextafter(1., 2.) - 1.; + + if ((y > 0) && (y <= xbig)) + { + if (y <= eps) + res = -log(y); + else if (y <= 1.5) + { + if (y < PNT68) + { + corr = -log(y); + xm1 = y; + } + else + { + corr = 0; + xm1 = (y - 0.5) - 0.5; + } + + if ((y <= 0.5) || (y >= PNT68)) + { + xden = 1; + xnum = 0; + for (i = 0; i < 8; i++) + { + xnum = xnum*xm1 + p1[i]; + xden = xden*xm1 + q1[i]; + } + res = corr + (xm1 * (D1 + xm1*(xnum/xden))); + } + else + { + xm2 = (y - 0.5) - 0.5; + xden = 1; + xnum = 0; + for (i = 0; i < 8; i++) + { + xnum = xnum*xm2 + p2[i]; + xden = xden*xm2 + q2[i]; + } + res = corr + xm2 * (D2 + xm2*(xnum/xden)); + } + } + else if (y <= 4) + { + xm2 = y - 2; + xden = 1; + xnum = 0; + for (i = 0; i < 8; i++) + { + xnum = xnum*xm2 + p2[i]; + xden = xden*xm2 + q2[i]; + } + res = xm2 * (D2 + xm2*(xnum/xden)); + } + else if (y <= 12) + { + xm4 = y - 4; + xden = -1; + xnum = 0; + for (i = 0; i < 8; i++) + { + xnum = xnum*xm4 + p4[i]; + xden = xden*xm4 + q4[i]; + } + res = D4 + xm4*(xnum/xden); + } + else + { + res = 0; + if (y <= frtbig) + { + res = c[6]; + ysq = y * y; + for (i = 0; i < 6; i++) + res = res / ysq + c[i]; + } + res = res/y; + corr = log(y); + res = res + SQRTPI - 0.5*corr; + res = res + y*(corr-1); + } + } + else if (y < 0 && __builtin_floor (y) != y) + { + /* lgamma(y) = log(pi/(|y*sin(pi*y)|)) - lgamma(-y) + For abs(y) very close to zero, we use a series expansion to + the first order in y to avoid overflow. */ + if (y > -1.e-100) + res = -2 * log (fabs (y)) - lgamma (-y); + else + res = log (PI / fabs (y * sin (PI * y))) - lgamma (-y); + } + else + res = xinf; + + return res; +} +#endif + + +#if defined(HAVE_TGAMMA) && !defined(HAVE_TGAMMAF) +#define HAVE_TGAMMAF 1 +extern float tgammaf (float); + +float +tgammaf (float x) +{ + return (float) tgamma ((double) x); +} +#endif + +#if defined(HAVE_LGAMMA) && !defined(HAVE_LGAMMAF) +#define HAVE_LGAMMAF 1 +extern float lgammaf (float); + +float +lgammaf (float x) +{ + return (float) lgamma ((double) x); +} +#endif |