sumBinomMpfr {Rmpfr}R Documentation

(Alternating) Binomial Sums via Rmpfr

Description

Compute (alternating) binomial sums via high-precision arithmetic. If sBn(f,n):=sBn(f, n) :=sumBinomMpfr(n, f), (default alternating is true, and n0 = 0),

sBn(f,n)=k=n0n(1)(nk)(nk)f(k)=Δnf,sBn(f,n) = \sum_{k = n0}^n (-1)^(n-k) {n \choose k}\cdot f(k) = \Delta^n f,

see Details for the nn-th forward difference operator Δnf\Delta^n f. If alternating is false, the (1)(nk)(-1)^(n-k) factor is dropped (or replaced by 11) above.

Such sums appear in different contexts and are typically challenging, i.e., currently impossible, to evaluate reliably as soon as nn is larger than around 507050--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 kk 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)sB(f,n) := sumBinom(n, f, n0=0) is equal to the nn-th forward difference operator Δnf\Delta^n f,

sB(f,n)=Δnf,sB(f,n) = \Delta^n f,

where

Δnf=k=0n(1)nk(nk)f(k),\Delta^n f = \sum_{k=0}^{n} (-1)^{n-k}{n \choose k}\cdot f(k),

is the nn-fold iterated forward difference Δf(x)=f(x+1)f(x)\Delta f(x) = f(x+1) - f(x) (for x=0x = 0).

The current implementation might be improved in the future, notably for the case where sB(f,n)=sB(f,n)=sumBinomMpfr(n, f, *) is to be computed for a whole sequence n=1,,Nn = 1,\dots,N.

Value

an mpfr number of precision precBits. ss. If alternating is true (as per default),

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

if alternating is false, the (1)k(-1)^k factor is dropped (or replaced by 11) 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")

[Package Rmpfr version 0.9-5 Index]