Rmpfr (version 0.5-3)

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);$$ 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)))

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

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.

concept

• Rice integral
• Forward Difference

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

## "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
stopifnot(
    all.equal(resMpfr[1:10], res.R[1:10], tol=1e-12),
    all.equal(resMpfr[1:20], res.R[1:20], tol=1e-9 ),
    all.equal(resMpfr[1:30], res.R[1:30], tol=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")