Running 'arith-ex.R' [174s]
  Running 'binomial-etc.R' [1s]
  Running 'bit-repr.R' [1s]
  Comparing 'bit-repr.Rout' to 'bit-repr.Rout.save' ... OK
  Running 'create.R' [2s]
  Running 'functionals.R' [13s]
  Running 'lowlevel.R' [1s]
  Running 'matrix-ex.R' [156s]
  Running 'special-fun-ex.R' [173s]
  Running 'tstHexBin.R' [1s]
Running the tests in 'tests/arith-ex.R' failed.
Complete output:
  > 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 object is masked from 'package:gmp':
  
      outer
  
  The following objects are masked from 'package:stats':
  
      dbinom, dgamma, dnbinom, dnorm, dpois, dt, pnorm
  
  The following objects are masked from 'package:base':
  
      cbind, pmax, pmin, rbind
  
  > ## includes ("gmp")# want to check "mixed arithmetic" too  __ TODO __
  > 
  > `%=N=%` <- function(x,y) (x == y) | (is.na(x) & is.na(y))
  > all.EQ <- function(x,y, tolerance = 2^-98, ...) # very small tol. for MPFR
  +     all.equal(x, y, tolerance=tolerance, ...)
  > warningI <- function(...) warning(..., immediate. = TRUE)
  > 
  > unlist(.Platform)
       OS.type     file.sep   dynlib.ext          GUI       endian      pkgType 
     "windows"          "/"       ".dll"      "RTerm"     "little" "win.binary" 
      path.sep       r_arch 
           ";"        "x64" 
  > 
  > ## Check that we got the "which.*" methods also from "bigq":
  > bcl <- c("ANY", "bigq", "bigz", "mpfr")
  > ##if(packageVersion("gmp") >= "0.5-8") {
  > stopifnot(identical(bcl,
  + 		    sort(unlist(findMethods("which.max")@signatures))),
  + 	  identical(bcl,
  + 		    sort(unlist(findMethods("which.min")@signatures))))
  > ##}
  > 
  > options(warn = 1)# warnings *immediately*
  > (doExtras <- Rmpfr:::doExtras())
  [1] FALSE
  > eps2 <-   2 * .Machine$double.eps
  > eps8 <-   8 * .Machine$double.eps
  > eps32 <- 32 * .Machine$double.eps
  > 
  > ## must take the *larger* of the two precisions:
  > stopifnot(substr(format(mpfr(1, 60)/mpfr(7, 160)), 1,51) == # format() may show more digits
  +           "0.1428571428571428571428571428571428571428571428571")# again has extra "2" at end
  > 
  > (x <- mpfr(0:7, 100) / 7)
  8 'mpfr' numbers of precision  100   bits 
  [1]                                  0 0.14285714285714285714285714285711
  [3] 0.28571428571428571428571428571423 0.42857142857142857142857142857154
  [5] 0.57142857142857142857142857142846 0.71428571428571428571428571428538
  [7] 0.85714285714285714285714285714308                                  1
  > ix <- x^-1000
  > iX <- asNumeric(ix)
  > 
  > stopifnot( mpfrIs0(x - x), # badly failed on 64-bit
  + 	  identical(-x, 0-x),# testing "- x"
  +           all.equal(ix, (1/x)^1000, tol= 1e-25),
  +           is.numeric(iX), iX[1:4] == Inf, # failed previously as we used RNDD (downward rounding)
  +           all.equal(log(iX[5:8]), c(559.6157879, 336.4722366, 154.1506798, 0),
  +                     tol = 1e-9))
  > 
  > ## checking hexadecimal input :
  > stopifnot(mpfr("0xFFFFFFFFFFFFFFFFFFFF", base=16) + 1 == 2^80,
  + ## sign(0) == 0:
  +           identical(sign(as(-1:1, "mpfr")), -1:1 + 0))
  > 
  > stopifnot(all.equal(as.numeric(x+ 1L),
  +                     as.numeric(x)+1L, tol = eps2),
  + 	  as.integer(  x [x < 1]) == 0,# was *wrong* {we round()ed; previously "down"!}
  + 	  as.integer((-x)[x < 1]) == 0,#  (ditto)
  +           (3 * x)/3 <= x,
  +           all.equal(as.numeric(x * 2L),
  +                     as.numeric(x + x), tol = 0))
  > 
  > u <- mpfr(0:17, 128)/17
  > two <- mpfr(2,100)
  > stopifnot(all.EQ(u ^ two, u ^ 2),
  +           identical(u ^ 2, u ^ 2L),
  +           all.EQ(two ^ u, 2 ^ u),
  +           identical(2 ^ u, 2L ^ u),
  +           floor  (3*u) == floor  (3/17*(0:17)),
  +           ceiling(u*5) == ceiling(5/17*(0:17))
  +           )
  > 
  > i7  <- mpfr(0:7,  200)/ 7
  > i17 <- mpfr(0:17, 300)/17
  > stopifnot(all.equal(as.numeric(x+1),
  + 		    as.numeric(x)+1),
  + 	  all.equal(round(x,2), round(asNumeric(x), 2), tol=1e-15),
  + 	  all.equal(round(mpfr(1.152, 80), 2), 1.15), # was wrong {as.integer() bug}
  + 	  all.equal(0:7,   7 * round ( i7, 25), tol = 2e-25),
  + 	  all.equal(0:7,   7 * round ( i7, 50), tol = 2e-50),
  + 	  all.equal(0:17, 17 * signif(i17,100), tol = 2e-100),
  + 	  all.equal(0:17, 17 * signif(i17, 20), tol = 2e-20)
  +           )
  > 
  > ## When we compute with 100 bits,
  > ## we should compare relative errors with  2^-100 :
  > del <- abs((x+pi)-pi - x) / 2^-100
  > stopifnot(del <= 4) ## <= 2 already
  > (fd <- format(del, drop0 = TRUE))
  [1] "0"    "1.75" "0.5"  "1"    "1"    "1"    "2"    "0"   
  > stopifnot(all.equal(as.numeric(del),
  + 		    as.numeric(fd), tol = 1e-15))
  > if(print(Sys.info()[["machine"]]) == "x86_64")
  +     stopifnot(fd %in% as.character(c(0:2, c(2,7)/4)))
  [1] "x86-64"
  > 
  > 
  > checkPmin <- function(x, nx = as(x, "numeric")) {
  +     rx <- if(is(x,"mpfr")) round(x, 25) else x
  +     isZ <- is(x, "bigz") || is(nx, "bigz")
  +     M.X <- max(x, na.rm=TRUE)
  +     m.x <- min(x, na.rm=TRUE)
  +     stopifnot(all.equal(x, nx),
  + 	      pmin(x, x, M.X) %=N=% x, x %=N=% pmax(x, m.x, x),
  + 	      all.equal(x, pmin(x, nx, x, M.X)),
  + 	      all.equal(x, pmax(m.x, nx, x, rx, m.x)),
  + 	      if(isZ)TRUE else all.equal(pmin(x, 0.75), pmin(nx, 0.75)),
  + 	      if(isZ)TRUE else all.equal(pmax(x, 0.25), pmax(nx, 0.25)))
  + }
  > 
  > x <- mpfr(0:7, 100) / 7
  > checkPmin(x)
  > 
  > nx <- (0:7)/7
  > (qx <- as.bigq(0:7, 7))
  Big Rational ('bigq') object of length 8:
  [1] 0   1/7 2/7 3/7 4/7 5/7 6/7 1  
  >  x[c(2,5)] <- NA
  > nx[c(2,5)] <- NA
  > qx[c(2,5)] <- NA
  > 
  > Z <- as.bigz(1:7)
  > mZ <- mpfr(Z, 64)
  > stopifnot(Z == mZ, mZ == Z)
  > 
  > checkPmin(x, nx)
  >     cat("checking pmin(. bigq ): ")
  checking pmin(. bigq ): >     ## FIXME checkPmin(x, qx); cat("[Ok]\n")
  >     ##
  >     print( base::pmin(Z, Z, max(Z)) )# via  gmp:: rep.bigz(x, length.out = *)
  Big Integer ('bigz') object of length 7:
  [1] 1 2 3 4 5 6 7
  >     cat("checking pmin(. bigz )
  +  [currently with lots of pmin() and pmax(...) warnings 'incompatible methods]:\n ")
  checking pmin(. bigz )
   [currently with lots of pmin() and pmax(...) warnings 'incompatible methods]:
   >     checkPmin(Z); cat("[Ok]\n") # via gmp:: all.equal.bigz()
  [Ok]
  > 
  > stopifnot(all.equal( round(x, 10),  round(nx, 10)),
  +           all.equal(signif(x, 10), signif(nx, 10)))
  > 
  > ## L & x ,  x & L  failed in Rmpfr 0.2* and 0.4-2
  > stopifnot(identical(L <- x > 0.5, L & x),
  + 	  identical(L, x & L),
  + 	  identical(x > 0, x | L))
  > 
  > ## Summary() methods {including NA
  > stopifnot(exprs = {
  +     is.na(min(x))
  +     is.na(max(x))
  +     is.na(range(x))
  +     is.na(sum(x))
  +     is.na(prod(x))
  +     min(x, na.rm=TRUE) == 0
  +     max(x, na.rm=TRUE) == 1
  +     range(x, na.rm=TRUE) == 0:1
  +     all.equal(sum (x, na.rm=TRUE)*7,  2+3+5+6+7,  tolerance = 1e-28) # 1.0975e-30
  +     prod(x, na.rm=TRUE) == 0
  +     all.equal(180, prod(x[-1], na.rm=TRUE)*7^4, tol = 1e-15) # 1.579e-16
  +     ##
  +     ## all(), any()  had memory bug [PROTECT missing, but more, somehow]
  +     !all(x)
  +     is.na( all(x[-1]) )
  +     any(x)
  +     is.na(any(x[c(2,5)]))
  +     ## do these *twice* {that triggered R-forge bug  #6764 }
  +     ! all(x, na.rm=TRUE)
  +       any(x, na.rm=TRUE)
  +     ##
  +     ! all(x, na.rm=TRUE)
  +       any(x, na.rm=TRUE)
  + })
  > 
  > ##-------------- Modulo and "integer division" -------------
  > 
  > ## R's	?Arithmetic :
  > ##
  > ##   ‘%%’ indicates ‘x mod y’ and ‘%/%’ indicates integer division.  It
  > ##   is guaranteed that ‘x == (x %% y) + y * ( x %/% y )’ (up to
  > ##   rounding error) unless ‘y == 0’ where the result of ‘%%’ is
  > ##   ‘NA_integer_’ or ‘NaN’ (depending on the ‘typeof’ of the
  > ##   arguments).
  > ##
  > ## and has 'details' about  how	 non-integer 'y' works
  > ##
  > (N <- if(doExtras) 1000 else 200)
  [1] 200
  > (todays.seed <- eval(parse(text=Sys.Date())))# so this is reproducible
  [1] 19805
  >                                         # (and constant within one day)
  > set.seed(todays.seed)
  > mm <- c(-4:4, sample(50, N-9, replace=TRUE))
  > for(n in seq_len(N)) {
  +     cat("."); if(n %% 50 == 0) cat(n,"\n")
  +     m <- mm[n]
  +     prec <- sample(52:200, 1)# "high precision" ==> can use small tol
  +     x <- sample(100, 50) - 20
  +     for(kind in c('int','real')) {
  + 	if(kind == "real") {
  + 	    m <- jitter(m)
  + 	    x <- jitter(x)
  + 	    tol.1 <- eps32 * pmax(1, 1/abs(m))
  + 	    EQ <- function(x,y, tol = tol.1)
  + 		isTRUE(all.equal(x, as.numeric(y), tol=tol))
  + 	    EQ2 <- function(x,y, tol = tol.1) {
  + 		## for the DIV - MOD identity, a small x leads to cancellation
  + 		all((x %=N=% y) | abs(x - y) < tol*pmax(abs(x), 1)) ||
  + 		    isTRUE(all.equal(x, as.numeric(y), tol=tol))
  + 	    }
  + 	} else { ## "integer"
  + 	    EQ2 <- EQ <- function(x,y, tol) all(x %=N=% y)
  + 	}
  + 	i.m <- mpfr(x, prec) %% mpfr(m, prec)
  + 	if(!EQ2(x %% m, i.m)) {
  + 	    cat("\n -- m = ",m,"  (prec = ",prec,")\n")
  + 	    rE <- range(rel.E <- as.numeric(1 - (x %% m)/i.m))
  + 	    print(cbind(x, 'R.%%' = x %% m, rel.E))
  + 	    MSG <- if(max(abs(rE)) < 1e-10) warningI else stop
  + 	    MSG(sprintf("not all equal: range(rel.Err.) = [%g, %g]", rE[1],rE[2]))
  + 	}
  + 	##
  + 	if(m != 0) {
  + 	    ##---Check the    x	 ==  (x %% m) +	 m * ( x %/% m )  assertion ------
  + 	    ##
  + 	    if(EQ2(x, (x %% m) + m*( x %/% m ), tol = 1e-12)) { ## ok for R
  + 		## --> also ok for mpfr ?
  + 		iDm <- mpfr(x, prec) %/% mpfr(m, prec)
  + 		rhs <- i.m + m*iDm
  + 		if(!EQ2(x, i.m + m*iDm)) {
  + 		    cat("\n -- m = ",m,"  (prec = ",prec,")\n")
  + 		    print(cbind(x,' MPFR[ x%%m + m(x %/% m) ]' = as.numeric(rhs), rel.E))
  + 		    MSG <- if(max(abs(rE)) < 1e-10) warningI else stop
  + 		    MSG(sprintf("Identity(MOD - DIV) not all eq.: range(rel.Err.) = [%g, %g]",
  + 				rE[1],rE[2]))
  + 		}
  + 	    } else {
  + 		cat("\n hmm.. the basic	 %% <-> %/%  assertion 'fails' in *R* :\n")
  + 		rhs <- (x %% m)	 +  m * ( x %/% m )
  + 		rel.E <- (1 - rhs/x)
  + 		print(cbind(x, 'x%%m + m(x %/% m)' = rhs, rel.E))
  + 	    }
  + 	}
  +     }
  + }
  ..................................................50 
  ..................................................100 
  ..................................................150 
  ..................................................200 
  > 
  > ## mpfr  o  <number>  now implemented, for  '%%', too :
  > r <- as.double(i <- -10:20)
  > 
  > stopifnot(
  +     ## %% -------------------------------------
  +     mpfr(i, prec=99) %% 7  == i %% 7
  +     , ##
  +     mpfr(i, prec=99) %% 7  ==
  +     mpfr(i, prec=99) %% 7L
  +     , ##
  +     i %% mpfr(27, prec=99) == i %% 27
  +     , ##
  +     r %% mpfr(27, prec=99) == r %% 27
  +     , ## %/% -------------------------------------
  +     mpfr(i, prec=99) %/% 7  == i %/% 7
  +     , ##
  +     mpfr(i, prec=99) %/% 7  ==
  +     mpfr(i, prec=99) %/% 7L
  +     , ##
  +     mpfr(i, prec=99) %/% mpfr(27, prec=99) == i %/% 27
  +     , ##
  +     i %/% mpfr(27, prec=99) == i %/% 27
  +     , ##
  +     i %/% mpfr(27, prec=99) ==
  +     r %/% mpfr(27, prec=99)
  +     , TRUE ##
  + )
  > 
  > cat('Time elapsed: ', proc.time(),'\n') # "stats"
  Time elapsed:  4.56 0.14 4.68 NA NA 
  > 
  > ## Was reproducible BUG in Rmpfr-addition (on Linux, MPFR 4.x.y) -- ## but the bug was Rmpfr,
  > ## in ../src/Ops.c, detecting if *integer*, i.e., long can be used
  > dn <- 1e20
  > dOO <- 9223372036854775808; formatC(dOO) # "9.2...e18"
  [1] "9.223e+18"
  > (r <- dn / (dn + dOO)) # 0.915555  (double prec arithmetic)
  [1] 0.915555
  > ## but *so* strange when switching to Rmpfr :  addition accidentally *subtract*!!
  > n <- mpfr(dn, precBits = 99)
  > (rM <- n / (n + dOO)) # wrongly gave " 1 'mpfr' .... 99 bits; 1.101605140483951.....
  1 'mpfr' number of precision  99   bits 
  [1] 0.9155549598510603318013907353116
  > stopifnot(exprs = {
  +     all.equal(n + dOO, dn + dOO)
  +     all.equal(n / (n + dOO), r)
  + })
  > 
  > ##  log(., base)  :
  > (ten40 <- as.bigz(10)^40)
  Big Integer ('bigz') :
  [1] 10000000000000000000000000000000000000000
  > ten40m <- mpfr(ten40)
  > (lt40 <- log(ten40m, 10)) # gave Error in ... : base != exp(1) is not yet implemented
  1 'mpfr' number of precision  133   bits 
  [1] 40
  > ##  'mpfr' ..  133 bits \\ [1] 40
  > stopifnot(exprs = {
  +     grepl("^40[.]000+$", print(format(lt40, digits = 60)))
  +     identical(lt40, log10(ten40m))
  +     identical(log(ten40m, 2), log2(ten40m))
  +     inherits(Pi <- Const("pi", 140), "mpfr")
  +     all.equal(show(log(ten40m, Pi)),
  +                    log(ten40m)/log(Pi), tol = 1e-40)
  + })
  [1] "40.0000000000000000000000000000000000000000000000000000000000"
  1 'mpfr' number of precision  133   bits 
  [1] 80.45863470352243755058888189154834786778
  > 
  > 
  > 
  > ###------Standard Statistics Functions --------------------------------------------------------
  > 
  > x <- c(del, 1000)
  > stopifnot(identical(mean(x), mean(x, trim=0)))
  > for(tr in (0:8)/16)
  +     stopifnot(all.equal(mean(          x,  trim = tr),
  +                         mean(asNumeric(x), trim = tr), tol=1e-15))
  > 
  > cat('Time elapsed: ', proc.time(),'\n') # "stats"
  Time elapsed:  4.65 0.14 4.78 NA NA 
  > 
  > proc.time()
     user  system elapsed 
     4.65    0.14    4.78 
Running the tests in 'tests/matrix-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 object is masked from 'package:gmp':
  
      outer
  
  The following objects are masked from 'package:stats':
  
      dbinom, dgamma, dnbinom, dnorm, dpois, dt, pnorm
  
  The following objects are masked from 'package:base':
  
      cbind, pmax, pmin, rbind
  
  > 
  > x <- mpfr(0:7, 64)/7
  > mx <- x
  > dim(mx) <- c(4,2)
  > (m. <- mx) # "print"
  'mpfrMatrix' of dim(.) =  (4, 2) of precision  64   bits 
       [,1]                    [,2]                   
  [1,]                      0. 0.571428571428571428564
  [2,] 0.142857142857142857141 0.714285714285714285691
  [3,] 0.285714285714285714282 0.857142857142857142873
  [4,] 0.428571428571428571436  1.00000000000000000000
  > m.[,2] <- Const("pi", 80)
  > m.[,] <- exp(mpfr(1, 90))
  > stopifnot(is(mx, "mpfrMatrix"), dim(mx) == c(4,2),
  +           is(m., "mpfrMatrix"), dim(m.) == dim(mx),
  + 	  dim(is.finite(mx)) == dim(mx),
  + 	  dim(is.nan(mx)) == dim(mx),
  +           getPrec(m.) == 90)
  > 
  > xx <- (0:7)/7
  > m.x <- matrix(xx, 4,2)
  > m2 <- mpfr(xx, 64); dim(m2) <- dim(m.x)
  > ##
  > u <- 10*(1:4)
  > y <- 7 * mpfr(1:12, 80)
  > my <- y
  > dim(my) <- 3:4
  > m.y <- asNumeric(my)
  > stopifnot(all.equal(m2, mpfr(m.x, 64), tol=0), # not identical(..)
  + 	  my[2,2] == 35,
  + 	  my[,1] == 7*(1:3))
  > 
  > .N <- function(x) { if(!is.null(dim(x))) as(x,"array") else as(x,"numeric") }
  > noDN <- function(.) { dimnames(.) <- NULL ; . }
  > allEQ <- function(x,y) all.equal(x,y, tol=1e-15)
  > 
  > ## FIXME write "functions" that take  x -> {mx , m.x}  and run the following as *function*
  > ## ----  then use previous case *and* cases with NA's !
  > ##       and use  higher precision via  fPrec = 2 etc ...
  > 
  > stopifnot(allEQ(m.x, noDN(.N(mx))),
  + 	  allEQ(m.y, noDN(.N(my))),
  + 	  allEQ(noDN(.N(my %*% mx)), m.y %*% m.x),
  + 	  allEQ(noDN(.N(crossprod(mx, t(my)))), crossprod(m.x, t(m.y))),
  + 	  allEQ(noDN(.N(tcrossprod(my, t(mx)))),
  + 			tcrossprod(m.y, t(m.x))),
  + 	  ##
  + 	  identical(mx, t(t(mx))),
  + 	  identical(my, t(t(my))),
  + 	  ## matrix o vector .. even  vector o vector
  + 	  identical(noDN(.N(my %*% 1:4)), m.y %*% 1:4 ),
  + 	  identical(noDN(.N(my %*% my[2,])), m.y %*% .N(my[2,])),
  + 	  identical( crossprod(1:3, my), 1:3 %*%   my),
  + 	  identical(tcrossprod(1:4, my), 1:4 %*% t(my)),
  + 	  identical(crossprod(y), t(y) %*% y),
  + 	  identical(tcrossprod(y), y %*% t(y)),
  + 	  identical(noDN(.N( crossprod(y))), crossprod(7 * 1:12)),
  + 	  identical(noDN(.N(tcrossprod(y))),tcrossprod(7 * 1:12)),
  + 	  identical(tcrossprod(1:3, u), noDN(.N(tcrossprod(1:3, as(u,"mpfr")))))
  + 	  )
  > 
  > mx[3,1] <- Const("pi", 64)
  > stopifnot(allEQ(sum(mx[,1]), pi + 4/7))
  > m2 <- mx[c(1,4),]
  > stopifnot(dim(m2) == c(2,2), sum(m2) == 2)
  > 
  > ## "mpfrArray" o "mpfr" :
  > Tmx <- array(TRUE, dim(mx), dimnames=dimnames(mx))
  > stopifnot(identical(Tmx, mx == (mx - mpfr(0, 10))),
  + 	  identical(Tmx, mx - mpfr(1, 10) * mx == 0))
  Note: method with signature 'mpfrArray#mpfr' chosen for function '==',
   target signature 'mpfrMatrix#mpfrMatrix'.
   "mpfr#mpfrArray" would also be valid
  > ## subassignment, many kinds
  > mx[5] <- pi
  > mx[6] <- Const("pi",100)
  > stopifnot(validObject(mx), allEQ(mx[5], mx[6]),
  + 	  getPrec(mx) == c(rep(64,5), 100, 64,64))
  > 
  > ## %*% with vectors on LHS, ...
  > y <- t(2:4) # 1 x 3 matrix
  > m1 <-     (0:10)     %*% y
  > m2 <- mpfr(0:10, 50) %*% y
  > stopifnot((d <- m1 - m2) == 0, identical(dim(m1), dim(d)),
  +           m2 == m1, m1 == m2)
  > 
  > r <- 10*(0:4)
  > y <- t(2:6)
  > m1 <- 1:3 %*% y  %*% r
  > y. <- t(mpfr(2:6, 20))
  > m2 <- 1:3 %*% y. %*% r
  > stopifnot(m1 == m2, m1 - m2 == 0, identical(dim(m1), dim(m2)))
  > 
  > ### Array (non-matrix) ---- indexing & sub-assignment :
  > A <- mpfrArray(1:24, prec = 96, dim = 2:4)
  > a <-     array(1:24,            dim = 2:4)
  > a.1 <- as(A[,,1], "array")
  > a1. <- as(A[1,,], "array")
  > A1. <- as(A[1,,], "mpfr")
  > 
  > stopifnot(all.equal(noDN(a.1), a[,,1], tol=0),
  + 	  identical(A1., as.vector(A[1,,])),
  + 	  ## arithmetic, subsetting etc:
  + 	  allEQ(noDN(.N(A / A1.)), a/c(a1.)),
  + 	  allEQ(noDN(.N(a / A1.)), a/c(a1.)),
  +           identical(noDN(A == 23), a == 23),
  +           identical(noDN(10 >= A), 10 >= a),
  +           identical(noDN(A <=  2), a <=  2),
  +           identical(noDN(A < 2.5), a < 2.5),
  +           identical(noDN(A !=  5), a !=  5),
  +           identical(A != 3, !(3 == A)),
  +           identical(-1 > A, A == 100),
  +           identical(noDN(A <= 0), a == pi)
  + 	  )
  > 
  > A[1,2,3] <- Const("pi")
  > A[1, ,2] <- 1 / A[1,,2]
  > A
  'mpfrArray' of dim(.) =  (2, 3, 4) of precision  96 .. 120  bits 
  , , 1
  
       [,1]                            [,2]                           
  [1,] 1.00000000000000000000000000000 3.00000000000000000000000000000
  [2,] 2.00000000000000000000000000000 4.00000000000000000000000000000
       [,3]                           
  [1,] 5.00000000000000000000000000000
  [2,] 6.00000000000000000000000000000
  
  , , 2
  
       [,1]                             [,2]                            
  [1,] 0.142857142857142857142857142858 0.111111111111111111111111111111
  [2,]  8.00000000000000000000000000000  10.0000000000000000000000000000
       [,3]                             
  [1,] 0.0909090909090909090909090909098
  [2,]   12.0000000000000000000000000000
  
  , , 3
  
       [,1]                            [,2]                                   
  [1,] 13.0000000000000000000000000000 3.1415926535897932384626433832795028847
  [2,] 14.0000000000000000000000000000         16.0000000000000000000000000000
       [,3]                           
  [1,] 17.0000000000000000000000000000
  [2,] 18.0000000000000000000000000000
  
  , , 4
  
       [,1]                            [,2]                           
  [1,] 19.0000000000000000000000000000 21.0000000000000000000000000000
  [2,] 20.0000000000000000000000000000 22.0000000000000000000000000000
       [,3]                           
  [1,] 23.0000000000000000000000000000
  [2,] 24.0000000000000000000000000000
  
  > 
  > ## check that A is "==" a  where
  > a <- array(1:24, 2:4); a[1,2,3] <- pi; a[1,,2] <- 1/a[1,,2]
  > stopifnot(allEQ(noDN(.N(A)), a),
  +           ## check aperm() methods :
  +           allEQ(noDN(.N(aperm(A))), aperm(a)),
  +           {p <- c(3,1:2); allEQ(noDN(.N(aperm(A,p))), aperm(a,p))},
  +           {p <- c(2:1,3); allEQ(noDN(.N(aperm(A,p))), aperm(a,p))})
  > 
  > ## cbind() / rbind():
  > options(warn = 2)## no warnings here - ("exact recycling"):
  > validObject(m0 <- cbind(pi=pi, i = 1:6))
  [1] TRUE
  > validObject(m1 <- cbind(a=Const("pi",60),i = 1:6, "1/mp" = 1/mpfr(1:3,70)))
  [1] TRUE
  > validObject(m2 <- cbind(pi=pi, i = 1:2, 1/mpfr(1:6,70)))
  [1] TRUE
  > validObject(n2 <- rbind(pi=pi, i = 1:2, 1/mpfr(1:6,70)))
  [1] TRUE
  > stopifnot(is(m0,"matrix"), is(m1, "mpfrMatrix"), is(m2, "mpfrMatrix"),
  +           dim(m0) == c(6,2), dim(m1) == c(6,3), dim(m2) == c(6,3))
  > options(warn = 1)
  > suppressWarnings(eval(ex <- quote(m3 <- cbind(I=10, 1:3, inv=1/mpfr(2:3,80)))))
  > validObject(suppressWarnings(     n3 <- rbind(I=10, 1:3, inv=1/mpfr(2:3,80))))
  [1] TRUE
  > stopifnot(identical(t(n2), m2),
  +           identical(t(n3), m3), validObject(m3),
  +           is(tryCatch(eval(ex), warning=function(.).), "warning"),
  + 	  identical(cbind("A", "c"), matrix(c("A", "c"), 1,2)),
  + 	  identical(rbind("A", 2),   matrix(c("A", "2"), 2,1)) )
  > 
  > ## head() / tail() :
  > stopifnot(all.equal(c(21, 12),
  + 		    dim(mm3 <- m3[rep(1:3, each=7), rep(3:1, 4)])),
  + 	  all.equal(dim(t3 <- tail(mm3)), c(6, 12)),
  + 	  all.equal(head(mm3), mm3[1:6,]))
  > 
  > ## matrix(<mpfr>) works since 2015-02-28:
  > x <- mpfr(pi,64)*mpfr(2,64)^(2^(0:19))
  > (mx <- matrix(x, 4,5))
  'mpfrMatrix' of dim(.) =  (4, 5) of precision  64   bits 
       [,1]                   [,2]                      
  [1,] 6.28318530717958623200     205887.416145660681650
  [2,] 12.5663706143591724640     13493037704.5220184326
  [3,] 50.2654824574366898560     57952155664616980480.0
  [4,] 804.247719318987037695 1.06902858406496670619e+39
       [,3]                        [,4]                        
  [1,]  3.63771576891766310100e+77 3.28104443733842154527e+1233
  [2,] 4.21218708934506345775e+154 3.42668632977872042983e+2466
  [3,] 5.64761954589225685791e+308 3.73765172555867836946e+4932
  [4,] 1.01526868859647402324e+617 4.44680197657346404992e+9864
       [,5]                          
  [1,]  6.29427491061341831593e+19728
  [2,]  1.26107681736228617526e+39457
  [3,]  5.06212903659356755184e+78913
  [4,] 8.15673870189463167140e+157826
  > 
  > stopifnot(is(mx, "mpfrMatrix"),
  +     all.equal(matrix(0:19, 4,5),
  +               asNumeric(log2(log2(mx) - log2(Const("pi")))),
  +               tol = 1e-15)) # 64b-lnx: see 8.1e-17
  > 
  > 
  > ## Ensure that apply() continues to work with 'bigz'/'bigq'
  > A <- matrix(2^as.bigz(1:12), 3,4)
  > Aq <- as.bigq(A)
  > mA <- as(A, "mpfr") # failed {as dim(A) is "double", not "integer"}
  > (Qm <- A / (A^2 - 1)) # big rational matrix
  Big Rational ('bigq') 3 x 4 matrix:
       [,1] [,2]    [,3]       [,4]         
  [1,] 2/3  16/255  128/16383  1024/1048575 
  [2,] 4/15 32/1023 256/65535  2048/4194303 
  [3,] 8/63 64/4095 512/262143 4096/16777215
  > MQ <- mpfr(Qm, precBits = 64)
  > stopifnot(exprs = {
  +     mA == A
  +     mA == Aq
  +     is.bigq(Aq)
  +     mA == mpfr(A, precBits=16)
  +     mA == asNumeric(A)
  +     is.bigq(Qm)
  +     is(MQ, "mpfrMatrix")
  +     all.equal(Qm, MQ, tol = 1e-18)
  +     identical(dim(mA), dim(A))
  +     identical(dim(mA), dim(Qm))
  +     identical(asNumeric(apply(A,  1, min)),
  +               apply(asNumeric(A), 1, min))
  +     identical(asNumeric(apply(A,  1, max)),
  +               apply(asNumeric(A), 1, max))
  +     identical(asNumeric(apply(A,  2, max)),
  +               apply(asNumeric(A), 2, max))
  +     identical(asNumeric(apply(A,  2, min)),
  +               apply(asNumeric(A), 2, min))
  + })
  Warning in getPrec(x) :
    default precision for 'bigq' arbitrarily chosen as 128
  > ## mA etc, failed up to Rmpfr 0.8-1; the apply() parts failed up to Rmpfr 0.6.0
  > 
  > if(FALSE) ## Bug in gmp :  apply(*, range) does *not* return matrix
  + stopifnot(
  +           identical(asNumeric(apply(A,  1, range)),
  +                     apply(asNumeric(A), 1, range))
  + )
  > 
  > if(FALSE) ## Bug in gmp : --- no mean method for bigz, just mean.bigq
  + stopifnot(
  +           all.equal(asNumeric(apply(A,  1, mean)),
  +                     apply(asNumeric(A), 1, mean))
  +           ,
  +           all.equal(asNumeric(apply(A,  2, mean)),
  +                     apply(asNumeric(A), 2, mean))
  + )
  > 
  > cat('Time elapsed: ', proc.time(),'\n') # "stats"
  Time elapsed:  1.29 0.07 1.37 NA NA 
  > if(!interactive()) warnings()
  > 
  > proc.time()
     user  system elapsed 
     1.29    0.07    1.37 
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 object is masked from 'package:gmp':
  
      outer
  
  The following objects are masked from 'package:stats':
  
      dbinom, dgamma, dnbinom, dnorm, dpois, dt, pnorm
  
  The following objects are masked from 'package:base':
  
      cbind, pmax, pmin, rbind
  
  > (doExtras <- Rmpfr:::doExtras())
  [1] FALSE
  > options(nwarnings = 50000, width = 99)
  > 
  > (do.pdf <- !dev.interactive(orNone = TRUE))
  [1] TRUE
  > if(do.pdf) {
  +     pdf.options(width = 8.5, height = 6) # for all pdf plots
  +     pdf("special-fun.pdf")
  + }
  > 
  > 
  > ## to enhance  |rel.Err| plots:  {also in ~/R/Pkgs/DPQ/tests/pow-tst.R }
  > drawEps.h <- function(p2 = -(53:51), side = 4, lty=3, lwd=2, col=adjustcolor(2, 1/2)) {
  +     abline(h = 2^p2, lty=lty, lwd=lwd, col=col)
  +     axis(side, las=2, line=-1, at = 2^p2,
  +          labels = as.expression(lapply(p2, function(p) substitute(2^E, list(E=p)))),
  +          col.axis = col, col=NA, col.ticks=NA)
  + }
  > mtextVersion <- function(adj = 1, col = 1) {
  +     mtext(osVersion, line=1, col=col, adj=adj)
  +     mtext(sfsmisc::shortRversion(spaces=FALSE), col=col, adj=adj)
  + }
  > 
  > 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.0010000000000000000000000000000008
  [3] 0.0020000000000000000000000000000015 0.0030000000000000000000000000000007
  [5]  0.004000000000000000000000000000003 0.0050000000000000000000000000000022
  > 
  > stopifnot(exprs = {
  +     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 
     2.78    0.10    2.87 
  > 
  > ## Jerry Lewis - Aug 2, 2019
  > ## Contrast the results of pnorm with double and mpfr inputs
  > x <- c(1:9, 5*(2:9), 10*(5:20)) ; x <- c(-rev(x), 0, x)
  > pdL <- pnorm(x, log.p=TRUE)
  > pdU <- pnorm(x, log.p=TRUE, lower.tail=FALSE)
  > stopifnot(exprs = {
  +     !is.unsorted(x)
  +     35 %in% x
  +     x == -rev(x) # exactly
  +     pdL == rev(pdU) # even exactly, currently
  + })
  > mx <- mpfr(x, precBits = 128)
  > pmL <- pnorm(mx, log.p=TRUE)
  > pmU <- pnorm(mx, log.p=TRUE, lower.tail=FALSE)
  > stopifnot(exprs = {
  +     pmL < 0 # not true for 'pdL' which underflows
  +     pmL == rev(pmU) # even exactly, currently
  +     all.equal(pmL, pdL, tol=4e-16) # 'tol=0' shows 4.46e-17
  + })
  > ## some explorations :
  > dlp <- diff(log(-pmL))/diff(x)
  > n <- length(x)
  > x.1 <- (x[-1] + x[-n])/2
  > plot(x.1, dlp, type="b", ylab = "d/dx  log(-pnorm(., log=TRUE))"); mtextVersion()
  > plot(x.1[-1], diff(dlp)/diff(x.1), type="b", ylab = "d^2/dx^2  log(-pnorm(., log=TRUE))")
  > stopifnot(exprs = {
  +     -1 < (d2 <- diff(dlp)/diff(x.1))
  +     d2 < 0
  +     diff(d2) < 0
  + })
  > x.3 <- x.1[-c(1L,n-1L)]
  > plot(x.3, -diff(d2)/ diff(x.1)[-1], type="o", log="y")
  > 
  > 
  > 
  > 
  > ### 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] "57721566490153286060651209008240243104215933593992359880576723488486772677766467093694706329174674948"
  
  $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)
  > 
  > if(do.pdf) { dev.off(); pdf("special-fun-beta.pdf") }
  > 
  > 
  > x <- -10:10 + 0; X <- mpfr(x, 128L)
  > stopifnot(exprs = {
  +     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)
  + })
  > 
  > ## (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), allow.logical0 = TRUE)
  + }
  -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)))
  > C1 <- replicate(10, system.time(bb <<- beta(vx, vx+2)))
  > C2 <- replicate(10, system.time(b2 <<-    B(vx, vx+2)))
  > summary(1000*C1[1,]) ##  80.3 {cmath-5, 2009}
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
       10      10      20      18      20      30 
  > summary(1000*C2[1,]) ## 125.1 { " }
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
       10      20      25      24      30      40 
  > stopifnot(all.equal(bb, b2))
  > ## and for a single number, the speedup is a factor 3:
  > x1 <- vx[1]; x2 <- x1+2
  > system.time(for(i in 1:100) bb <- beta(x1, x2))# .27
     user  system elapsed 
     0.07    0.00    0.06 
  > system.time(for(i in 1:100) b2 <-    B(x1, x2))# .83
     user  system elapsed 
     0.17    0.00    0.17 
  > 
  > ## a+b is integer <= 0, but a and b are not integer:
  > a <- b <- .5 + -10:10
  > ab <- data.matrix(expand.grid(a=a, b=b, KEEP.OUT.ATTRS=FALSE))
  > ab <- mpfr(ab[rowSums(ab) <= 0, ], precBits = 128)
  > stopifnot( beta(ab[,"a"], ab[,"b"]) == 0,
  + 	  lbeta(ab[,"a"], ab[,"b"]) == -Inf)
  > ## was  NaN  in Rmpfr <= 0.5-2
  > 
  > stopifnot(all.equal(6 * beta(mpfr(1:3,99), -3.), c(-2,1,-2), tol=1e-20))
  > ## add more checks, notably for b (> 0)  above and below the "large_b" in
  > ## ../src/utils.c :
  > bb <- beta(mpfr(1:23, 128), -23)
  > stopifnot(all.equal(bb, Bi1(1:23, -23), tol=1e-7))
  >                                         # Bi1() does not get high prec for small b
  > ## can be written via rationals:  N / D :
  > bn <- c(330, -360, 468, -728, 1365, -3120, 8840, -31824,
  +         151164, -1007760, 10581480, -232792560)
  > bn <- c(rev(bn[-1]), bn)
  > bd <- 24* as.bigz(2 * 3 * 5 * 7 * 11) * 13 * 17 * 19 * 23
  > stopifnot(all.equal(bb, as(bn/bd,"mpfr"), tol=0))
  > 
  > stopifnot(all.equal(6 * beta(mpfr(1:3,	99), -3.),
  + 			     c(-2,1,-2),	    tol=1e-20),
  + 	  all.equal(   lbeta(mpfr(1:3, 128), -3.),
  + 		    log(mpfr(c( 2,1, 2), 128) / 6), tol=1e-20))
  > 
  > ## add more checks, notably for b (> 0)  above and below the "large_b" in
  > ## ../src/utils.c :
  > bb <- beta(mpfr(1:23, 128), -23)
  > stopifnot(all.equal(bb, Bi1(1:23, -23), tol=1e-7))
  > 					# Bi1() does not get high prec for small b
  > ## can be written via rationals:  N / D :
  > bn <- c(330, -360, 468, -728, 1365, -3120, 8840, -31824,
  + 	151164, -1007760, 10581480, -232792560)
  > bn <- c(rev(bn[-1]), bn)
  > bd <- 24* as.bigz(2 * 3 * 5 * 7 * 11) * 13 * 17 * 19 * 23
  > stopifnot(all.equal(bb, as(bn/bd,"mpfr"), tol=0))
  > 
  > ## 2) add check for 'b' >  maximal unsigned int {so C code uses different branch}
  > two <- mpfr(2, 128)
  > for(b in list(mpfr(9, 128), mpfr(5, 128)^10, two^25, two^26, two^100)) {
  +     a <- -(b+ (1:7))
  +     stopifnot(a+b == -(1:7), # just ensuring that there was no cancellation
  + 	      is.finite( B <-  beta(a,b)), ## was NaN ..
  + 	      is.finite(lB <- lbeta(a,b)), ## ditto
  + 	      all.equal(log(abs(B)), lB),
  + 	      TRUE)
  + }
  > 
  > ee <- c(10:145, 5*(30:59), 10*(30:39), 25*(16:30))
  > b <- mpfr(2, precBits = 10 + max(ee))^ee # enough precision {now "automatic"}
  > stopifnot((b+4)-b == 4, # <==> enough precision above
  + 	  b == (b. <- as(as(b,"bigz"),"mpfr")))
  > (pp <- getPrec(b.))# shows why b. is not *identical* to b.
    [1]  12  12  16  16  16  16  20  20  20  20  24  24  24  24  28  28  28  28  32  32  32  32  36
   [24]  36  36  36  40  40  40  40  44  44  44  44  48  48  48  48  52  52  52  52  56  56  56  56
   [47]  60  60  60  60  64  64  64  64  68  68  68  68  72  72  72  72  76  76  76  76  80  80  80
   [70]  80  84  84  84  84  88  88  88  88  92  92  92  92  96  96  96  96 100 100 100 100 104 104
   [93] 104 104 108 108 108 108 112 112 112 112 116 116 116 116 120 120 120 120 124 124 124 124 128
  [116] 128 128 128 132 132 132 132 136 136 136 136 140 140 140 140 144 144 144 144 148 148 152 156
  [139] 164 168 172 176 184 188 192 196 204 208 212 216 224 228 232 236 244 248 252 256 264 268 272
  [162] 276 284 288 292 296 304 312 324 332 344 352 364 372 384 392 404 428 452 476 504 528 552 576
  [185] 604 628 652 676 704 728 752
  > system.time(Bb <- beta(-b-4, b))# 0.334 sec
     user  system elapsed 
     0.07    0.00    0.08 
  > if(dev.interactive())
  +     plot(ee, asNumeric(log(Bb)), type="o",col=2)
  > lb <- asNumeric(log(Bb))
  > ## using  coef(lm(lb ~ ee))
  > stopifnot(all.equal(lb, 3.175933 -3.46571851*ee, tol = 1e-5))# 4.254666 e-6
  > 
  > 
  > bb <- beta(           1:4,   mpfr(2,99))
  > stopifnot(identical(bb, beta(mpfr(2,99), 1:4)),
  + 	  all.equal((2*bb)*cumsum(1:4), rep(1, 4), tol=1e-20),
  + 	  getPrec(bb) == 128)
  > 
  > 
  > ##-- The d*() density functions from ../R/special-fun.R  |  ../man/distr-etc.Rd ---
  > 
  > if(do.pdf) { dev.off(); pdf("special-fun-density.pdf") }
  > 
  > dx <- 1400+ 0:10
  > mx <- mpfr(dx, 120)
  > nx <- sort(c(c(-32:32)/2, 50*(-8:8)))
  > 
  > xL <- 2^(989+(0:139)/4) # "close" to double.xmax
  > dnbD   <- dnbinom(xL, prob=1-1/4096, size=1e307, log=TRUE)# R's own
  > iF <- -(130:140) # index of finite dnbD[]
  > dnbx8  <- dnbinom(xL, prob=1-mpfr(2, 2^ 8)^-12, size=1e307, log=TRUE)
  > dnbx10 <- dnbinom(xL, prob=1-mpfr(2, 2^10)^-12, size=1e307, log=TRUE)
  > dnbx13 <- dnbinom(xL, prob=1-mpfr(2, 2^13)^-12, size=1e307, log=TRUE)
  > 
  > stopifnot(exprs = {
  +     all.equal(dpois(dx, 1000), dpois(mx, 1000), tol = 3e-13) # 64b Lnx: 7.369e-14
  +     all.equal(dbinom(0:16, 16, pr = 4 / 5),
  +               dbinom(0:16, 16, pr = 4/mpfr(5, 128)) -> db, tol = 5e-15)# 64b Lnx: 4.3e-16
  +     all.equal(dnorm(     -3:3,       m=10, s=1/4),
  +               dnorm(mpfr(-3:3, 128), m=10, s=1/4), tol = 1e-15) # 64b Lnx: 6.45e-17
  +     all.equal(dnorm(nx), dnorm(mpfr(nx, 99)), tol = 1e-15)
  +     all.equal(dnorm(     nx,      m = 4, s = 1/4),
  +               dnorm(mpfr(nx, 99), m = 4, s = 1/4), tol = 1e-15)
  +     all.equal(dnorm(     nx,      m = -10, s = 1/4, log=TRUE),
  +               dnorm(mpfr(nx, 99), m = -10, s = 1/4, log=TRUE), tol = 1e-15)
  +     ## t-distrib. :
  +     all.equal(dt(nx, df=3), dt(mpfr(nx, 99), df=3), tol = 1e-15)
  +     all.equal(dt(     nx,      df = 0.75),
  +               dt(mpfr(nx, 99), df = 0.75), tol = 1e-15)
  +     all.equal(dt(     nx,      df = 2.5, log=TRUE),
  +               dt(mpfr(nx, 99), df = 2.5, log=TRUE), tol = 1e-15)
  +     ## negative binomial  dnbinom():
  +     all.equal(dnbx13, dnbx10, tol = 2^-999) # see 2^-1007, but not 2^-1008
  +     all.equal(dnbx13, dnbx8,  tol = 2^-238) # see 2^-239,  but not 2^-240
  +     all.equal(dnbx10[iF], dnbD[iF], tol = 6e-16) # R's *is* accurate here (seen 2.9e-16)
  + })
  > 
  > 
  > ## plot dt() "error" of R's implementation
  > nx <- seq(-100, 100, by=1/8)
  > dtd <- dt(     nx,        df= .75)
  > dtM <- dt(mpfr(nx,  256), df= .75)
  > if(doExtras) withAutoprint({
  +  system.time(
  +   dtMx <- dt(mpfr(nx, 2048), df= .75) ) # 2.5 sec
  +  stopifnot(all.equal(dtMx, dtM, tol = 2^-254)) # almost all of dtM's 256 bits are correct
  + })
  > relE <- asNumeric(dtd/dtM - 1)
  > plot(relE ~ nx,      type="l", col=2); mtextVersion()
  > plot(abs(relE) ~ nx, type="l", col=2, log="y", ylim=c(5e-17, 1.5e-15))
  > 
  > ## ============== even smaller 'df' such that lgamma1p(df) is better than lgamma(1+df) ====
  > 
  > require(sfsmisc)# -> eaxis(); relErrV()
  Loading required package: sfsmisc
  
  Attaching package: 'sfsmisc'
  
  The following objects are masked from 'package:gmp':
  
      factorize, is.whole
  
  > 
  > u <- sort(outer(10^-(20:1), c(1,2,5))) # *not* "exact" on purpose
  > ## .. unfinished .. exploring *when* dt() would suffer from inaccurate stirlerr()  -- would it?
  > 
  > nu <- 2^-(70:1)
  > dt10  <- dt(     10,        df=nu)
  > dt10M <- dt(mpfr(10, 1024), df=nu)
  > re10 <- asNumeric(relErrV(dt10M, dt10))
  > 
  > plot(re10 ~ nu, type="l", lwd=2, log="x", main = quote(rel.Err( dt(10, df==nu) )),
  +      xaxt="n"); eaxis(1, nintLog=20)
  > mtextVersion()
  > abline(h = (-1:1)*2^-53, lty=4, col=adjustcolor("blue", 1/2))
  > 
  > plot(abs(re10) ~ nu, type="l", lwd=2, log="xy",
  +      xlab = quote(df == nu), ylab = quote(abs(relE)),
  +      main = quote(abs(rel.Err( dt(10, df==nu) ))), xaxt="n", yaxt="n")
  > eaxis(1, nintLog=20); eaxis(2); drawEps.h()
  > 
  > x0 <- c(0, 10^(-5:10)) # only >= 0 should be sufficient; x0 <- c(-rev(x0),0,x0)
  > stopifnot(!is.unsorted(nu), # just for plotting ..
  +           !is.unsorted(x0))
  > xnu <- expand.grid(x=x0, df=nu)
  > dt2  <- with(xnu, dt(     x,       df=df))
  > dtM2 <- with(xnu, dt(mpfr(x, 512), df=df))
  > str(relE2 <- `attributes<-`(asNumeric(relErrV(dtM2, dt2)),
  +                             attr(xnu, "out.attrs")))
   num [1:17, 1:70] 1.78e-15 1.08e-15 1.11e-15 7.85e-16 1.73e-15 ...
   - attr(*, "dimnames")=List of 2
    ..$ x : chr [1:17] "x=0e+00" "x=1e-05" "x=1e-04" "x=1e-03" ...
    ..$ df: chr [1:70] "df=8.470329e-22" "df=1.694066e-21" "df=3.388132e-21" "df=6.776264e-21" ...
  > 
  > ## consistency check that with() etc was fine:
  > stopifnot(identical(re10, unname(relE2[which(x0 == 10), ])))
  > 
  > filled.contour(x=log10(1e-7+x0), y=log10(nu), z = relE2)
  > filled.contour(x=log10(1e-7+x0), y=log10(nu), z = abs(relE2))
  > ## around nu = 10^-16 is the most critical place
  > 
  > (pch <- c(1L:9L, 0L, letters, LETTERS)[1:ncol(relE2)])
   [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "0" "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m"
  [24] "n" "o" "p" "q" "r" "s" "t" "u" "v" "w" "x" "y" "z" "A" "B" "C" "D" "E" "F" "G" "H" "I" "J"
  [47] "K" "L" "M" "N" "O" "P" "Q" "R" "S" "T" "U" "V" "W" "X" "Y" "Z" NA  NA  NA  NA  NA  NA  NA 
  [70] NA 
  > 
  > matplot(x0+1e-7, relE2, type="b", log="x", main="rel.err{  dt(x, df=df) }")
  Warning message:
  In matplot(x0 + 1e-07, relE2, type = "b", log = "x", main = "rel.err{  dt(x, df=df) }") :
    default 'pch' is smaller than number of columns and hence recycled
  > legend("topright", legend = paste0("df=",formatC(nu,wid=3)), ncol=7,
  +        bty="n", lwd=1, pch=pch, col=1:6, lty=1:5, cex = 0.8)
  > abline(h = c(-4:4)*2^-53, lty=3, col="gray")
  > 
  > matplot(nu, t(relE2), type="b", log="x", main="rel.err{  dt(x, df=df) }")
  > legend("topright", legend = paste0("x=",formatC(x0,wid=3)), ncol=7,
  +        bty="n", lwd=1, pch=pch, col=1:6, lty=1:5, cex = 0.8)
  > abline(h = c(-4:4)*2^-53, lty=3, col="gray")
  > 
  > matplot(nu, pmax(abs(t(relE2)), 1e-19), type="b", log="xy", axes=FALSE, ylab = quote(abs("rel Err")),
  +         ylim = c(7e-17, max(abs(relE2))), main="|rel.err{ dt(x, df=df)}|")
  > eaxis(1, nintLog=22) ; eaxis(2, line=-1/2); drawEps.h()
  > legend("topright", legend = paste0("x=",formatC(x0,wid=3)), ncol=7,
  +        bty="n", lwd=1, pch=pch, col=1:6, lty=1:5, cex = 0.8)
  > 
  > 
  > 1
  [1] 1
  > ## dnbinom() -- has mode as expected, but with huge size, the scales are "off reality" ..
  > 
  > ### ..... TODO !
  > 
  > ### dgamma(): ----------------------------------------------------
  > if(do.pdf) { dev.off(); pdf("special-fun-dgamma.pdf") }
  > 
  > xe <- c(-2e5, -1e5, -2e4, -1e4, -2000, -1000, -500, -200, -100, -50, -20, -10)
  > (xe <- c(xe, -8:8, -rev(xe)))
   [1] -2e+05 -1e+05 -2e+04 -1e+04 -2e+03 -1e+03 -5e+02 -2e+02 -1e+02 -5e+01 -2e+01 -1e+01 -8e+00
  [14] -7e+00 -6e+00 -5e+00 -4e+00 -3e+00 -2e+00 -1e+00  0e+00  1e+00  2e+00  3e+00  4e+00  5e+00
  [27]  6e+00  7e+00  8e+00  1e+01  2e+01  5e+01  1e+02  2e+02  5e+02  1e+03  2e+03  1e+04  2e+04
  [40]  1e+05  2e+05
  > two <- mpfr(2, 64)
  > ## For centering at E[.], will use xP(x, shp) :
  > xP <- function(x, d) x - d*(x > d)
  > aEQformat <- function(xy, ...) format(xy, digits = 7, ...)
  > allEQ_0 <- function (target, current, ...)
  +     all.equal(target, current, tolerance = 0, formatFUN = aEQformat, ...)
  > stopIfNot <-
  +     if("allow.logical0" %in% names(formals(stopifnot))) { # experimental (MM only)
  +         stopifnot
  +     } else function(exprs, allow.logical0) stopifnot(exprs=exprs)
  > 
  > for(shp in c(2^c(-20, -3, -1:1, 4, 10, 50))) {
  +     cat("shape = 2^", log2(shp), ":\n-------------\n")
  +     d.dg  <- dgamma(xP(2 ^ xe, shp), shape=shp)
  +     m.dg  <- dgamma(xP(two^xe, shp), shape=shp)
  +     m.ldg <- dgamma(xP(two^xe, shp), shape=shp, log=TRUE)
  +     stopIfNot(exprs = {
  +         !is.unsorted(xe)
  +         is.finite(m.dg)
  +         m.dg >= 0
  +         shp > 1  || all(diff(m.dg) <= 0)
  +         shp > 100|| all((m.dg > 0) >= (d.dg > 0))
  +         any(fin.d <- is.finite(d.dg))
  +         m.dg[!fin.d] > 1e300
  +         { cat("all.EQ(<mpfr>, <doubl>):", allEQ_0(m.dg[fin.d], d.dg[fin.d]), "\n")
  +           shp > 100  ||                   all.equal(m.dg[fin.d], d.dg[fin.d],
  +                                                     tol = 1e-13) # 2.063241e-14
  +         }
  +         ## compare with log scale :
  +         if(any(pos.d <- m.dg > 0)) {
  +             cat("all.EQ(log(d), d*(log)):",
  +               allEQ_0  (log(m.dg[pos.d]), m.ldg[pos.d]),"\n")
  +               all.equal(log(m.dg[pos.d]), m.ldg[pos.d], tol = 1e-14)
  +         }
  +     }, allow.logical0 = TRUE)
  + }
  shape = 2^ -20 :
  -------------
  all.EQ(<mpfr>, <doubl>): Mean relative difference: 5.766884e-16 
  all.EQ(log(d), d*(log)): Mean relative difference: 4.826256e-16 
  shape = 2^ -3 :
  -------------
  all.EQ(<mpfr>, <doubl>): Mean relative difference: 1.723515e-16 
  all.EQ(log(d), d*(log)): Mean relative difference: 8.342404e-20 
  shape = 2^ -1 :
  -------------
  all.EQ(<mpfr>, <doubl>): Mean relative difference: 2.536642e-17 
  all.EQ(log(d), d*(log)): Mean relative difference: 7.917217e-20 
  shape = 2^ 0 :
  -------------
  all.EQ(<mpfr>, <doubl>): Mean relative difference: 2.432573e-17 
  all.EQ(log(d), d*(log)): Mean relative difference: 1.509670e-19 
  shape = 2^ 1 :
  -------------
  all.EQ(<mpfr>, <doubl>): Mean relative difference: 1.133389e-16 
  all.EQ(log(d), d*(log)): TRUE 
  shape = 2^ 4 :
  -------------
  all.EQ(<mpfr>, <doubl>): Mean relative difference: 1.911332e-16 
  all.EQ(log(d), d*(log)): TRUE 
  shape = 2^ 10 :
  -------------
  all.EQ(<mpfr>, <doubl>): Mean relative difference: 5.728002e-17 
  all.EQ(log(d), d*(log)): TRUE 
  shape = 2^ 50 :
  -------------
  all.EQ(<mpfr>, <doubl>): Mean relative difference: 0.001523136 
  all.EQ(log(d), d*(log)): TRUE 
  > 
  > cat('Time elapsed: ', proc.time(),'\n') # "stats"
  Time elapsed:  8.87 0.17 9.06 NA NA 
  > if(!interactive()) warnings()
  Warning message:
  In matplot(x0 + 1e-07, relE2, type = "b", log = "x", main = "rel.err{  dt(x, df=df) }") :
    default 'pch' is smaller than number of columns and hence recycled
  > 
  > proc.time()
     user  system elapsed 
     8.87    0.17    9.06