Mercurial > cpython-lru_cache-weakref

comparison Modules/mathmodule.c @ 57655:cb22eaee86c9 legacy-trunk

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Issue #3366: Add lgamma function to math module.
author Mark Dickinson <dickinsm@gmail.com>
date Fri, 11 Dec 2009 17:29:33 +0000
parents 7e02371aa506
children b206ad8cf5bf
comparison
equal deleted inserted replaced
57650:ae9396260835 57655:cb22eaee86c9
318 } 318 }
319 if (Py_IS_INFINITY(r)) 319 if (Py_IS_INFINITY(r))
320 errno = ERANGE; 320 errno = ERANGE;
321 return r; 321 return r;
322 } 322 }
323
324 /*
325 lgamma: natural log of the absolute value of the Gamma function.
326 For large arguments, Lanczos' formula works extremely well here.
327 */
328
329 static double
330 m_lgamma(double x)
331 {
332 double r, absx;
333
334 /* special cases */
335 if (!Py_IS_FINITE(x)) {
336 if (Py_IS_NAN(x))
337 return x; /* lgamma(nan) = nan */
338 else
339 return Py_HUGE_VAL; /* lgamma(+-inf) = +inf */
340 }
341
342 /* integer arguments */
343 if (x == floor(x) && x <= 2.0) {
344 if (x <= 0.0) {
345 errno = EDOM; /* lgamma(n) = inf, divide-by-zero for */
346 return Py_HUGE_VAL; /* integers n <= 0 */
347 }
348 else {
349 return 0.0; /* lgamma(1) = lgamma(2) = 0.0 */
350 }
351 }
352
353 absx = fabs(x);
354 /* tiny arguments: lgamma(x) ~ -log(fabs(x)) for small x */
355 if (absx < 1e-20)
356 return -log(absx);
357
358 /* Lanczos' formula */
359 if (x > 0.0) {
360 /* we could save a fraction of a ulp in accuracy by having a
361 second set of numerator coefficients for lanczos_sum that
362 absorbed the exp(-lanczos_g) term, and throwing out the
363 lanczos_g subtraction below; it's probably not worth it. */
364 r = log(lanczos_sum(x)) - lanczos_g +
365 (x-0.5)*(log(x+lanczos_g-0.5)-1);
366 }
367 else {
368 r = log(pi) - log(fabs(sinpi(absx))) - log(absx) -
369 (log(lanczos_sum(absx)) - lanczos_g +
370 (absx-0.5)*(log(absx+lanczos_g-0.5)-1));
371 }
372 if (Py_IS_INFINITY(r))
373 errno = ERANGE;
374 return r;
375 }
376
323 377
324 /* 378 /*
325 wrapper for atan2 that deals directly with special cases before 379 wrapper for atan2 that deals directly with special cases before
326 delegating to the platform libm for the remaining cases. This 380 delegating to the platform libm for the remaining cases. This
327 is necessary to get consistent behaviour across platforms. 381 is necessary to get consistent behaviour across platforms.
637 FUNC1(floor, floor, 0, 691 FUNC1(floor, floor, 0,
638 "floor(x)\n\nReturn the floor of x as a float.\n" 692 "floor(x)\n\nReturn the floor of x as a float.\n"
639 "This is the largest integral value <= x.") 693 "This is the largest integral value <= x.")
640 FUNC1A(gamma, m_tgamma, 694 FUNC1A(gamma, m_tgamma,
641 "gamma(x)\n\nGamma function at x.") 695 "gamma(x)\n\nGamma function at x.")
696 FUNC1A(lgamma, m_lgamma,
697 "lgamma(x)\n\nNatural logarithm of absolute value of Gamma function at x.")
642 FUNC1(log1p, log1p, 1, 698 FUNC1(log1p, log1p, 1,
643 "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n" 699 "log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n"
644 "The result is computed in a way which is accurate for x near zero.") 700 "The result is computed in a way which is accurate for x near zero.")
645 FUNC1(sin, sin, 0, 701 FUNC1(sin, sin, 0,
646 "sin(x)\n\nReturn the sine of x (measured in radians).") 702 "sin(x)\n\nReturn the sine of x (measured in radians).")
1373 {"gamma", math_gamma, METH_O, math_gamma_doc}, 1429 {"gamma", math_gamma, METH_O, math_gamma_doc},
1374 {"hypot", math_hypot, METH_VARARGS, math_hypot_doc}, 1430 {"hypot", math_hypot, METH_VARARGS, math_hypot_doc},
1375 {"isinf", math_isinf, METH_O, math_isinf_doc}, 1431 {"isinf", math_isinf, METH_O, math_isinf_doc},
1376 {"isnan", math_isnan, METH_O, math_isnan_doc}, 1432 {"isnan", math_isnan, METH_O, math_isnan_doc},
1377 {"ldexp", math_ldexp, METH_VARARGS, math_ldexp_doc}, 1433 {"ldexp", math_ldexp, METH_VARARGS, math_ldexp_doc},
1434 {"lgamma", math_lgamma, METH_O, math_lgamma_doc},
1378 {"log", math_log, METH_VARARGS, math_log_doc}, 1435 {"log", math_log, METH_VARARGS, math_log_doc},
1379 {"log1p", math_log1p, METH_O, math_log1p_doc}, 1436 {"log1p", math_log1p, METH_O, math_log1p_doc},
1380 {"log10", math_log10, METH_O, math_log10_doc}, 1437 {"log10", math_log10, METH_O, math_log10_doc},
1381 {"modf", math_modf, METH_O, math_modf_doc}, 1438 {"modf", math_modf, METH_O, math_modf_doc},
1382 {"pow", math_pow, METH_VARARGS, math_pow_doc}, 1439 {"pow", math_pow, METH_VARARGS, math_pow_doc},