(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

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.

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

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(Rmpfr)
Loading required package: gmp

Attaching package: 'gmp'

The following objects are masked from 'package:base':

    %*%, apply, crossprod, matrix, tcrossprod

C code of R package 'Rmpfr': GMP using 64 bits per limb


Attaching package: 'Rmpfr'

The following objects are masked from 'package:stats':

    dbinom, dnorm, dpois, pnorm

The following objects are masked from 'package:base':

    cbind, pmax, pmin, rbind

> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/Rmpfr/sumBinomMpfr.Rd_%03d_medium.png", width=480, height=480)
> ### Name: sumBinomMpfr
> ### Title: (Alternating) Binomial Sums via Rmpfr
> ### Aliases: sumBinomMpfr
> ### Keywords: arith
> 
> ### ** 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..
   user  system elapsed 
  0.000   0.000   0.001 
> system.time(resMpfr <- sumBin.all  (n.set, f = sqrt)) ## ~ 0.6 seconds
   user  system elapsed 
  0.188   0.000   0.186 
> ## Don't show: 
> stopifnot(
+     all.equal(resMpfr[1:10], res.R[1:10], tolerance=1e-12),
+     all.equal(resMpfr[1:20], res.R[1:20], tolerance=1e-9 ),
+     all.equal(resMpfr[1:30], res.R[1:30], tolerance=1e-6 ))
> ## End(Don't show)
> 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")
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>