igamma: Incomplete Gamma Function

Description

For MPFR version >= 3.2.0, the following MPFR library function is provided: mpfr_gamma_inc(a,x), the R interface of which is igamma(a,x), where igamma(a,x) is the “upper” incomplete gamma function $$\Gamma(a,x) :=: \Gamma(a) - \gamma(a,x),$$ where $$\gamma(a,x) := \int_0^x t^{a-1} e^{-t} dt,$$ and hence $$\Gamma(a,x) := \int_x^\infty t^{a-1} e^{-t} dt,$$ and $$\Gamma(a) := \gamma(a, \infty).$$

As R's pgamma(x,a) is $$\code{pgamma(x, a)} := \gamma(a,x) / \Gamma(a),$$ we get

        igamma(a,x) ==  gamma(a) * pgamma(x, a, lower.tail=FALSE)

Usage

igamma(a, x, rnd.mode = c("N", "D", "U", "Z", "A"))

Arguments

a, x

an object of class mpfr or numeric.

rnd.mode

a 1-letter string specifying how rounding should happen at C-level conversion to MPFR, see mpfr.

Value

a numeric vector of “common length”, recyling along a and x.

References

NIST Digital Library of Mathematical Functions, section 8.2. http://dlmf.nist.gov/8.2.i

Wikipedia (2019). Incomplete gamma function; https://en.wikipedia.org/wiki/Incomplete_gamma_function

Examples

Run this code

# NOT RUN {
## show how close pgamma() is :
x <- c(seq(0,20, by=1/4), 21:50, seq(55, 100, by=5))
if(mpfrVersion() >= "3.2.0") { print(
all.equal(igamma(Const("pi", 80), x),
          pgamma(x, pi, lower.tail=FALSE) * gamma(pi),
          tol=0, formatFUN = function(., ...) format(., digits = 7)) #-> 3.13e-16
)
## and ensure *some* closeness:
stopifnot(exprs = {
   all.equal(igamma(Const("pi", 80), x),
             pgamma(x, pi, lower.tail=FALSE) * gamma(pi),
             tol = 1e-15)
})
} # only if MPFR version >= 3.2.0
# }