Issue #3366: Add lgamma function to math module.
diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c
index 350034f..b1e4521 100644
--- a/Modules/mathmodule.c
+++ b/Modules/mathmodule.c
@@ -322,6 +322,60 @@
}
/*
+ lgamma: natural log of the absolute value of the Gamma function.
+ For large arguments, Lanczos' formula works extremely well here.
+*/
+
+static double
+m_lgamma(double x)
+{
+ double r, absx;
+
+ /* special cases */
+ if (!Py_IS_FINITE(x)) {
+ if (Py_IS_NAN(x))
+ return x; /* lgamma(nan) = nan */
+ else
+ return Py_HUGE_VAL; /* lgamma(+-inf) = +inf */
+ }
+
+ /* integer arguments */
+ if (x == floor(x) && x <= 2.0) {
+ if (x <= 0.0) {
+ errno = EDOM; /* lgamma(n) = inf, divide-by-zero for */
+ return Py_HUGE_VAL; /* integers n <= 0 */
+ }
+ else {
+ return 0.0; /* lgamma(1) = lgamma(2) = 0.0 */
+ }
+ }
+
+ absx = fabs(x);
+ /* tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x */
+ if (absx < 1e-20)
+ return -log(absx);
+
+ /* Lanczos' formula */
+ if (x > 0.0) {
+ /* we could save a fraction of a ulp in accuracy by having a
+ second set of numerator coefficients for lanczos_sum that
+ absorbed the exp(-lanczos_g) term, and throwing out the
+ lanczos_g subtraction below; it's probably not worth it. */
+ r = log(lanczos_sum(x)) - lanczos_g +
+ (x-0.5)*(log(x+lanczos_g-0.5)-1);
+ }
+ else {
+ r = log(pi) - log(fabs(sinpi(absx))) - log(absx) -
+ (log(lanczos_sum(absx)) - lanczos_g +
+ (absx-0.5)*(log(absx+lanczos_g-0.5)-1));
+ }
+ if (Py_IS_INFINITY(r))
+ errno = ERANGE;
+ return r;
+}
+
+
+/*
wrapper for atan2 that deals directly with special cases before
delegating to the platform libm for the remaining cases. This
is necessary to get consistent behaviour across platforms.
@@ -639,6 +693,8 @@
"This is the largest integral value <= x.")
FUNC1A(gamma, m_tgamma,
"gamma(x)\n\nGamma function at x.")
+FUNC1A(lgamma, m_lgamma,
+ "lgamma(x)\n\nNatural logarithm of absolute value of Gamma function at x.")
FUNC1(log1p, log1p, 1,
"log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n"
"The result is computed in a way which is accurate for x near zero.")
@@ -1375,6 +1431,7 @@
{"isinf", math_isinf, METH_O, math_isinf_doc},
{"isnan", math_isnan, METH_O, math_isnan_doc},
{"ldexp", math_ldexp, METH_VARARGS, math_ldexp_doc},
+ {"lgamma", math_lgamma, METH_O, math_lgamma_doc},
{"log", math_log, METH_VARARGS, math_log_doc},
{"log1p", math_log1p, METH_O, math_log1p_doc},
{"log10", math_log10, METH_O, math_log10_doc},