chooseMpfr: Binomial Coefficients and Pochhammer Symbol aka Rising... In Rmpfr: R MPFR - Multiple Precision Floating-Point Reliable

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.

`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 ```

Rmpfr documentation built on July 31, 2018, 1 a.m.