Running ‘arith-ex.R’ [7s/9s]
  Running ‘binomial-etc.R’ [2s/3s]
  Running ‘bit-repr.R’ [1s/2s]
  Comparing ‘bit-repr.Rout’ to ‘bit-repr.Rout.save’ ... OK
  Running ‘create.R’ [2s/3s]
  Running ‘functionals.R’ [16s/24s]
  Running ‘lowlevel.R’ [141s/198s]
  Running ‘matrix-ex.R’ [2s/3s]
  Running ‘special-fun-ex.R’ [5s/9s]
Running the tests in ‘tests/special-fun-ex.R’ failed.
Complete output:
  > stopifnot(require("Rmpfr"))
  Loading required package: 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
  
  > (doExtras <- Rmpfr:::doExtras())
  [1] FALSE
  > 
  > all.eq.finite <- function(x,y, ...) {
  +     ## x = 'target'   y = 'current'
  +     if(any(is.finite(y[!(fx <- is.finite(x))])))
  + 	return("current has finite values where target has not")
  +     if(any(is.finite(x[!(fy <- is.finite(y))])))
  + 	return("target has finite values where current has not")
  +     ## now they have finite values at the same locations
  +     all.equal(x[fx], y[fy], ...)
  + }
  > n <- 1000
  > head(x <- mpfr(0:n, 100) / n)
  6 'mpfr' numbers of precision  100   bits 
  [1]                                   0 0.001000000000000000000000000000001
  [3] 0.002000000000000000000000000000002 0.003000000000000000000000000000001
  [5] 0.004000000000000000000000000000003 0.005000000000000000000000000000002
  > 
  > stopifnot(range(x) == 0:1
  + 	  ,all.equal(as.numeric(j0(x)),
  + 		     besselJ(as.numeric(x), 0), tol = 1e-14)
  + 	  ,all.equal(as.numeric(j1(x)),
  + 		     besselJ(as.numeric(x), 1), tol = 1e-14)
  + 	  ,all.equal(as.numeric(y0(x)),
  + 		     besselY(as.numeric(x), 0), tol = 1e-14)
  + 	  ,all.equal(as.numeric(y1(x)),
  + 		     besselY(as.numeric(x), 1), tol = 1e-14)
  + 	  )
  > 
  > ### pnorm() -> erf() :
  > u <- 7*x - 2
  > stopifnot(all.equal(pnorm(as.numeric(u)),
  + 		    as.numeric(pnorm(u)), tol = 1e-14))
  > ## systematic random input testing:
  > set.seed(101)
  > if(doExtras) {
  +     nSim <- 50
  +     n2 <- 100
  + } else {
  +     nSim <- 10
  +     n2 <- 64
  + }
  > for(n in 1:nSim) {
  +     N <- rpois(1, lambda=n2)
  +     N3 <- N %/% 3
  +     x <- c(rnorm(N-N3), 10*rt(N3, df=1.25))# <- some large values
  +     m <- rnorm(N, sd = 1/32)
  +     s <- rlnorm(N, sd = 1/8)
  +     cEps <- .Machine$double.eps
  +     for(LOG in c(TRUE,FALSE))
  + 	for(L.T in c(TRUE,FALSE)) {
  + 	    p. <- pnorm( x, m=m,sd=s, log.p=LOG, lower.tail=L.T)
  + 	    stopifnot(all.equal(p., pnorm(mpfr(x, precBits= 48), m=m,sd=s,
  +                                           log.p=LOG, lower.tail=L.T),
  + 				tol = 128 * cEps))
  + 	    stopifnot(all.equal(p., pnorm(mpfr(x, precBits= 60), m=m,sd=s,
  +                                           log.p=LOG, lower.tail=L.T),
  + 				tol = 2 * cEps))
  + 	}
  +     cat(".")
  + };cat("\n")
  ..........
  > proc.time()
     user  system elapsed 
    3.375   0.087   6.098 
  > 
  > 
  > ### Riemann's Zeta function:
  > 
  > ## -- integer arguments --
  > stopifnot(all(mpfrIs0(zeta(-2*(1:100)))))
  > 
  > k.neg <- 2*(-100:0) - 1
  > Z.neg <- zeta(k.neg)
  > plot(k.neg, abs(as.numeric(Z.neg)), type = "l", log="y")
  > 
  > Pi <- Const("pi", 128L)
  > 
  > ## confirm published value of Euler's gamma to 100 digits
  > pub.g <-
  +     paste("0.5772156649", "0153286060", "6512090082", "4024310421", "5933593992",
  + 	  "3598805767", "2348848677", "2677766467", "0936947063", "2917467495",
  + 	  sep="")
  > 
  > ## almost =
  > our.g <- Const("gamma", log2(10) * 100) # 100 digits
  > (ff.g <- .mpfr2str(our.g))
  $str
  [1] "5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495"
  
  $exp
  [1] 0
  
  $finite
  [1] TRUE
  
  $is.0
  [1] FALSE
  
  > 
  > 
  > M <- function(x) mpfr(x, 128L)
  > stopifnot(all.equal(zeta( 0), -1/2,      tol = 2^-100)
  + 	  , all.equal(zeta(-1), -1/M(12),  tol = 2^-100)
  + 	  , all.equal(zeta(-3),  1/M(120), tol = 2^-100)
  + 	  ## positive ones :
  + 	  , all.equal(zeta(2),  Pi^2/6,   tol = 2^-100)
  + 	  , all.equal(zeta(4),  Pi^4/90,  tol = 2^-100)
  + 	  , all.equal(zeta(6),  Pi^6/945, tol = 2^-100)
  + 	  )
  > 
  > ### Exponential Integral Ei(.)
  > curve(Ei, 0,5, n=5001)
  > if(mpfrVersion() >= "3") { ## only available since MPFR 3.0.0
  +   ### Airy function Ai(.)
  +   curve(Ai, -10, 5, n=5001); abline(h=0,v=0, col="gray", lty=3)
  + }
  > 
  > ### Utilities  hypot(), atan2() : --- TODO !
  > 
  > ## beta(), lbeta()
  > ## ---------------
  > ## The simplistic "slow" versions:
  > B  <- function(a,b) { a <- as(a, "mpfr"); b <- as(b, "mpfr"); gamma(a)*gamma(b) / gamma(a+b) }
  > lB <- function(a,b) { a <- as(a, "mpfr"); b <- as(b, "mpfr"); lgamma(a)+lgamma(b) - lgamma(a+b) }
  > 
  > ## For partly *integer* arguments
  > Bi1 <- function(a,b) 1/(a*chooseMpfr(a+b-1, a)) # a must be integer >= 0
  > Bi2 <- function(a,b) 1/(b*chooseMpfr(a+b-1, b)) # b must be integer >= 0
  > 
  > x <- 1:10 + 0 ; (b10 <- mpfr(x, 128L))
  10 'mpfr' numbers of precision  128   bits 
   [1]  1  2  3  4  5  6  7  8  9 10
  > 
  > stopifnot(all.equal(	B(1,b10),  1/x),
  + 	  all.equal(	B(2,b10),  1/(x*(x+1))),
  + 	  all.equal( beta(1,b10),  1/x),
  + 	  all.equal( beta(2,b10),  1/(x*(x+1))),
  + 	  TRUE)
  > 
  > x <- -10:10 + 0; X <- mpfr(x, 128L)
  > stopifnot(Bi1(1,X) == (B1x <- Bi2(X,1)),
  + 	  Bi1(2,X) == (B2x <- Bi2(X,2)),
  + 	  Bi1(3,X) == (B3x <- Bi2(X,3)),
  + 	  all.equal(B1x,  1/x,               tol= 4e-16) ,
  + 	  all.equal(B2x,  1/(x*(x+1)),       tol= 8e-16) ,
  + 	  all.equal(B3x,  2/(x*(x+1)*(x+2)), tol=16e-16) ,
  + 	  ## these the "poles" are all odd i.e. result in { +Inf / -Inf / NaN}
  + 	  ## are all "ok" {e.g. 1/(x*(x+1)) gives (-Inf, Inf) for x = -1:0 }
  + 	  all.eq.finite(beta(1,X),  1/x) ,
  + 	  all.eq.finite(beta(X,2),  1/(x*(x+1))) ,
  + 	  all.eq.finite(beta(3,X),  2/(x*(x+1)*(x+2)), tol=16e-16) ,
  + 	  TRUE)
  > 
  > ## (a,b)  *both* integer, one negative:
  > for(i in (-20):(-1)) {
  +     cat(i,":\n")
  +     a <- mpfr(i, 99)
  +     i1 <- i+1
  +     b. <- seq_len(-i1)
  +     Bab <- beta(a, b.)
  +     stopifnot(is.nan(beta(a, (i1:0))), is.nan(lbeta(a, (i1:0))),
  + 	      all.equal(Bab, Bi2(a, b.),             tol=1e-20),
  + 	      all.equal(lbeta(a, b.), log(abs(Bab)), tol=1e-20))
  + }
  -20 :
  -19 :
  -18 :
  -17 :
  -16 :
  -15 :
  -14 :
  -13 :
  -12 :
  -11 :
  -10 :
  -9 :
  -8 :
  -7 :
  -6 :
  -5 :
  -4 :
  -3 :
  -2 :
  -1 :
  > 
  > ## (a,b) all positive
  > c10 <- b10 + 0.25
  > for(a in c(0.1, 1, 1.5, 2, 20)) {
  +     stopifnot(all.equal( B(a,b10), (bb <-  beta(a, b10))),
  + 	      all.equal(lB(a,b10), (lb <- lbeta(a, b10))),
  + 	      all.equal(lb, log(bb)),
  + 	      all.equal( B(a,c10), (bb <-  beta(a, c10))),
  + 	      all.equal(lB(a,c10), (lb <- lbeta(a, c10))),
  + 	      all.equal(lb, log(bb)),
  + 	      TRUE)
  + }
  > 
  > ## However, the speedup is *not* much (50%) when applied to vectors:
  > stopifnot(validObject(xx <- outer(b10, runif(20))),
  + 	  dim(xx) == c(length(b10), 20),
  + 	  validObject(vx <- as(xx, "mpfr")), class(vx) == "mpfr", is.null(dim(vx)))
  Error in validObject(xx <- outer(b10, runif(20))) : 
    requires numeric/complex matrix/vector arguments
  Calls: stopifnot ... tryCatch -> tryCatchList -> tryCatchOne -> <Anonymous>
  Execution halted