traverseda
/
cpython
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Minor refactoring in lgamma code, for clarity.

This commit is contained in:
Mark Dickinson 2010-07-07 16:17:31 +00:00
parent 7e4a6ebd42
commit 9c91eb844c
1 changed files with 10 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
Modules/mathmodule.c
View File

@ -69,6 +69,7 @@ extern double copysign(double, double);
static const double pi = 3.141592653589793238462643383279502884197; static const double pi = 3.141592653589793238462643383279502884197;
static const double sqrtpi = 1.772453850905516027298167483341145182798; static const double sqrtpi = 1.772453850905516027298167483341145182798;
static const double logpi = 1.144729885849400174143427351353058711647;
static double static double
sinpi(double x) sinpi(double x)
@ -356,20 +357,15 @@ m_lgamma(double x)
if (absx < 1e-20) if (absx < 1e-20)
return -log(absx); return -log(absx);
/* Lanczos' formula */ /* Lanczos' formula. We could save a fraction of a ulp in accuracy by
if (x > 0.0) { having a second set of numerator coefficients for lanczos_sum that
/* we could save a fraction of a ulp in accuracy by having a absorbed the exp(-lanczos_g) term, and throwing out the lanczos_g
second set of numerator coefficients for lanczos_sum that subtraction below; it's probably not worth it. */
absorbed the exp(-lanczos_g) term, and throwing out the r = log(lanczos_sum(absx)) - lanczos_g;
lanczos_g subtraction below; it's probably not worth it. */ r += (absx - 0.5) * (log(absx + lanczos_g - 0.5) - 1);
r = log(lanczos_sum(x)) - lanczos_g + if (x < 0.0)
(x-0.5)*(log(x+lanczos_g-0.5)-1); /* Use reflection formula to get value for negative x. */
} r = logpi - log(fabs(sinpi(absx))) - log(absx) - r;
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)) if (Py_IS_INFINITY(r))
errno = ERANGE; errno = ERANGE;
return r; return r;