R Under development (unstable) (2023-11-26 r85638 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. > library(gmp) Attaching package: 'gmp' The following objects are masked from 'package:base': %*%, apply, crossprod, matrix, tcrossprod > > ## From ~/R/Pkgs/Matrix/inst/test-tools-1.R -- only for R <= 3.0.1 -- > ##' @title Ensure evaluating 'expr' signals an error > ##' @param expr > ##' @return the caught error, invisibly > ##' @author Martin Maechler > assertError <- function(expr, verbose=getOption("verbose")) { + d.expr <- deparse(substitute(expr)) + t.res <- tryCatch(expr, error = function(e) e) + if(!inherits(t.res, "error")) + stop(d.expr, "\n\t did not give an error", call. = FALSE) + if(verbose) cat("Asserted Error:", conditionMessage(t.res),"\n") + invisible(t.res) + } > > Z1 <- as.bigz(1) ; Z1[FALSE] bigz(0) > Q1 <- as.bigq(1) ; Q1[FALSE] bigq(0) > stopifnot(0 == length(z0 <- as.bigz(0[FALSE])),# failed earlier + 0 == length(q0 <- as.bigq(0[FALSE])),# ditto + is.bigz(Z1), is.bigz(z0), !is.bigz(1L), !is.bigz(1), !is.bigz(Q1), + is.bigq(Q1), is.bigq(q0), !is.bigq(1L), !is.bigq(1/2), !is.bigq(Z1)) > > Z1[integer()] <- 2 # segfaulted earlier > Q1[integer()] <- 2 # ditto > assertError(Z1[1] <- list(1)) # segfaulted > assertError(Q1[1] <- list(1)) # " > assertError(Z1[1] <- NULL ) # " > assertError(Q1[1] <- NULL ) # " > > stopifnot(identical(Z1, as.bigz(1L)), identical(Q1, as.bigq(1L)), + identical(1L, as.integer(Z1)), + identical(1L, as.integer(Q1)),## failed earlier + identical(as.bigz(1[FALSE]), Z1[FALSE]), + identical(as.bigz(1[-1]), Z1[-1]), + identical(Z1[-1], rep(Z1, 0)) + , ##----------- bigq ------------- + identical(as.bigq(1[FALSE]), Q1[-1]), + identical(Q1[FALSE], Q1[-1]), + identical(Q1[-1], rep(Q1, 0)), + identical(q0, rep(Q1, 0)) + ) > > stopifnot(length(1[0]) == 0, 0 == length(Z1[0])) > Z <- as.bigz(I <- 2^(5*0:5)); mZ <- as.bigz(mI <- matrix(I, 2,3)) > Q <- Z / 4 ; mQ <- matrix(Q, 2,3) > > ii <- c(3:2,0:2,1:0,0:2) > i. <- c(2:0,1:0,1); j. <- ii[1:7] > i <- i.[i. != 0] > j <- j.[j. != 0] > I[ii] ; mI[i.,j.] [1] 1024 32 1 32 1 1 32 [,1] [,2] [,3] [,4] [,5] [1,] 33554432 32768 32 32768 32 [2,] 1048576 1024 1 1024 1 [3,] 1048576 1024 1 1024 1 [4,] 1048576 1024 1 1024 1 > stopifnot(all.equal( Z[ii], I[ii], tol=0), + all.equal(4*Q[ii], I[ii], tol=0), + identical(mI[i,j], mI[i.,j.]), + identical(mZ[i,j], mZ[i.,j.]), + identical(mQ[i,j], mQ[i.,j.])) > stopifnot(all.equal(asNumeric(mZ[i,j]), mI[i,j], tol=0), + all.equal( 4*mQ[i,j], mI[i,j], tol=0)) > > ## Outside indexing for *matrices* now gives an error: > assertError(mI[1,4]); assertError(mZ[1,4]); assertError(mQ[1,4]) > assertError(mI[3,2]); assertError(mZ[3,2]); assertError(mQ[3,2]) > ## whereas outside indexing of vectors should give NA: > stopifnot(identical(I[8:5], asNumeric(Z[8:5])), + identical(I[8:5], asNumeric(Q[8:5] * 4))) > > ## "basics", including as.matrix(), as.array(), as.list() : > i <- 1:9 > (x <- as.bigz(i, mod = 3)) [1] (1 %% 3) (2 %% 3) (0 %% 3) (1 %% 3) (2 %% 3) (0 %% 3) (1 %% 3) (2 %% 3) [9] (0 %% 3) > mx <- as.matrix(x) ## used to "bomb" badly: > ## (terminate called after throwing an instance of 'std::bad_alloc') > lx <- as.list(x) > stopifnot(5*x == (5*i) %% 3, + identical(as.bigz(x), x), # was not the case in gmp 0.5-14 + identical(mx, as.array(x)), + is(mx, "bigz"), dim(mx) == c(9,1), + is.list(lx), + identical(unlist(lx), + unlist(lapply(x, unclass)))) > > ## remove modulus "the new way" (NULL did fail): > modulus(x) <- NULL > Q <- x / 2 > mq <- as.matrix(Q) > lq <- as.list(Q) > stopifnot(identical(x, as.bigz(i %% 3)), + identical(mq, as.array(Q)), + is(mq, "bigq"), dim(mq) == c(9,1), + is.list(lq), + identical(unlist(lq), + unlist(lapply(Q, unclass)))) > > ## Check that as.bigq() is exact *and* asNumeric() is its inverse -------------- > set.seed(47) > summary(x1 <- rt(10000, df = 0.5)) # really long tailed Min. 1st Qu. Median Mean 3rd Qu. Max. -1.881e+11 -2.000e+00 0.000e+00 -1.880e+07 1.000e+00 9.875e+07 > summary(x2 <- rlnorm(10000, 200, 100)) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.000e+00 2.624e+57 7.518e+86 2.256e+250 2.840e+116 2.256e+254 > x <- c(x1, x2) > qx <- as.bigq(x) > nx <- asNumeric(qx) ## asNumeric()'s method for "bigq" is internal .bigq2num() > stopifnot(identical(x, nx), + identical(nx, gmp:::.bigq2num(qx)) + ) > > ## duplicated(), unique() : ---------------------- > q7 <- as.bigq(-5:7, 7) > if(FALSE)# not yet {well, *HARD* / impossible(?) without S4 } + Q <- q7^2 * as.bigz(77)^10 > Q <- q7^2 * as.bigq(77, 2)^10 > (uQ <- unique(Q)) Big Rational ('bigq') object of length 8: [1] 3738102281931735025/1024 149524091277269401/64 1345716821495424609/1024 [4] 149524091277269401/256 149524091277269401/1024 0 [7] 1345716821495424609/256 7326680472586200649/1024 > (sDup <- sum(duplicated(Q))) # = 5 [1] 5 > stopifnot(!duplicated(uQ), + sDup + length(uQ) == length(Q)) > nQ <- asNumeric(Q) > > stopifnot( identical(duplicated(Q), duplicated(nQ)) + , all.equal(unique(Q), unique(nQ)) + , sort(asNumeric(unique(denominator(Q)))) == 4^c(0, 3:5) + , TRUE) > > ## _ TODO _ rep() [times, length.out, each] > checkRep <- function(x) { + if((n <- length(x)) < 2) stop("'length(x)' must at least be 2, for these checks") + ii <- seq_len(n) + n1 <- pmin(.9*n, n-1) + stopifnot(identical(rep(x, 1), x), + identical(rep(x, 3), c(x,x,x)), + identical(rep(x, length.out=n1), x[1:n1]) + , + identical(rep(x, length.out=n+2), x[c(ii,1:2)]) + , ## times is *not* considered when 'length.out' is specified: + identical(rep(x, 4, length.out=n+2), x[c(ii,1:2)]) + , + identical(rep(x, 2, length.out=n1), x[1:n1]) + , + identical(x, rep(x, each=2)[2*ii]) + ) + } > > checkRep(Q) > checkRep(q7) > (Nu <- numerator(uQ)) Big Integer ('bigz') object of length 8: [1] 3738102281931735025 149524091277269401 1345716821495424609 [4] 149524091277269401 149524091277269401 0 [7] 1345716821495424609 7326680472586200649 > checkRep(Nu) > > ##------ Now check that base :: pmin() / pmax() works *in simple cases* for bigz > ##------ (because rep(., length.out) works: > ## {{MM: compare with ~/R/Pkgs/Rmpfr/tests/arith-ex.R }} > (x <- as.bigz(ix <- 2^(3* 0:7))) Big Integer ('bigz') object of length 8: [1] 1 8 64 512 4096 32768 262144 2097152 > (x9 <- pmin(x,9)) Big Integer ('bigz') object of length 8: [1] 1 8 9 9 9 9 9 9 > xp123 <- pmax(x, 123) > stopifnot(x9 == c(1,8, rep(9,6)), + xp123[1:3] == 123, + xp123[-(1:3)] > 123) > > chk.pmin <- function(x) { + message(deparse(sys.call()),": ") + x9 <- pmin(x, 9) + xp123 <- pmax(x, 123) + stopifnot( + identical(x, pmin(x, Inf)), + identical(x9, pmin(x, 23, Inf, 9)) + , identical(dim(x9), dim(x)) + , identical(dim(xp123), dim(x)) + ) + } > chk.pmin(x) chk.pmin(x): > mx <- matrix(x, nrow=3) # with correct warning Warning message: In matrix.bigz(x, nrow = 3) : data length [8] is not a sub-multiple or multiple of the number of rows [3] in matrix > chk.pmin(mx) chk.pmin(mx): > qq <- x / 47 > Mq <- matrix(qq, nrow=3) # with correct warning Warning message: In matrix.bigq(qq, nrow = 3) : data length [8] is not a sub-multiple or multiple of the number of rows [3] in matrix > if(FALSE) { ## FIXME: pmin() / pmax() are completely wrong for "bigq" !! + chk.pmin(qq) + chk.pmin(Mq) + } > > ## [<- : Used to return a *matrix* -- not what we want! > chk.subassign <- function(x, i, value) { + x0 <- x + x[i] <- value + stopifnot(identical(dim(x0), dim(x)), # only when not indexing *outside* + all(x[i] == value))# not always identical() + invisible(x) + } > > x. <- chk.subassign(x , 1, -1) > q. <- chk.subassign(qq, 1, -1) > q. <- chk.subassign(Mq, 1, -1) > x. <- chk.subassign(mx, 1, -1) > > if(require("Rmpfr") && packageVersion("Rmpfr") >= "0.5-2") { + stopifnot( + all.equal(pmin(14, x, 9), + pmin(14, ix, 9), tol=0) + , + all.equal(mq <- pmin(14, x/3, 9), ## numbers + bigq + pmin(14, ix/3, 9), tol= 1e-15) + , + is.bigq(mq)) + ## + ## Now, does pmin etc still work for bigz {it did fail!} + chk.pmin(x) + if(FALSE) ## FIXME: "Rmpfr's pmin / pmax methods destroy this ==> Fix Rmpfr! + chk.pmin(mx) + if(FALSE) { ## FIXME: pmin() / pmax() are completely wrong for "bigq" !! + chk.pmin(qq) + chk.pmin(Mq) + } + ## + ## Ditto for "[<-" : + x. <- chk.subassign(x , 1, -1) + q. <- chk.subassign(qq, 1, -1) + q. <- chk.subassign(Mq, 1, -1) + x. <- chk.subassign(mx, 1, -1) + ## + } else + message("{Rmpfr + gmp} checks __not__ done") Loading required package: Rmpfr 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 chk.pmin(x): > > ##--------------------------- order(), sort.list() -------------------------- > x <- as.bigz("0x123456789abcdef") # my secret message > B <- x + as.bigz(2)^(110:100) > (dB <- diff(B)) # now works Big Integer ('bigz') object of length 10: [1] -649037107316853453566312041152512 -324518553658426726783156020576256 [3] -162259276829213363391578010288128 -81129638414606681695789005144064 [5] -40564819207303340847894502572032 -20282409603651670423947251286016 [7] -10141204801825835211973625643008 -5070602400912917605986812821504 [9] -2535301200456458802993406410752 -1267650600228229401496703205376 > stopifnot(dB < 0, + log2(-dB) == 109:100 # 2^{n+1} - 2^n == 2^n + ) > rev(B) # is sorted Big Integer ('bigz') object of length 11: [1] 1267650600228311387025919692271 2535301200456540788522622897647 [3] 5070602400912999591516029308399 10141204801825917197502842129903 [5] 20282409603651752409476467772911 40564819207303422833423719058927 [7] 81129638414606763681318221630959 162259276829213445377107226775023 [9] 324518553658426808768685237063151 649037107316853535551841257639407 [11] 1298074214633706989118153298791919 > is.unsorted(rev(B))# TRUE but should be FALSE [1] TRUE > if(FALSE) ## not yet + identical(sort(B), rev(B)) > > proc.time() user system elapsed 0.59 0.07 0.67