(Alternating) Binomial Sums via Rmpfr
Description
Compute (alternating) binomial sums via highprecision arithmetic.
If sBn(f, n) :=sumBinomMpfr(n, f)
, (default
alternating
is true, and n0 = 0
),
sBn(f,n) = sum(k = n0:n ; (1)^(nk) choose(n,k) * f(k)) = Δ^n f,
see Details for the nth forward difference operator
Δ^n f. If alternating
is false, the
(1)^(nk) factor is dropped (or replaced by 1) above.
Such sums appear in different contexts and are typically challenging, i.e., currently impossible, to evaluate reliably as soon as n is larger than around 5070.
Usage
1 2  sumBinomMpfr(n, f, n0 = 0, alternating = TRUE, precBits = 256,
f.k = f(mpfr(k, precBits=precBits)))

Arguments
n 
upper summation index (integer). 
f 

n0 
lower summation index, typically 
alternating 
logical indicating if the sum is alternating, see below. 
precBits 
the number of bits for MPFR precision, see

f.k 
can be specified instead of 
Details
The alternating binomial sum sB(f,n) := sumBinom(n, f, n0=0) is equal to the nth forward difference operator Δ^n f,
sB(f,n) = Δ^n f,
where
Delta^n f = sum(k = n0:n ; (1)^(nk) choose(n,k) * f(k)),
is the nfold iterated forward difference Δ f(x) = f(x+1)  f(x) (for x = 0).
The current implementation might be improved in the future, notably
for the case where
sB(f,n)=sumBinomMpfr(n, f, *)
is to be computed for a whole sequence
n = 1,...,N.
Value
an mpfr
number of precision precBits
.
s. If alternating
is true (as per default),
s = sum(k = n0:n ; (1)^k choose(n,k) * f(k)),
if alternating
is false, the (1)^k factor is dropped (or
replaced by 1) above.
Author(s)
Martin Maechler, after conversations with Christophe Dutang.
References
Wikipedia (2012) The N\"orlundRice integral, http://en.wikipedia.org/wiki/Rice_integral
Flajolet, P. and Sedgewick, R. (1995) Mellin Transforms and Asymptotics: Finite Differences and Rice's Integrals, Theoretical Computer Science 144, 101–124.
See Also
chooseMpfr
, chooseZ
from package gmp.
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  ## "naive" R implementation:
sumBinom < function(n, f, n0=0, ...) {
k < n0:n
sum( choose(n, k) * (1)^(nk) * f(k, ...))
}
## compute sumBinomMpfr(.) for a whole set of 'n' values:
sumBin.all < function(n, f, n0=0, precBits = 256, ...)
{
N < length(n)
precBits < rep(precBits, length = N)
ll < lapply(seq_len(N), function(i)
sumBinomMpfr(n[i], f, n0=n0, precBits=precBits[i], ...))
sapply(ll, as, "double")
}
sumBin.all.R < function(n, f, n0=0, ...)
sapply(n, sumBinom, f=f, n0=n0, ...)
n.set < 5:80
system.time(res.R < sumBin.all.R(n.set, f = sqrt)) ## instantaneous..
system.time(resMpfr < sumBin.all (n.set, f = sqrt)) ## ~ 0.6 seconds
matplot(n.set, cbind(res.R, resMpfr), type = "l", lty=1,
ylim = extendrange(resMpfr, f = 0.25), xlab = "n",
main = "sumBinomMpfr(n, f = sqrt) vs. R double precision")
legend("topleft", leg=c("double prec.", "mpfr"), lty=1, col=1:2, bty = "n")
