chooseMpfr: Binomial Coefficients and Pochhammer Symbol aka Rising...

In Rmpfr: Interface R to MPFR - Multiple Precision Floating-Point Reliable

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. 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)\cdots(x+n-1)= \frac{(x+n-1)!}{(x-1)!} = \frac{\Gamma(x+n)}{\Gamma(x)}.

Usage

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\cdot C_{n,m} for m \in k0:n.

Note

Currently this works via a (C level) for(i in 1:n)-loop which really slow for large n, say 10^6, with computational cost O(n^2). In such cases, if you need high precision choose(a,n) (or Pochhammer(a,n)) for large n, preferably 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

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

## If you want to measure more:
measureMore <- TRUE
measureMore <- FALSE
if(measureMore) { ## takes ~ 2 minutes (on "lynne", Intel i7-7700T, ~2019)
  n.s <- 2^(5:20)
  r <- lapply(n.s, function(n) {
      N <- ceiling(10000/n)
      cat(sprintf("n=%9g => N=%d: ",n,N))
      ct <- system.time(C <- replicate(N, chooseMpfr(n, n/2)))
      cat("[Ok]\n")
      list(C=C, ct=ct/N)
  })
  print(ct.n <- t(sapply(r, `[[`, "ct")))
  hasSfS <- requireNamespace("sfsmisc")
  plot(ct.n[,"user.self"] ~ n.s, xlab=quote(n), ylab="system.time(.) [s]",
       main = "CPU Time for  chooseMpfr(n, n/2)",
       log ="xy", type = "b", axes = !hasSfS)
  if(hasSfS) for(side in 1:2) sfsmisc::eaxis(side)
  summary(fm <- lm(log(ct.n[,"user.self"]) ~ log(n.s), subset = n.s >= 10^4))
  ## --> slope ~= 2  ==>  It's  O(n^2)
  nn <- 2^seq(11,21, by=1/16) ; Lcol <- adjustcolor(2, 1/2)
  bet <- coef(fm)
  lines(nn, exp(predict(fm, list(n.s = nn))), col=Lcol, lwd=3)
  text(500000,1, substitute(AA %*% n^EE,
                            list(AA = signif(exp(bet[1]),3),
                                 EE = signif(    bet[2], 3))), col=2)
} # measure more

Rmpfr documentation built on Nov. 18, 2024, 3:01 p.m.