Rmpfr: sumBinomMpfr – R documentation

Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!

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) choose(n,k) * f(k)) = Δ^n f,

see Details for the n-th forward difference operator Δ^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 Δ^n f,

sB(f,n) = Δ^n f,

where

Delta^n f = sum(k = n0:n ; (-1)^(n-k) choose(n,k) * f(k)),

is the n-fold 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.

Value

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.

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.

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

R MPFR - Multiple Precision Floating-Point Reliable

v0.8-4

GPL (>= 2)

Authors

Martin Maechler [aut, cre] (<https://orcid.org/0000-0002-8685-9910>), Richard M. Heiberger [ctb] (formatHex(), *Bin, *Dec), John C. Nash [ctb] (hjkMpfr(), origin of unirootR()), Hans W. Borchers [ctb] (optimizeR(*, "GoldenRatio"); origin of hjkMpfr())

Initial release

2021-04-08