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.
1 2  sumBinomMpfr(n, f, n0 = 0, alternating = TRUE, precBits = 256,
f.k = f(mpfr(k, precBits=precBits)))

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 
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.
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.
Martin Maechler, after conversations with Christophe Dutang.
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.
chooseMpfr
, chooseZ
from package gmp.
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")

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.