R Under development (unstable) (2023-06-27 r84609 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 > > ## for reference (==> *not* using *.Rout.save here!) > sessionInfo() R Under development (unstable) (2023-06-27 r84609 ucrt) Platform: x86_64-w64-mingw32/x64 Running under: Windows Server 2022 x64 (build 20348) Matrix products: default locale: [1] LC_COLLATE=C LC_CTYPE=German_Germany.utf8 [3] LC_MONETARY=C LC_NUMERIC=C [5] LC_TIME=C time zone: Europe/Berlin tzcode source: internal attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] gmp_0.7-2 loaded via a namespace (and not attached): [1] compiler_4.4.0 > packageDescription("gmp") Package: gmp Version: 0.7-2 Date: 2023-06-28 Title: Multiple Precision Arithmetic Author: Antoine Lucas, Immanuel Scholz, Rainer Boehme , Sylvain Jasson , Martin Maechler Maintainer: Antoine Lucas Description: Multiple Precision Arithmetic (big integers and rationals, prime number tests, matrix computation), "arithmetic without limitations" using the C library GMP (GNU Multiple Precision Arithmetic). Depends: R (>= 3.5.0) Imports: methods Suggests: Rmpfr, MASS, round SystemRequirements: gmp (>= 4.2.3) License: GPL (>= 2) BuildResaveData: no LazyDataNote: not available, as we use data/*.R *and* our classes NeedsCompilation: yes URL: https://forgemia.inra.fr/sylvain.jasson/gmp Packaged: 2023-06-29 15:47:47 UTC; antoine Built: R 4.4.0; x86_64-w64-mingw32; 2023-06-29 18:36:46 UTC; windows Archs: x64 -- File: D:/RCompile/CRANincoming/R-devel/lib/gmp/Meta/package.rds > > ##' an (x == y) which gives TRUE also when both are NA: > isEQ <- function(x,y) (x == y) | (is.na(x) & is.na(y)) > > ## want to test all these > (ops <- sapply(getGroupMembers("Ops"), getGroupMembers)) $Arith [1] "+" "-" "*" "^" "%%" "%/%" "/" $Compare [1] "==" ">" "<" "!=" "<=" ">=" $Logic [1] "&" "|" > > N. <- as.bigz(NA) > Nq <- as.bigq(NA) > stopifnot(identical(Nq, as.bigq(N.)), + identical(N., as.bigz(Nq)))# used to fail > > xx <- c(NaN, NA, -Inf, -123:-121, -1:2, 7:8, Inf) > (xxI <- as.bigz(xx))# Inf's and NaN's do not exist ==> very large integers for +/- Inf Big Integer ('bigz') object of length 13: [1] NA [2] NA [3] 0 [4] -123 [5] -122 [6] -121 [7] -1 [8] 0 [9] 1 [10] 2 [11] 7 [12] 8 [13] 173766203193809456599982445949435627061939786100117250547173286503262376022458008465094333630120854338003194362163007597987225472483598640843335685441710193966274131338557192586399006789292714554767500194796127964596906605976605873665859580600161998556511368530960400907199253450604168622770350228527124626728538626805418833470107651091641919900725415994689920112219170907023561354484047025713734651608777544579846111001059482132180956689444108315785401642188044178788629853592228467331730519810763559577944882016286493908631503101121166109571682295769470379514531105239965209245314082665518579335511291525230373316486697786532335206274149240813489201828773854353041855598709390675430960381072270432383913542702130202430186637321862331068861776780211082856984506050024895394320139435868484643843368002496089956046419964019877586845530207748994394501505588146979082629871366088121763790555364513243984244004147636040219136443410377798011608722717131323621700159335786445601947601694025107888293017058178562647175461026384343438874861406516767158373279032321096262126551620255666605185789463207944391905756886829667520553014724372245300878786091700563444079107099009003380230356461989260377273986023281444076082783406824471703499844642915587790146384758051663547775336021829171033411043796977042190519657861762804226147480755555085278062866268677842432851421790544407006581148631979148571299417963950579210719961422405768071335213324842709316205032078384168750091017964584060285240107161561019930505687950233196051962261970932008838279760834318101044311710769457048672103958655016388894770892065267451228938951370237422841366052736174160431593023473217066764172949768821843606479073866252864377064398085101223216558344281956767163876579889759124956035672317578122141070933058555310274598884089982879647974020264495921703064439532898207943134374576254840272047075633856749514044298135927611328433323640657533550512376900773273703275329924651465759145114579174356770593439987135755889403613364529029604049868233807295134382284730745937309910703657676103447124097631074153287120040247837143656624045055614076111832245239612708339272798262887437416818440064925049838443370805645609424314780108030016683461562597569371539974003402697903023830108053034645133078208043917492087248958344081026378788915528519967248989338592027124423914083391771884524464968645052058218151010508471258285907685355807229880747677634789376 > (x <- c(NA, xx[is.finite(xx)])) [1] NA -123 -122 -121 -1 0 1 2 7 8 > xI <- as.bigz(x) > xQ <- as.bigq(xI) > stopifnot(identical(xI, as.bigz(xQ)), + identical(numerator(xQ), xI)) # numerator( ) > > stopifnot(isEQ(x, as.integer(x)), isEQ(x, xI), isEQ(x, xQ), + identical(xQ, as.bigq(x)), + identical(is.na(x), is.na(xI)), identical(is.na(x), is.na(xQ)), + identical(is.finite(x), is.finite(xI)), + identical(is.finite(x), is.finite(xQ)), + identical(is.infinite(x), is.infinite(xI)), + identical(is.infinite(x), is.infinite(xQ)), + ## The next 4 all failed till 2012-05-05: + isEQ(x, as.integer(xI)), + isEQ(x, as.integer(xQ)), + isEQ(x, as.numeric(xI)), + isEQ(x, as.numeric(xQ)), + TRUE) > > ## Finally (2020-06-06): mixed arithmetic works : > stopifnot(exprs = { + isEQ(xI - xQ, c(NA, rep(0, 9))) + isEQ(xI + xQ, 2*xI) + isEQ(xI * xQ, x^2) + all.equal(xQ^xI, x^x) + ## as do mixed comparisons + (xI == xQ)[-1] + !(xI < xQ)[-1] + !(xI > xQ)[-1] + (xI >= xQ)[-1] + }) > > ## double precision factorial() is exact up to n=22 > stopifnot(factorialZ(0:22) == factorial(0:22)) > > ## factorialZ() etc must also work when passed a bigz instead of an integer; > ## till Jan.2014, they silently produced nonsense. > N <- as.bigz(n <- 3:8) > stopifnot(identical(factorialZ(N), factorialZ(n)), factorialZ (n) == factorial(n), + identical(chooseZ(12, N), chooseZ(12, n)), chooseZ(12,n) == choose(12,n), + identical(fibnum (N), fibnum (n)), + identical(fibnum2(N), fibnum2(n)), + identical(lucnum (N), lucnum (n)), + identical(lucnum2(N), lucnum2(n))) > > > ## This one does *NOT* distinguish NA and NaN -- that's wanted here > EQ1 <- function(x,y) { + (abs(x-y) <= 1e-13*(abs(x)+abs(y)) & !(nx <- is.na(x)) & !(ny <- is.na(y))) | + (nx & ny) + } > stopifnot(EQ1(x, xI)) > EQ <- function(x,y) mapply(EQ1, x, y, USE.NAMES=FALSE) > > ## a version of outer() that should work also with these objects > mOuter <- function(X, Y=X, FUN, ...) { + lapply(seq_along(X), function(i) FUN(X[i], Y, ...)) + } > > matOuter <- function(X, Y=X, FUN, ...) { + t(array(unlist(mOuter(X, Y, FUN, ...)), + dim = c(length(Y), length(X)))) + } > > ##' @title > ##' @param OP an arithmetic OPerator such +, *,.. as R function > ##' @param u numeric vector > ##' @param uI a bigz/biginteger vector, "the same" as 'u'. > ##' @return a logical n x n matrix, say R, R[i,j] := TRUE iff > ##' u[i] OP v[j] are all the same when u,v vary in {u, uI}. > ##' @author Martin Maechler > opEQ <- function(OP, u, uI=as.bigz(u), eq=TRUE) { + stopifnot(length(u) == length(uI)) + if(eq) stopifnot(isEQ(u, uI)) # should be the case when result should be all TRUE + ## + ## choose only some on the RHS: + iR <- + if(no0.R <- (identical(OP, `/`) || identical(OP, `%/%`) || identical(OP, `%%`))) { + ## no zero on the RHS i.e., 2nd operand + is.na(u) | u != 0 + } else TRUE + ## choose only some on the LHS: + iL <- + if(no0.L <- (identical(OP, `^`))) { + ## no zero on the LHS i.e., 1st operand + is.na(u) | u != 0 + } else TRUE + ## + EQ(mOuter(u [iL],u [iR], OP) -> R, + mOuter(uI[iL],uI[iR], OP)) & + EQ(mOuter(u [iL],uI[iR], OP) -> S, + mOuter(uI[iL], u[iR], OP)) & + EQ(R, S) + } > > ## "Compare" - works "out of the box > eqC <- lapply(sapply(ops$Compare, get), + function(op) opEQ(op, x, xI)) > stopifnot(do.call(all, eqC)) > > opsA <- ops$Arith > > eqA <- lapply(sapply(opsA, get), function(op) opEQ(op, x, xI)) > > op6 <- c("+","-", "*", "/", "%/%", "^")## << are fine - now including "^" _and_ %/% ! > stopifnot(sapply(eqA, all)[op6]) > ## The others: now (2014-07): only %% is left: has several "wrong": > lapply(eqA[is.na(match(names(eqA), op6))], symnum) $`%%` [1,] | | | | | | | | | | [2,] | | . . . | . . . . [3,] | . | . . | . . . . [4,] | . . | . | . . . . [5,] | | | | | | | | | | [6,] | | | | | | | | | | [7,] | | | | | | | | | | [8,] | | | | | | | | | | [9,] | | | | | | | | | | > > ## For example: > symnum(opEQ(`%%`, x, xI))# not all TRUE, since, e.g., [1,] | | | | | | | | | | [2,] | | . . . | . . . . [3,] | . | . . | . . . . [4,] | . . | . | . . . . [5,] | | | | | | | | | | [6,] | | | | | | | | | | [7,] | | | | | | | | | | [8,] | | | | | | | | | | [9,] | | | | | | | | | | > x [3] %% x [1] NA -122 0 -1 0 NaN 0 0 4 6 > x [3] %% xI ## (negative turned into >= 0; warning 'division by zero') Big Integer ('bigz') object of length 10: [1] NA 1 0 120 0 NA 0 0 4 6 Warning message: In `%%.bigz`(x[3], xI) : biginteger division by zero: returning NA > > x %% x [3] [1] NA -1 0 -121 -1 0 -121 -120 -115 -114 > xI %% x [3] ## (no negatives ..) Big Integer ('bigz') object of length 10: [1] NA 121 0 1 121 0 1 2 7 8 > > > ##-- "^" ------------ > z1i <- 0:1 > z1n <- as.double(z1i) > c(NA^0, NA^0L, z1i^NA, z1n^NA)# <- in R (<= 2011), the first and last are 1 [1] 1 1 NA 1 NA 1 > stopifnot(isEQ(c(N.^0, N.^0L, z1i^N.), c(1,1,NA,1)), + isEQ(c(Nq^0, Nq^0L, z1i^Nq), c(1,1,NA,1))) > > ## need non-negative values: > x.po0 <- x >= 0 > stopifnot(M.pow <- opEQ(`^`, x[x.po0], xI[x.po0])) > if(FALSE)# FIXME + stopifnot(M.powQ <- opEQ(`^`, x[x.po0], xQ[x.po0])) > if(FALSE)# FIXME {z - q} + M.poIQ <- opEQ(`^`,xI[x.po0], xQ[x.po0]) > > ## Modulo arithmetic > i <- as.bigz(-5:10, 16); i <- i[i != 0]; i [1] (11 %% 16) (12 %% 16) (13 %% 16) (14 %% 16) (15 %% 16) (1 %% 16) [7] (2 %% 16) (3 %% 16) (4 %% 16) (5 %% 16) (6 %% 16) (7 %% 16) [13] (8 %% 16) (9 %% 16) (10 %% 16) > stopifnot(identical(as.integer(i), c(11:15, 1:10))) > (Ii <- 1/i )## BUG: in all versions of gmp up to 0.5-5 -- now 7 warnings pow(x, -|n|) [1] (3 %% 16) NA (5 %% 16) NA (15 %% 16) (1 %% 16) [7] NA (11 %% 16) NA (13 %% 16) NA (7 %% 16) [13] NA (9 %% 16) NA Warning messages: 1: In `/.bigz`(1, i) : pow(x, -|n|) returning NA as x has no inverse wrt modulus 2: In `/.bigz`(1, i) : pow(x, -|n|) returning NA as x has no inverse wrt modulus 3: In `/.bigz`(1, i) : pow(x, -|n|) returning NA as x has no inverse wrt modulus 4: In `/.bigz`(1, i) : pow(x, -|n|) returning NA as x has no inverse wrt modulus 5: In `/.bigz`(1, i) : pow(x, -|n|) returning NA as x has no inverse wrt modulus 6: In `/.bigz`(1, i) : pow(x, -|n|) returning NA as x has no inverse wrt modulus 7: In `/.bigz`(1, i) : pow(x, -|n|) returning NA as x has no inverse wrt modulus > I2 <- i^(-1)## BUG: not considering (mod) // segmentation fault in gmp 0.5-1 {now: 7 warn..} Warning messages: 1: In `^.bigz`(i, (-1)) : pow(x, -|n|) returning NA as x has no inverse wrt modulus 2: In `^.bigz`(i, (-1)) : pow(x, -|n|) returning NA as x has no inverse wrt modulus 3: In `^.bigz`(i, (-1)) : pow(x, -|n|) returning NA as x has no inverse wrt modulus 4: In `^.bigz`(i, (-1)) : pow(x, -|n|) returning NA as x has no inverse wrt modulus 5: In `^.bigz`(i, (-1)) : pow(x, -|n|) returning NA as x has no inverse wrt modulus 6: In `^.bigz`(i, (-1)) : pow(x, -|n|) returning NA as x has no inverse wrt modulus 7: In `^.bigz`(i, (-1)) : pow(x, -|n|) returning NA as x has no inverse wrt modulus > stopifnot(identical(Ii, I2), + is.na(Ii[c(2, 4, 7, 9, 11, 13, 15)]), + identical(Ii[c(1,3)], as.bigz(c(3,5), 16))) > (Iz <- 1/(z <- as.bigz(1:12, 13))) [1] (1 %% 13) (7 %% 13) (9 %% 13) (10 %% 13) (8 %% 13) (11 %% 13) [7] (2 %% 13) (5 %% 13) (3 %% 13) (4 %% 13) (6 %% 13) (12 %% 13) > stopifnot(identical(Iz, z^-1), + Iz == c(1, 7, 9, 10, 8, 11, 2, 5, 3, 4, 6, 12), + identical(modulus(Iz), as.bigz(13))) > ## The first two of course give fractions: > (r1 <- as.bigz(3) / 1:12) Big Rational ('bigq') object of length 12: [1] 3 3/2 1 3/4 3/5 1/2 3/7 3/8 1/3 3/10 3/11 1/4 > r2 <- as.bigz(3) / as.bigz(1:12) > stopifnot(identical(r1, r2)) > > ## Now, the new scheme : > (iLR <- as.bigz(3, 13) / as.bigz(1:12, 13)) [1] (3 %% 13) (8 %% 13) (1 %% 13) (4 %% 13) (11 %% 13) (7 %% 13) [7] (6 %% 13) (2 %% 13) (9 %% 13) (12 %% 13) (5 %% 13) (10 %% 13) > ## [1] (3 %% 13) (8 %% 13) (1 %% 13) (4 %% 13) (11 %% 13) (7 %% 13) > ## [7] (6 %% 13) (2 %% 13) (9 %% 13) (12 %% 13) (5 %% 13) (10 %% 13) > iL <- as.bigz(3, 13) / as.bigz(1:12) > iLi <- as.bigz(3, 13) / 1:12 > iR <- as.bigz(3) / as.bigz(1:12, 13) > iiR <- 3 / as.bigz(1:12, 13) > stopifnot(identical(iL, iLi) + , identical(iR, iiR) + , identical(iR, iLR) + , identical(iL, iR)) ## failed until recently... > > ## whereas these two always use divq.bigz : > (q <- as.bigz(3, 13) %/% as.bigz(1:12)) [1] (3 %% 13) (1 %% 13) (1 %% 13) (0 %% 13) (0 %% 13) (0 %% 13) (0 %% 13) [8] (0 %% 13) (0 %% 13) (0 %% 13) (0 %% 13) (0 %% 13) > ## [1] (3 %% 13) (1 %% 13) (1 %% 13) (0 %% 13) (0 %% 13) (0 %% 13) > ## [7] (0 %% 13) (0 %% 13) (0 %% 13) (0 %% 13) (0 %% 13) (0 %% 13) > stopifnot(identical(q, divq.bigz(as.bigz(3, 13), 1:12)), + ## --------- + identical(q, 3 %/% as.bigz(1:12, 13)), + q == c(3, 1, 1, rep(0,9))) > s <- as.bigz(3, 13) / as.bigz(1:12, 17) > ## used to give > ## Big Integer ('bigz') object of length 12: > ## [1] 3 1 1 0 0 0 0 0 0 0 0 0 > ## but now, really just `` drops the contradicting "mod" '' ==> uses rational: > stopifnot(identical(s, r1)) > > ##----- Z^e (modulo m) --------------- > z12 <- as.bigz(1:12,12) > stopifnot(identical(z12^1, z12), z12^0 == 1, + identical(z12^2, as.bigz(rep(c(1,4,9,4,1,0), 2), 12)), + identical(z12^3, + as.bigz(c(1,8,3:5,0,7:9,4,11,0), 12)), + identical(z12^4, z12^2), + identical(z12^5, z12^3), + identical(z12^6, z12^2), + identical(z12^6, (1:12) ^ as.bigz(6, 12)) + ) > > for(E in 6:20) { + ir <- as.integer(r <- z12 ^ E) + stopifnot(identical(modulus(r), as.bigz(12)), + 0 <= ir, ir <= 11) + } > > z17 <- as.bigz(1:16, 17) > stopifnot(z17^0 == 1, identical(z17^1, z17), identical(z17^-1, iz <- 1/z17), + identical(z17^-2, iz^2), (iz^2) * (sq <- z17^2) == 1, + modulus(sq) == 17, unique(sq) == (1:8)^2 %% 17) > > > > ##--- Log()s ------------------------- > (ex <- c(outer(c(2,5,10), 10^(1:3))))# 20 .. 10'000 [1] 20 50 100 200 500 1000 2000 5000 10000 > stopifnot(dim(L <- outer(as.bigz(2:4), ex, `^`)) == c(3, length(ex))) > l2 <- array(log2(L), dim = dim(L)) > lnL <- log(L) > a.EQ <- function(x,y, tol=1e-15, ...) all.equal(x,y, tol=tol, ...) > stopifnot(a.EQ(l2[1,], ex), + a.EQ(l2[3,], 2*ex), + a.EQ(log(L, 8), lnL/log(8)), + a.EQ(c(l2), lnL/log(2))) > > > ###------------------ bigq -------------------------------- > > > xQ1 <- as.bigq(x, 1) > eqC <- lapply(sapply(ops$Compare, get), function(op) opEQ(op, x, xQ1)) > stopifnot(Reduce(`&`, eqC))## > ## > xQ <- as.bigq(x, 17) # == x/17 .. *are* not equal, i.e., not expecting all TRUE: > eqQ <- lapply(sapply(ops$Compare, get), + function(op) opEQ(op, x, xQ, eq=FALSE)) > lapply(eqQ, symnum)## <- symnum, for nice output $`==` [1,] | | | | | | | | | | [2,] | . | | | | | | | | [3,] | | . | | | | | | | [4,] | | | . | | | | | | [5,] | | | | . | | | | | [6,] | | | | | | | | | | [7,] | | | | | | . | | | [8,] | | | | | | | . | | [9,] | | | | | | | | . | [10,] | | | | | | | | | . $`>` [1,] | | | | | | | | | | [2,] | . . . | | | | | | [3,] | . . . | | | | | | [4,] | . . . | | | | | | [5,] | | | | . | | | | | [6,] | | | | | | | | | | [7,] | | | | | | . . . . [8,] | | | | | | . . . . [9,] | | | | | | . . . . [10,] | | | | | | . . . . $`<` [1,] | | | | | | | | | | [2,] | . . . | | | | | | [3,] | . . . | | | | | | [4,] | . . . | | | | | | [5,] | | | | . | | | | | [6,] | | | | | | | | | | [7,] | | | | | | . . . . [8,] | | | | | | . . . . [9,] | | | | | | . . . . [10,] | | | | | | . . . . $`!=` [1,] | | | | | | | | | | [2,] | . | | | | | | | | [3,] | | . | | | | | | | [4,] | | | . | | | | | | [5,] | | | | . | | | | | [6,] | | | | | | | | | | [7,] | | | | | | . | | | [8,] | | | | | | | . | | [9,] | | | | | | | | . | [10,] | | | | | | | | | . $`<=` [1,] | | | | | | | | | | [2,] | . . . | | | | | | [3,] | . . . | | | | | | [4,] | . . . | | | | | | [5,] | | | | . | | | | | [6,] | | | | | | | | | | [7,] | | | | | | . . . . [8,] | | | | | | . . . . [9,] | | | | | | . . . . [10,] | | | | | | . . . . $`>=` [1,] | | | | | | | | | | [2,] | . . . | | | | | | [3,] | . . . | | | | | | [4,] | . . . | | | | | | [5,] | | | | . | | | | | [6,] | | | | | | | | | | [7,] | | | | | | . . . . [8,] | | | | | | . . . . [9,] | | | | | | . . . . [10,] | | | | | | . . . . > > Fn <- gmp:::pow.bigq; q <- 2.3 > stopifnot(inherits(e1 <- tryCatch(Fn(q,q), error=identity), "error"), + inherits(e2 <- tryCatch(q ^ as.bigq(1,3), error=identity), "error"), + grepl("Rmpfr", e1$message), + identical(e1$message, e2$message)) > > > ## FIXME(2): %% and %/% do not work at all for bigq > (opsA4 <- opsA[opsA != "^" & !grepl("^%", opsA)]) [1] "+" "-" "*" "/" > eqA1 <- lapply(sapply(opsA4, get), function(op) opEQ(op, x, xQ1)) > sapply(eqA1, table) +.TRUE -.TRUE *.TRUE /.TRUE 100 100 100 90 > ## .TRUE -.TRUE *.TRUE /.TRUE > ## 100 100 100 90 > ## ^^^^ (90: was 81) [not dividing by 0] > > ## xQ *is* different from x (apart from x[6] (and, NA x[1])) > eqA <- lapply(sapply(opsA4, get), function(op) opEQ(op, x, xQ, eq=FALSE)) > lapply(eqA, symnum) $`+` [1,] | | | | | | | | | | [2,] | . . . . . . . . . [3,] | . . . . . . . . . [4,] | . . . . . . . . . [5,] | . . . . . . . . . [6,] | . . . . | . . . . [7,] | . . . . . . . . . [8,] | . . . . . . . . . [9,] | . . . . . . . . . [10,] | . . . . . . . . . $`-` [1,] | | | | | | | | | | [2,] | . . . . . . . . . [3,] | . . . . . . . . . [4,] | . . . . . . . . . [5,] | . . . . . . . . . [6,] | . . . . | . . . . [7,] | . . . . . . . . . [8,] | . . . . . . . . . [9,] | . . . . . . . . . [10,] | . . . . . . . . . $`*` [1,] | | | | | | | | | | [2,] | . . . . | . . . . [3,] | . . . . | . . . . [4,] | . . . . | . . . . [5,] | . . . . | . . . . [6,] | | | | | | | | | | [7,] | . . . . | . . . . [8,] | . . . . | . . . . [9,] | . . . . | . . . . [10,] | . . . . | . . . . $`/` [1,] | | | | | | | | | | [2,] | . . . . | . . . . [3,] | . . . . | . . . . [4,] | . . . . | . . . . [5,] | . . . . | . . . . [6,] | . . . . | . . . . [7,] | . . . . | . . . . [8,] | . . . . | . . . . [9,] | . . . . | . . . . > > > ## round(x, digits) -- should work *and* be vectorized in both (x, digits) > x1 <- as.bigq((-19:19), 10) > stopifnot(round(x1, 1) == x1) > > half <- as.bigq(1, 2) > i1 <- (-19:29) > x <- half + i1 > cbind(x, round(x)) Big Rational ('bigq') 49 x 2 matrix: [,1] [,2] [1,] -37/2 -18 [2,] -35/2 -18 [3,] -33/2 -16 [4,] -31/2 -16 [5,] -29/2 -14 [6,] -27/2 -14 [7,] -25/2 -12 [8,] -23/2 -12 [9,] -21/2 -10 [10,] -19/2 -10 [11,] -17/2 -8 [12,] -15/2 -8 [13,] -13/2 -6 [14,] -11/2 -6 [15,] -9/2 -4 [16,] -7/2 -4 [17,] -5/2 -2 [18,] -3/2 -2 [19,] -1/2 0 [20,] 1/2 0 [21,] 3/2 2 [22,] 5/2 2 [23,] 7/2 4 [24,] 9/2 4 [25,] 11/2 6 [26,] 13/2 6 [27,] 15/2 8 [28,] 17/2 8 [29,] 19/2 10 [30,] 21/2 10 [31,] 23/2 12 [32,] 25/2 12 [33,] 27/2 14 [34,] 29/2 14 [35,] 31/2 16 [36,] 33/2 16 [37,] 35/2 18 [38,] 37/2 18 [39,] 39/2 20 [40,] 41/2 20 [41,] 43/2 22 [42,] 45/2 22 [43,] 47/2 24 [44,] 49/2 24 [45,] 51/2 26 [46,] 53/2 26 [47,] 55/2 28 [48,] 57/2 28 [49,] 59/2 30 > rx1 <- round(x/10, 1) > stopifnot(exprs = { + as.bigz(round(x)) %% 2 == 0 + identical(round(x) > x, i1 %% 2 == 1) + (rx1 - x/10) * 20 == c(1,-1) # {recycling up/down}: perfect rounding to even + (round(x/100, 2) - x/100) * 200 == c(1,-1) # (ditto) + }) > (drx1 <- asNumeric(rx1))# shows perfect round to *even* [1] -1.8 -1.8 -1.6 -1.6 -1.4 -1.4 -1.2 -1.2 -1.0 -1.0 -0.8 -0.8 -0.6 -0.6 -0.4 [16] -0.4 -0.2 -0.2 0.0 0.0 0.2 0.2 0.4 0.4 0.6 0.6 0.8 0.8 1.0 1.0 [31] 1.2 1.2 1.4 1.4 1.6 1.6 1.8 1.8 2.0 2.0 2.2 2.2 2.4 2.4 2.6 [46] 2.6 2.8 2.8 3.0 > ## but double precision rounding cannot be perfect (as numbers are not exact!): > dx <- asNumeric(x/10) > dx1 <- round(dx, 1) > dmat <- cbind(x=dx, r.x = dx1, rQx = drx1) > ## shows "the picture" a bit {see Martin's vignette in CRAN package 'round'}: > noquote(cbind(apply(dmat, 2, formatC), + ER = ifelse(abs(dx1 - drx1) > 1e-10, "*", ""))) x r.x rQx ER [1,] -1.85 -1.8 -1.8 [2,] -1.75 -1.8 -1.8 [3,] -1.65 -1.6 -1.6 [4,] -1.55 -1.5 -1.6 * [5,] -1.45 -1.4 -1.4 [6,] -1.35 -1.3 -1.4 * [7,] -1.25 -1.2 -1.2 [8,] -1.15 -1.1 -1.2 * [9,] -1.05 -1 -1 [10,] -0.95 -0.9 -1 * [11,] -0.85 -0.8 -0.8 [12,] -0.75 -0.8 -0.8 [13,] -0.65 -0.6 -0.6 [14,] -0.55 -0.5 -0.6 * [15,] -0.45 -0.4 -0.4 [16,] -0.35 -0.3 -0.4 * [17,] -0.25 -0.2 -0.2 [18,] -0.15 -0.1 -0.2 * [19,] -0.05 -0 0 [20,] 0.05 0 0 [21,] 0.15 0.1 0.2 * [22,] 0.25 0.2 0.2 [23,] 0.35 0.3 0.4 * [24,] 0.45 0.4 0.4 [25,] 0.55 0.5 0.6 * [26,] 0.65 0.6 0.6 [27,] 0.75 0.8 0.8 [28,] 0.85 0.8 0.8 [29,] 0.95 0.9 1 * [30,] 1.05 1 1 [31,] 1.15 1.1 1.2 * [32,] 1.25 1.2 1.2 [33,] 1.35 1.3 1.4 * [34,] 1.45 1.4 1.4 [35,] 1.55 1.5 1.6 * [36,] 1.65 1.6 1.6 [37,] 1.75 1.8 1.8 [38,] 1.85 1.8 1.8 [39,] 1.95 2 2 [40,] 2.05 2 2 [41,] 2.15 2.1 2.2 * [42,] 2.25 2.2 2.2 [43,] 2.35 2.3 2.4 * [44,] 2.45 2.4 2.4 [45,] 2.55 2.5 2.6 * [46,] 2.65 2.6 2.6 [47,] 2.75 2.8 2.8 [48,] 2.85 2.8 2.8 [49,] 2.95 2.9 3 * > > ## standard R: > rd <- round(pi*10^(-2:5), digits=7:0) > formatC(rd, digits=12, width=1) [1] "0.0314159" "0.314159" "3.14159" "31.4159" "314.159" "3141.59" [7] "31415.9" "314159" > ## bigq -- show we vectorize in both x, digits > (rQ <- round(as.bigq(pi*10^(-2:5)), digits=7:0)) Big Rational ('bigq') object of length 8: [1] 314159/10000000 314159/1000000 314159/100000 314159/10000 [5] 314159/1000 314159/100 314159/10 314159 > stopifnot(exprs = { + as.integer(numerator (rQ)) == 314159L + as.integer(denominator(rQ)) == 10^(7:0) + all.equal(asNumeric(rQ), rd, tol = 1e-15) + }) > > > > proc.time() user system elapsed 0.43 0.07 0.50