R: Binomial Coefficients and Pochhammer Symbol aka Rising...
chooseMpfr
R Documentation
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”,
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 ink0: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) (baseR) 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
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(Rmpfr)
Loading required package: gmp
Attaching package: 'gmp'
The following objects are masked from 'package:base':
%*%, apply, crossprod, matrix, tcrossprod
C code of R package 'Rmpfr': GMP using 64 bits per limb
Attaching package: 'Rmpfr'
The following objects are masked from 'package:stats':
dbinom, dnorm, dpois, pnorm
The following objects are masked from 'package:base':
cbind, pmax, pmin, rbind
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/Rmpfr/chooseMpfr.Rd_%03d_medium.png", width=480, height=480)
> ### Name: chooseMpfr
> ### Title: Binomial Coefficients and Pochhammer Symbol aka Rising Factorial
> ### Aliases: chooseMpfr chooseMpfr.all pochMpfr
> ### Keywords: arith
>
> ### ** Examples
>
> pochMpfr(100, 4) == 100*101*102*103 # TRUE
[1] TRUE
> a <- 100:110
> pochMpfr(a, 10) # exact (but too high precision)
11 'mpfr' numbers of precision 320 bits
[1] 154711039801002048000 170182143781102252800 187031860987151980800
[4] 205368317946676684800 225306989591985100800 246971123206599052800
[7] 270492182559608486400 296010312989760230400 323674828222448102400
[10] 353644719724526630400 386089189424024486400
> x <- mpfr(a, 70)# should be enough
> (px <- pochMpfr(x, 10)) # the same as above (needing only 70 bits)
11 'mpfr' numbers of precision 70 bits
[1] 154711039801002048000 170182143781102252800 187031860987151980800
[4] 205368317946676684800 225306989591985100800 246971123206599052800
[7] 270492182559608486400 296010312989760230400 323674828222448102400
[10] 353644719724526630400 386089189424024486400
> stopifnot(pochMpfr(a, 10) == px,
+ px[1] ==prod(mpfr(100:109, 100)))# used to fail
>
> (c1 <- chooseMpfr(1000:997, 60)) # -> automatic "correct" precision
4 'mpfr' numbers of precision 384 bits
[1] 19742748621859838806445290867590842097270339314978449118678002674641952503075171748033908989976400
[2] 18558183704548248478058573415535391571434118956079742171557322514163435352890661443151874450577816
[3] 17443578076647452773670671108296028714290928628387265164256582423222688484849180275395005114206776
[4] 16394865967830972646997083666915505945896684422271798320714102518018919638064660419158832462050056
> 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")
+ }
1 : [Ok]
2 : [Ok]
3 : [Ok]
4 : [Ok]
5 : [Ok]
6 : [Ok]
7 : [Ok]
8 : [Ok]
9 : [Ok]
10 : [Ok]
20 : [Ok]
50 : [Ok]
100 : [Ok]
200 : [Ok]
300 : [Ok]
500 : [Ok]
699 : [Ok]
999 : [Ok]
> matplot(n.set, cbind(C1,C2), type="b", log="xy",
+ xlab = "n", ylab = "system.time(.) [s]")
Warning messages:
1: In xy.coords(x, y, xlabel, ylabel, log = log) :
12 y values <= 0 omitted from logarithmic plot
2: In xy.coords(x, y, xlabel, ylabel, log) :
1 y value <= 0 omitted from logarithmic plot
> 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
>
>
>
>
>
> dev.off()
null device
1
>