R Under development (unstable) (2023-12-02 r85657 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 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. > 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] 19695 > # (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 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.87 0.09 4.96 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.96 0.09 5.06 NA NA > > proc.time() user system elapsed 4.96 0.09 5.06