Rmpfr (version 0.4-2)

chooseMpfr: Binomial Coefficients and Pochhammer Symbol aka Rising Factorial

Description

Compute binomial coefficients, chooseMpfr(a,n) being mathematically the same as choose(a,n), but using high precision (MPFR) arithmetic.

chooseMpfr.all(n) means the vector choose(n, 1:n), using enough bits for exact computation via MPFR.

pochMpfr() computes the Pochhammer symbol or rising factorial, also called the Pochhammer function, Pochhammer polynomial, ascending factorial, rising sequential product or upper factorial,

$$x^{(n)}=x(x+1)(x+2)\cdots(x+n-1)= \frac{(x+n-1)!}{(x-1)!} = \frac{\Gamma(x+n)}{\Gamma(x)}.$$

Usage

chooseMpfr (a, n)
chooseMpfr.all(n)
pochMpfr(a, n)

Arguments

a

a numeric or mpfr vector.

n

an integer vector; if not of length one, n and a are recycled to the same length.

Value

• an mpfr vector as a.

See Also

choose(n,m) computes the binomial coefficient $C_{n,m}$ which can also be expressed via Pochhammer symbol as $C_{n,m} = (n-m+1)^{(m)}/m!$.

factorialMpfr.

Examples

pochMpfr(100, 4) == 100*101*102*103 # TRUE
a <- 100:110
pochMpfr(a, 10) # exact (but too high precision)
x <- mpfr(a, 70)# should be enough
(px <- pochMpfr(x, 10)) # the same as above (needing only 70 bits)
stopifnot(pochMpfr(a, 10) == px,
          px[1] ==prod(mpfr(100:109, 100)))# used to fail

(c1 <- chooseMpfr(1000:997, 60)) # -> automatic "correct" precision
stopifnot(all.equal(c1, choose(1000:997, 60), tol=1e-12))

## --- Experimenting & Checking
n.set <- c(1:20,50:55, 100:105, 200:203, 500:503, 699:702, 999:1006)
C1 <- C2 <- numeric(length(n.set))
for(i.n in seq_along(n.set)) {
  cat(n <- n.set[i.n],":")
  C1[i.n] <- system.time(c.c <- chooseMpfr.all(n) )[1]
  C2[i.n] <- system.time(c.2 <- chooseMpfr(n, 1:n))[1]
  stopifnot(is.whole(c.c), c.c == c.2,
            if(n > 60) TRUE else all.equal(c.c, choose(n, 1:n), tol = 1e-15))
  cat("[Ok]
")
}
matplot(n.set, cbind(C1,C2), type="b", log="xy",
        xlab = "n", ylab = "system.time(.)  [s]")
legend("topleft", c("chooseMpfr.all(n)", "chooseMpfr(n, 1:n)"),
       pch=as.character(1:2), col=1:2, lty=1:2, bty="n")

## Currently, chooseMpfr.all() is faster only for large n (~= 600)
## That would change if we used C-code for the *.all() version