DPQ (version 0.5-9)

gam1d: 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)

Value

a numeric-alike vector like a.

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.

Author

Martin Maechler building on C code of TOMS 708

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)

References

TOMS 708, see pbeta

See Also

gamma.

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
}