lgammaP11: Log Gamma(p) for Positive 'p' by Pugh's Method (11 Terms)
In DPQ: Density, Probability, Quantile ('DPQ') Computations

lgammaP11

R Documentation

Log Gamma(p) for Positive `p` by Pugh's Method (11 Terms)

Description

Computes \log(\Gamma(p)) for p > 0, where \Gamma(p) = \int_0^\infty s^{p-1} e^{-s}\; ds.

The evaluation hapenss via Pugh's method (approximation with 11 terms), which is a refinement of the Lanczos method (with 6 terms).

This implementation has been used in and for Abergel and Moisan (2020)'s TOMS 1006 algorithm, see the reference and our dltgammaInc(). Given the examples below and example(dltgammaInc), their use of “accurate” does not apply to R's “standard” accuracy.

Usage

lgammaP11(x)

Arguments

`x`	a numeric vector

Value

a numeric vector "as" x (with no attributes, currently).

Author(s)

C source from Remy Abergel and Lionel Moisan, see the reference; original is in ‘1006.zip’ at https://calgo.acm.org/.

Interface to R in DPQ: Martin Maechler

References

Rémy Abergel and Lionel Moisan (2020) Algorithm 1006: Fast and accurate evaluation of a generalized incomplete gamma function, ACM Transactions on Mathematical Software 46(1): 1–24. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1145/3365983")}

Examples

(r <- lgammaP11(1:20))
stopifnot(all.equal(r, lgamma(1:20)))
summary(rE <- sfsmisc::relErrV(r, lgamma(1:20)))
cbind(x=1:20, r, lgamma(1:20))
## at x=1 and 2 it is *not* fully accurate but gives -4.44e-16  (= - 2*EPS )

if(requireNamespace("Rmpfr")) withAutoprint({
    asN  <- Rmpfr::asNumeric
    mpfr <- Rmpfr::mpfr
    x <- seq(1, 20, by = 1/8)
    lgamM <- lgamma(x = mpfr(x, 128)) # uses Rmpfr::Math(.) -> Rmpfr:::.Math.codes
    lgam  <- lgammaP11(x)
    relErrV <- sfsmisc::relErrV
    relE <- asN(relErrV(target = lgamM, current = lgam))
    summary(relE)
    summary(relE.R <- asN(relErrV(lgamM, base::lgamma(x))))
    plot(x, relE, type = "b", col=2, # not very good: at  x ~= 3, goes up to  1e-14 !
         main = "rel.Error(lgammaP11(x)) vs  x")
    abline(h = 0, col = adjustcolor("gray", 0.75), lty = 2)
    lines(x, relE.R, col = adjustcolor(4, 1/2), lwd=2)
    ## ## ==> R's lgamma() is much more accurate
    legend("topright", expression(lgammaP11(x), lgamma(x)),
           col = c(palette()[2], adjustcolor(4, 1/2)), lwd = c(1, 2), pch = c(1, NA), bty = "n")

    ## |absolute rel.Err| on log-scale:
    cEps <- 2^-52 ; lc <- adjustcolor("gray", 0.75)
    plot(x, abs(relE), type = "b", col=2, log = "y", ylim = cEps * c(2^-4, 2^6), yaxt = "n",
         main = "|rel.Error(lgammaP11(x))| vs  x -- log-scale")
    sfsmisc::eaxis(2, cex.axis = 3/4, col.axis = adjustcolor(1, 3/4))
    abline(h= cEps, col = lc, lty = 2); axis(4, at=cEps, quote(epsilon[C]), las=1, col.axis=lc)
    lines(x, pmax(2^-6*cEps, abs(relE.R)), col = adjustcolor(4, 1/2), lwd=2)
    ## ## ==> R's lgamma() is *much* more accurate
    legend("topright", expression(lgammaP11(x), lgamma(x)), bg = adjustcolor("gray", 1/4),
           col = c(palette()[2], adjustcolor(4, 1/2)), lwd = c(1, 2), pch = c(1, NA), bty = "n")
    mtext(sfsmisc::shortRversion(), cex = .75, adj = 1)
})

DPQ documentation built on July 18, 2025, 3:01 p.m.