Binomial Coefficients and Pochhammer Symbol aka Rising Factorial

Share:

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. However, chooseMpfr.all() is now deprecated in favor of chooseZ from package gmp, as that is now vectorized.

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)...(x+n-1) = (x+n-1)! / (x-1)! = Gamma(x+n) / Gamma(x).

Usage

1
2
3
chooseMpfr (a, n, rnd.mode = c("N","D","U","Z","A"))
chooseMpfr.all(n, precBits=NULL, k0=1, alternating=FALSE)
pochMpfr(a, n, rnd.mode = c("N","D","U","Z","A"))

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.

rnd.mode

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

precBits

integer or NULL for increasing the default precision of the result.

k0

integer scalar

alternating

logical, for chooseMpfr.all(), indicating if alternating sign coefficients should be returned, see below.

Value

For

chooseMpfr(), pochMpfr():

an mpfr vector of length max(length(a), length(n));

chooseMpfr.all(n, k0):

a mpfr vector of length n-k0+1, of binomial coefficients C[n,m] or, if alternating is true, (-1)^m * C[n,m] for m in k0:n.

Note

If you need high precision choose(a,n) (or Pochhammer(a,n)) for large n, maybe better work with the corresponding factorial(mpfr(..)), or gamma(mpfr(..)) terms.

See Also

choose(n,m) (base R) computes the binomial coefficient C[n,m] which can also be expressed via Pochhammer symbol as C[n,m] = (n-m+1)^(m) / m!.

chooseZ from package gmp; for now, factorialMpfr.

For (alternating) binomial sums, directly use sumBinomMpfr, as that is potentially more efficient.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
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), tolerance=1e-12))

## --- Experimenting & Checking
n.set <- c(1:10, 20, 50:55, 100:105, 200:203, 300:303, 500:503,
           699:702, 999:1001)
if(!Rmpfr:::doExtras()) { ## speed up: smaller set
  n. <- n.set[-(1:10)]
  n.set <- c(1:10, n.[ c(TRUE, diff(n.) > 1)])
}
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), tolerance = 1e-15))
  cat(" [Ok]\n")
}
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 (>= 300)
## That would change if we used C-code for the *.all() version

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.