# sumBinomMpfr: (Alternating) Binomial Sums via Rmpfr In Rmpfr: R MPFR - Multiple Precision Floating-Point Reliable

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

 ```1 2``` ```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, 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.

`chooseMpfr`, `chooseZ` from package gmp.
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26``` ```## "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") ```