lgamma1pM: Compute log( Gamma(x+1) ) Arbitrarily (MPFR) Accurately

In DPQmpfr: DPQ (Density, Probability, Quantile) Distribution Computations using MPFR

Compute log( Gamma(x+1) ) Arbitrarily (MPFR) Accurately

Description

Computes \log \Gamma(x+1) accurately notably when |x| \ll 1. For "mpfr" numbers, the precision is increased intermediately such that a+1 should not lose precision.

R's "own" double prec version is soon available in package in DPQ, under the name gamln1() (from TOMS 708).

Usage

lgamma1pM(a, usePr = NULL, DPQmethod = c("lgamma1p", "algam1"))

Arguments

a	a numeric or numeric-alike vector, typically inheriting from class "mpfr".
usePr	positive integer specifying the precision in bits, or NULL when a smart default will be used.
DPQmethod	a character string; must be the name of an lgamma1p()-alike function from package DPQ. It will be called in case of is.numeric(a) (and when DPQ is available).

Value

a numeric-alike vector like a.

Author(s)

Martin Maechler

References

TOMS 708, see pbeta

See Also

lgamma() (and gamma() (same page)), and our algdivM(); further, package DPQ's lgamma1p() and (if already available) gamln1().

Examples

## Package {DPQ}'s  lgamma1p():
lgamma1p <- DPQ::lgamma1p
lg1 <- function(u) lgamma(u+1) # the simple direct form
u <- seq(-.5, 1.5, by=1/16); set.seed(1); u <- sample(u) # permuted (to check logic)

g11   <- vapply(u, lgamma1p, numeric(1))
lgamma1p. <- lgamma1p(u)
all.equal(lg1(u), g11, tolerance = 0) # see 3.148e-16
stopifnot(exprs = {
    all.equal(lg1(u), g11, tolerance = 2e-15)
    identical(g11, lgamma1p.)
})

## Comparison using Rmpfr; slightly extending the [-.5, 1.5] interval:
u <- seq(-0.525, 1.525, length.out = 2001)
lg1p  <- lgamma1pM(   u)
lg1pM <- lgamma1pM(Rmpfr::mpfr(u, 128))
asNumeric <- Rmpfr::asNumeric
relErrV   <- sfsmisc::relErrV
if(FALSE) { # DPQ "latest" version __FIXME__
lng1  <- DPQ::lngam1(u)
relE <- asNumeric(relErrV(lg1pM, cbind(lgamma1p = lg1p, lngam1 = lng1)))
} else {
relE <- asNumeric(relErrV(lg1pM, cbind(lgamma1p = lg1p)))#, lngam1 = lng1)))
}

## FIXME: lgamma1p() is *NOT* good around u =1. -- even though it should
##        and the R-only vs (not installed) *does* "work" (is accurate there) ?????
## --> ~/R/Pkgs/DPQ/TODO_R_versions_gam1_etc.R
if(FALSE) {
matplot(u, relE, type="l", ylim = c(-1,1) * 2.5e-15,
     main = expression("relative error of " ~~ lgamma1p(u) == log( Gamma(u+1) )))
} else {
plot(relE ~ u, type="l", ylim = c(-1,1) * 2.5e-15,
     main = expression("relative error of " ~~ lgamma1p(u) == log( Gamma(u+1) )))
}
grid(lty = 3); abline(v = c(-.5, 1.5), col = adjustcolor(4, 1/2), lty=2, lwd=2)

## what about the direct formula -- how bad is it really ?
relED <- asNumeric(relErrV(lg1pM, lg1(u)))
lines(relED ~ u, col = adjustcolor(2, 1/2), lwd = 2)

DPQmpfr documentation built on Sept. 11, 2024, 8:54 p.m.