gam1d: Compute 1/Gamma(x+1) - 1 Accurately
In DPQ: Density, Probability, Quantile ('DPQ') Computations

gam1d

R Documentation

Compute 1/Gamma(x+1) - 1 Accurately

Description

Computes 1/\Gamma(a+1) - 1 accurately in [-0.5, 1.5] for numeric argument a; For "mpfr" numbers, the precision is increased intermediately such that a+1 should not lose precision.

FIXME: "Pure-R" implementation is in ‘ ~/R/Pkgs/DPQ/TODO_R_versions_gam1_etc.R ’

Usage


gam1d(a, warnIf = TRUE, verbose = FALSE)

Arguments

`a`	a numeric or numeric-alike, typically inheriting from `class` `"mpfr"`.
`warnIf`	logical if a `warning` should be signalled when `a` is not in the “proper” range `[-0.5, 1.5]`.
`verbose`	logical indicating if some output from C code execution should be printed to the console.

Details

https://dlmf.nist.gov/ states the well-know Taylor series for

\frac{1}{\Gamma(z)} = \sum_{k=1}^\infty c_k z^k

with c_1 = 1, c_2 = \gamma, (Euler's gamma, \gamma = 0.5772..., with recursion c_k = (\gamma c_{k-1} - \zeta(2) c_{k-2} ... +(-1)^k \zeta(k-1) c_1) /(k-1).

Hence,

\frac{1}{\Gamma(z+1)} = z+1 + \sum_{k=2}^\infty c_k (z+1)^k

\frac{1}{\Gamma(z+1)} -1 = z + \gamma*(z+1)^2 + \sum_{k=3}^\infty c_k (z+1)^k

Consequently, for \zeta_k := \zeta(k), c_3 = (\gamma^2 - \zeta_2)/2, c_4 = \gamma^3/6 - \gamma \zeta_2/2 + \zeta_3/3.

  gam <- Const("gamma", 128)
  z <- Rmpfr::zeta(mpfr(1:7, 128))
  (c3 <- (gam^2 -z[2])/2)                       # -0.655878071520253881077019515145
  (c4 <- (gam*c3 - z[2]*c2 + z[3])/3)           # -0.04200263503409523552900393488
  (c4 <- gam*(gam^2/6 - z[2]/2) + z[3]/3)
  (c5 <- (gam*c4 - z[2]*c3 + z[3]*c2 - z[4])/4) # 0.1665386113822914895017007951
  (c5 <- (gam^4/6 - gam^2*z[2] + z[2]^2/2 + gam*z[3]*4/3 - z[4])/4)

Value

a numeric-alike vector like a.

Author(s)

Martin Maechler building on C code of TOMS 708

References

TOMS 708, see pbeta

Examples

g1 <- function(u) 1/gamma(u+1) - 1
u <- seq(-.5, 1.5, by=1/16); set.seed(1); u <- sample(u) # permuted (to check logic)


g11   <- vapply(u, gam1d, 1)
gam1d. <- gam1d(u)
stopifnot( all.equal(g1(u), g11) )
stopifnot( identical(g11, gam1d.) )

## Comparison of g1() and gam1d(), slightly extending the [-.5, 1.5] interval:
u <- seq(-0.525, 1.525, length.out = 2001)
mg1 <- cbind(g1 = g1(u), gam1d = gam1d(u))
clr <- adjustcolor(1:2, 1/2)
matplot(u, mg1, type = "l", lty = 1, lwd=1:2, col=clr) # *no* visual difference
## now look at *relative* errors
relErrV <- sfsmisc::relErrV
relE <- relErrV(mg1[,"gam1d"], mg1[,"g1"])
plot(u, relE, type = "l")
plot(u, abs(relE), type = "l", log = "y",
     main = "|rel.diff|  gam1d() vs 'direct' 1/gamma(u+1) - 1")

## now {Rmpfr} for "truth" :
if(requireNamespace("Rmpfr")) withAutoprint({
    asN  <- Rmpfr::asNumeric; mpfr <- Rmpfr::mpfr
    gam1M <- g1(mpfr(u, 512)) # "cheap": high precision avoiding "all" cancellation
    relE <- asN(relErrV(gam1M, gam1d(u)))
    plot(relE ~ u, type="l", ylim = c(-1,1) * 2.5e-15,
         main = expression("Relative Error of " ~~ gam1d(u) %~~% frac(1, Gamma(u+1)) - 1))
    grid(lty = 3); abline(v = c(-.5, 1.5), col = adjustcolor(4, 1/2), lty=2, lwd=2)
})





if(requireNamespace("Rmpfr") && FALSE) { 
  
## Comparison using Rmpfr; slightly extending the [-.5, 1.5] interval:
##	{relErrV(), mpfr(), asN() defined above}

u <- seq(-0.525, 1.525, length.out = 2001)
gam1M <- gam1(mpfr(u, 128))
relE <- asN(relErrV(gam1M, gam1d(u)))

plot(relE ~ u, type="l", ylim = c(-1,1) * 2.5e-15,
     main = expression("Relative Error of " ~~ gam1d(u) == frac(1, Gamma(u+1)) - 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 <- asN(relErrV(gam1M, g1(u)))

plot(relE ~ u, type="l", ylim = c(-1,1) * 1e-14,
     main = expression("Relative Error of " ~~ gam1d(u) == frac(1, Gamma(u+1)) - 1))
lines(relED ~ u, col = adjustcolor(2, 1/2), lwd = 2)
# mtext("comparing with direct formula   1/gamma(u+1) - 1")
legend("top", c("gam1d(u)", "1/gamma(u+1) - 1"), col = 1:2, lwd=1:2, bty="n")
## direct is clearly *worse* , but not catastrophical
}

DPQ documentation built on Nov. 3, 2024, 3 a.m.