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).*

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

`a` |
a numeric or |

`n` |
an integer vector; if not of length one, |

`rnd.mode` |
a 1-letter string specifying how |

`precBits` |
integer or NULL for increasing the default precision of the result. |

`k0` |
integer scalar |

`alternating` |
logical, for |

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`

.

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.

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

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.