

   SSppeecciiaall FFuunnccttiioonnss ooff MMaatthheemmaattiiccss

        beta(a, b)
        lbeta(a, b)
        gamma(x)
        lgamma(x)
        digamma(x)
        trigamma(x)
        tetragamma(x)
        pentagamma(x)
        choose(n,k)
        lchoose(n,k)

   DDeessccrriippttiioonn::

        The functions `beta' and `lbeta' return the beta func-
        tion and the natural logarithm of the beta function,

               B(a,b) = (Gamma(a)Gamma(b))/(Gamma(a+b)).

        The functions `gamma' and `lgamma' return the gamma
        function Gamma(x) and the natural logarithm of the
        absolute value of the gamma function.

        The functions `digamma', `trigamma', `tetragamma' and
        `pentagamma' return the first, second, third and fourth
        derivatives of the logarithm of the gamma function.

        `digamma(x)' = psi(x) = d/dx {ln Gamma(x)} = Gamma'(x) / Gamma(x)

        The functions `choose' and `lchoose' return binomial
        coefficients and their logarithms.

   RReeffeerreenncceess::

        Abramowitz, M. and Stegun, I. A. (1972).
        Handbook of Mathematical Functions, Dover, New York;
        Chapter 6: Gamma and Related Functions.

   SSeeee AAllssoo::

        `Arithmetic' for simple and `Math' for miscellaneous
        mathematical functions.

   EExxaammpplleess::

        choose(5, 2)
        for (n in 0:10) print(choose(n, k = 0:n))

        curve(gamma(x),-3,4, n=1001, ylim=c(-10,100),
              col="red", lwd=2, main="gamma(x)")
        abline(h=0,v=0, lty=3, col="midnightblue")

        x <- seq(.1, 4, length = 201); dx <- diff(x)[1]
        par(mfrow = c(2, 3))
        for (ch in c("", "l","di","tri","tetra","penta")) {
          is.deriv <- nchar(ch) >= 2
          if (is.deriv) dy <- diff(y) / dx
          nm <- paste(ch, "gamma", sep = "")
          y <- get(nm)(x)
          plot(x, y, type = "l", main = nm, col = "red")
          abline(h = 0, col = "lightgray")
          if (is.deriv) lines(x[-1], dy, col = "blue", lty = 2)
        }
        par(mfrow = c(2, 2))

