re PR libfortran/33583 (FAIL: gfortran.dg/gamma_1.f90)

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.

From-SVN: r130245
This commit is contained in:
Francois-Xavier Coudert 2007-11-16 22:31:28 +00:00 committed by François-Xavier Coudert
parent 4a58b9ad36
commit fb0a0e1591
6 changed files with 852 additions and 0 deletions

View File

@ -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.

View File

@ -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 */

View File

@ -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

462
libgfortran/configure vendored
View File

@ -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

View File

@ -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])])

View File

@ -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