sumBinomMpfr: (Alternating) Binomial Sums via Rmpfr

In Rmpfr: Interface R to MPFR - Multiple Precision Floating-Point Reliable

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

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

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.

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.

Author(s)

Martin Maechler, after conversations with Christophe Dutang.

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

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

Rmpfr documentation built on Nov. 18, 2024, 3:01 p.m.