lgamma1p: Accurate `log(gamma(a+1))`

Description

Compute $l Γ_{1} (a) := \log Γ (a + 1) = \log (a \cdot Γ (a)) = \log a + \log Γ (a),$ which is “in principle” the same as log(gamma(a+1)) or lgamma(a+1), accurately also for (very) small $a$ $(0 < a < 0.5)$ .

Usage

lgamma1p (a, tol_logcf = 1e-14, f.tol = 1, ...)	
lgamma1p.(a, cutoff.a = 1e-6, k = 3)		
lgamma1p_series(x, k)               		
lgamma1pC(x)

Value

a numeric vector with the same attributes as a.

Arguments

a, x: a numeric vector.
tol_logcf: for lgamma1p(): a non-negative number passed to logcf() (and log1pmx() which calls logcf()).
f.tol: numeric (factor) used in log1pmx(*, tol_logcf = f.tol * tol_logcf).
...: further optional arguments passed on to log1pmx().
cutoff.a: for lgamma1p.(): a positive number indicating the cutoff to switch from ...

k: an integer, the number of terms in the series expansion used internally.

Author

Morten Welinder (C code of Jan 2005, see R's bug issue tools:::Rd_expr_PR(7307)) for lgamma1p().

Martin Maechler, notably for lgamma1p_series() which works with package Rmpfr but otherwise may be much less accurate than Morten's 40 term series!

Details

lgamma1p() is an R translation of the function (in Fortran) in Didonato and Morris (1992) which uses a 40-degree polynomial approximation.

lgamma1p_series(x, k) is Taylor series approximation of order k, (derived via Maple), which is $- γ x + π^{2} x^{2} / 12 + O (x^{3})$ , where $γ$ is Euler's constant 0.5772156649. ...

lgamma1pC() is an interface to R C API (Mathlib / Rmath.h) function.

References

Didonato, A. and Morris, A., Jr, (1992) Algorithm 708: Significant digit computation of the incomplete beta function ratios. ACM Transactions on Mathematical Software, 18, 360--373; see also pbeta.

Examples

Run this code

curve(-log(x*gamma(x)), 1e-30, .8, log="xy", col="gray50", lwd = 3,
      axes = FALSE, ylim = c(1e-30,1))
sfsmisc::eaxis(1); sfsmisc::eaxis(2)
at <- 10^(1-4*(0:8))
abline(h = at, v = at, col = "lightgray", lty = "dotted")
curve(-lgamma( 1+x), add=TRUE, col="red2", lwd=1/2)# underflows even earlier
curve(-lgamma1p (x), add=TRUE, col="blue") -> lgxy
curve(-lgamma1p.(x), add=TRUE, col=adjustcolor("forest green",1/4),
      lwd = 5, lty = 2)
for(k in 1:7)
  curve(-lgamma1p_series(x, k=k), add=TRUE, col=paste0("gray",30+k*8), lty = 3)
stopifnot(with(lgxy, all.equal(y, -lgamma1pC(x))))

Run the code above in your browser using DataLab

Get 50% off unlimited learning