R Under development (unstable) (2023-04-20 r84291 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 (64-bit) 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. > 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() 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.32 0.14 1.4 NA NA > if(!interactive()) warnings() > > proc.time() user system elapsed 1.32 0.14 1.40