Rmpfr (version 0.9-5)

sumBinomMpfr: (Alternating) Binomial Sums via Rmpfr

Description

Compute (alternating) binomial sums via high-precision arithmetic. If \(sBn(f, n) :=\)sumBinomMpfr(n, f), (default alternating is true, and n0 = 0), $$sBn(f,n) = \sum_{k = n0}^n (-1)^(n-k) {n \choose k}\cdot f(k) = \Delta^n f,$$ see Details for the \(n\)-th forward difference operator \(\Delta^n f\). If alternating is false, the \((-1)^(n-k)\) 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 \(50--70\).

Usage

sumBinomMpfr(n, f, n0 = 0, alternating = TRUE, precBits = 256,
             f.k = f(mpfr(k, precBits=precBits)))

Value

an mpfr number of precision precBits.

\(s\). If alternating is true (as per default),

$$s = \sum_{k = n0}^n (-1)^k {n \choose k}\cdot f(k),$$

if alternating is false, the \((-1)^k\) factor is dropped (or replaced by \(1\)) above.

Arguments

n

upper summation index (integer).

f

function to be evaluated at \(k\) for k in n0:n (and which must return one value per k).

n0

lower summation index, typically 0 (= default) or 1.

alternating

logical indicating if the sum is alternating, see below.

precBits

the number of bits for MPFR precision, see mpfr.

f.k

can be specified instead of f and precBits, and must contain the equivalent of its default, f(mpfr(k, precBits=precBits)).

Author

Martin Maechler, after conversations with Christophe Dutang.

Details

The alternating binomial sum \(sB(f,n) := sumBinom(n, f, n0=0)\) is equal to the \(n\)-th forward difference operator \(\Delta^n f\), $$sB(f,n) = \Delta^n f,$$ where $$\Delta^n f = \sum_{k=0}^{n} (-1)^{n-k}{n \choose k}\cdot f(k),$$ is the \(n\)-fold iterated forward difference \(\Delta 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,\dots,N\).

References

Wikipedia (2012) The N\"orlund-Rice integral, https://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

## "naive" R implementation:
sumBinom <- function(n, f, n0=0, ...) {
  k <- n0:n
  sum( choose(n, k) * (-1)^(n-k) * 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
# \dontshow{
stopifnot(
    all.equal(resMpfr[1:10], res.R[1:10], tolerance=1e-12),
    all.equal(resMpfr[1:20], res.R[1:20], tolerance=1e-9 ),
    all.equal(resMpfr[1:30], res.R[1:30], tolerance=1e-6 ))
# }
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")