Nothing
require("Rmpfr")
## 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)
## 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())
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)
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))
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)))
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))
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 ): ")
## FIXME checkPmin(x, qx); cat("[Ok]\n")
##
print( base::pmin(Z, Z, max(Z)) )# via gmp:: rep.bigz(x, length.out = *)
cat("checking pmin(. bigz )
[currently with lots of pmin() and pmax(...) warnings 'incompatible methods]:\n ")
checkPmin(Z); cat("[Ok]\n") # via gmp:: all.equal.bigz()
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)
(todays.seed <- eval(parse(text=Sys.Date())))# so this is reproducible
# (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))
}
}
}
}
## 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"
## 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"
(r <- dn / (dn + dOO)) # 0.915555 (double prec arithmetic)
## 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.....
stopifnot(exprs = {
all.equal(n + dOO, dn + dOO)
all.equal(n / (n + dOO), r)
})
## log(., base) :
(ten40 <- as.bigz(10)^40)
ten40m <- mpfr(ten40)
(lt40 <- log(ten40m, 10)) # gave Error in ... : base != exp(1) is not yet implemented
## '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)
})
###------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"
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.