DPQmpfr (version 0.3-3)

lgamma1pM: 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"))

Value

a numeric-alike vector like a.

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).

Author

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)