#### Low level stuff - debugging etc #### ========= ========= require("Rmpfr") options(warn = 2)# warning -> error identical3 <- function(x,y,z) identical(x,y) && identical (y,z) identical4 <- function(a,b,c,d) identical(a,b) && identical3(b,c,d) ## sane state [when re-source()ing this file]: .mpfr_erange_set("Emin", -(2^30-1)) .mpfr_erange_set("Emax", +(2^30-1)) ###----- _1_ mpfr1 , import, xport etc ----------------------------------------- i8 <- mpfr(-2:5, 32) x4 <- mpfr(c(NA, NaN, -Inf, Inf), 32); x4 # NA -> NaN as well stopifnot(identical3(is.na(x4), is.nan(x4), c(T,T,F,F))) o1 <- as(x4[1], "mpfr1") stopifnot(is(o1, "mpfr1")) # failed previously validObject(o1) # ditto (failed on 64-bit only) stopifnot( getPrec("0xabc", base=16, doNumeric=FALSE) == 3*4, getPrec( "abc", base=16, doNumeric=FALSE) == 3*4, getPrec("0b1001", base=2, doNumeric=FALSE) == 4, getPrec( "1001", base=2, doNumeric=FALSE) == 4, identical3(mpfr("0b101", base= 2), mpfr( "101", base= 2), mpfr(5, precBits = 3)) , identical3(mpfr("0xabc", base=16), mpfr( "abc", base=16), mpfr(2748, base=16, precBits = 12)) ) ## save initial (Emin, Emax) eranges : erangesOrig <- .mpfr_erange() ###----- _2_ Debugging, changing MPFR defaults, .. ----------------------------- ## NB: Currently mostly *not* documented, not even .mpfr_erange() stopifnot(Rmpfr:::.mpfr_debug() == 0 # the default level ## Activate debugging level 1: , Rmpfr:::.mpfr_debug(1) == 0 # the previous level ## and check it : , Rmpfr:::.mpfr_debug() == 1 # the current level ) r <- mpfr(7, 100)^-1000 r ## (same as without debugging) ## where as this does print info: -- notably the very large values [3..6]: .eranges <- function() sapply(.mpfr_erange_kinds, .mpfr_erange, USE.NAMES=FALSE) ## now, mpfr_erange() works with a *vector* of args: .erange2 <- function() .mpfr_erange(.mpfr_erange_kinds) ## now returning *double* - which loses some precision [ending in '04' instead of '03']: formatC(.eranges(), format="fg") stopifnot(identical(.eranges(), .erange2())) .mpfr_minPrec() .mpfr_maxPrec()# debug printing shows the long integer (on 64 bit) ## Now, level 2 : stopifnot(Rmpfr:::.mpfr_debug(2) == 1) r ## with quite a bit of output if(FALSE) # on Winbuilder [2019-08-08, both 32 and 64 bit]: .mpfr_erange_set("Emax", 1073741823) r2 <- r^100 r2 L <- r^-100000 L3 <- L^3 str(L3, internal=TRUE) ## Class 'mpfr' [package "Rmpfr"] of length 1 and precision 100 ## internally @.Data: List of 1 ## $ :Formal class 'mpfr1' [package "Rmpfr"] with 4 slots ## .. ..@ prec: int 100 ## .. ..@ exp : int [1:2] 842206477 0 ## .. ..@ sign: int 1 ## .. ..@ d : int [1:4] 268435456 761715680 1492345294 -1000766770 str(L3) ## lots of debugging output, then ## 1.00989692356e+253529412 ## ^^~~~~~~~~~ 10 ^ 253'529'412 that is humongous if(!interactive()) # not seg.faulting, but printing a *huge* line [no longer!] show(L3) ## segmentation fault -- randomly; 2017-06: no longer see any problem, not even with if(FALSE) ## well, not really, definitely not interactively for now if(interactive()) for(i in 1:256) show(L3) ## ## quite platform dependent {valgrind ==> bug? even in mpfr/gmp/.. ?} str(.mpfr2list(x4)) ## slightly nicer ["uniformly not worse"] (still very similar) : str(x4, internal=TRUE) x4 ## "similar info" as .mpfr2list(.) ## Increase maximal exponent: tools:::assertWarning( .mpfr_erange_set("Emax", 5e18)) # too large {FIXME why only warning and not error ??} .mpfr_erange("Emax") # is unchanged if(4e18 < .mpfr_erange("max.emax")) { .mpfr_erange_set("Emax", 4e18) # now ok: stopifnot(.mpfr_erange("Emax") == 4e18) } ## revert to no debugging: stopifnot(Rmpfr:::.mpfr_debug(0) == 2) .mpfr_maxPrec() L / (r2^-1000)# 1.00000....448 (could be more accurate?) stopifnot(exprs = { all.equal(L, r2^-1000, tol= 1e-27) # why not more accurate? all.equal(log(L), -100000 * (-1000) * log(7), tol = 1e-15) }) ## Now, our experimental "transport vehicle": stopifnot(length(rv <- c(r, r2, L)) == 3) str(mpfrXport(rv)) str(mpfrXport(mpfr(2, 64)^(-3:3))) str(mpfrXport(Const("pi")* 2^(-3:3))) ## and a very large one mil <- mpfr(1025, 111) str(mm <- mpfrXport(xx <- mil^(2^25))) stopifnot(all.equal(log2(xx) * 2^-25, log2(mil), tol=1e-15)) ## even larger -- strictly needs extended erange: if(.mpfr_erange("min.emin") <= -2^40) { .mpfr_erange_set("Emin", - 2^40) show(xe <- 2^mpfr(-seq(1,70, by=3)*8e8, 64)) ## used to print wrongly {because of integer overflow in .mpfr2str()$exp}, ## with some exponents large positive stopifnot(exprs = { ! .mpfr_erange_is_int() # as 'exp's now are double (ee <- as.numeric(sub(".*e","", formatMpfr(xe)))) < -240e6 (diff(ee) + 722471990) %in% 0:1 }) } else { cat(sprintf( "Cannot set 'Emin' to -2^40 (= %g), as .mpfr_erange(\"min.emin\") is larger, namely %g.\n", - 2^40, .mpfr_erange("min.emin"))) } ## Bill Dunlap's example (with patch about convert S_alloc bug): ## (precision increases, then decreases) z <- c(mpfr(1,8)/19, mpfr(1,32)/19, mpfr(1,24)/19) cbind(fz <- format(z)) stopifnot(identical(fz, rev(format(rev(z))))) stopifnot(identical(fz, c("0.05273", "0.052631578947", "0.0526315793"))) # << smaller prec, again since 2019-08-09 e.xx. <- .mpfr2exp(xx) e.z. <- .mpfr2exp(z) ## revert to original 'erange' settings (which gives integer 'exp'): .mpfr_erange_set("Emax", erangesOrig[["Emax"]]) # typically 2^30 - 1 = 1073741823 .mpfr_erange_set("Emin", erangesOrig[["Emin"]]) e.xx <- .mpfr2exp(xx) e.z <- .mpfr2exp(z) stopifnot(exprs = { .mpfr_erange_is_int() e.xx == e.xx. e.xx == 335591572 e.z == e.z. e.z == -4 is.integer(e.xx) # but e.xx. is double is.integer(e.z) }) k1 <- mpfr( c(123, 1234, 12345, 123456), precBits=2) (N1 <- asNumeric(k1))# 128 1024 12288 131072 -- correct str(sk1 <- .mpfr2str(k1)) str(sk1. <- .mpfr2str(k1, maybe.full=TRUE)) str(sk1.2 <- .mpfr2str(k1, digits=2, base=2)) str(sk1.2F <- .mpfr2str(k1, maybe.full=TRUE, base=2)) stopifnot(exprs = { identical(sk1 [1:2], list(str = c("13", "10", "12", "13"), exp = 3:6)) identical(sk1.[1:2], list(str = c("128", "1024", "12288", "131072"), exp = 3:6)) identical(sk1.2, list(str = c("10", "10", "11", "10"), exp = c( 8L, 11L, 14L, 18L), finite = rep(TRUE, 4), is.0 = rep(FALSE, 4))) all.equal(sk1.2[2:4], .mpfr_formatinfo(k1), tol=0) # not identical(): int <-> double identical(formatMpfr(k1, base=2, digits=20, drop0trailing=TRUE), with(sk1.2, paste0(str, sapply(exp - nchar(str), strrep, x="0")))) identical(formatMpfr(k1, base=2, digits=2, exponent.plus=FALSE), c("1.0e7", "1.0e10", "1.1e13", "1.0e17")) }) ## MM: --> need_dig is fine but is not used in the string that is returned !! (fk1sF <- formatMpfr(k1, scientific=FALSE)) # "the bug" --- now fixed! ==> new "Bug" in new Rmpfr ???? ## was "128." "1024." "12288." "131072." , but now obeying internal precision gives ## "1.e2" "1.e3" "1.e4" "1.e5" (fk1 <- formatMpfr(k1, digits=6)) stopifnot(exprs = { N1 == as.numeric(fk1) ## FIXME: This should change again "1024" identical(format(k1, digits=3), c("128.", "1020.", "1.23e+4", "1.31e+5")) }) ## digs <- setNames(1:6, 1:6) ## Each of these are 4 x 6 matrices ffix <- sapply(digs, function(d) format(k1, digits = d, scientific = FALSE)) ## *not* good at all .. ## ==> need a maybe.full=TRUE even here ff <- sapply(digs, function(d) format(k1, digits = d))# sci..fic = NA -- digits=1 failing for '128' fsci <- sapply(digs, function(d) format(k1, digits = d, scientific = TRUE)) # perfect stopifnot(exprs = { length(dd <- dim(ff)) == 2 identical(dd, dim(ffix)) identical(dd, dim(fsci)) all.equal(asNumeric(fsci), asNumeric(ffix) -> dmat, tol=0) all.equal(asNumeric(ff), asNumeric(ffix), tol=0) }) rE <- 1 - dmat / asNumeric(k1) i <- 1:5 summary(fm <- lm(log10(colMeans(abs(rE)))[i] ~ i)) stopifnot(exprs = { rE[ cbind(FALSE, upper.tri(rE)[,-6]) ] == 0 abs(residuals(fm)) < 0.15 }) ## formatting / printing : tenth <- mpfr(-12:12, 52)/10 cents <- mpfr(-11:11, 64)/100 (kxi <- sort(c(k1, x4, i8, tenth, cents), na.last=FALSE)) mstr <- .mpfr2str (kxi) mfi <- .mpfr_formatinfo(kxi) es <- mstr$exp # base 10 ; with '0' when !is.finite or is0 ef <- mfi $exp # base 2 ; "undefined" when !is.finite or is0 j2 <- c("finite", "is.0") dxi <- cbind(x = asNumeric(kxi), prec = .getPrec(kxi), as.data.frame(mstr, stringsAsFactors = FALSE)) stopifnot(is.data.frame(dxi), identical(mstr$str, dxi[,"str"]), identical(mstr[j2], mfi[j2]), identical(ef, .mpfr2exp(kxi))) dxi ## 2019-08-09: again *varying* size of 'str' rather than only growing !! ## Show that *order* no longer matters: n <- length(ixk <- rev(kxi)) dix <- cbind(x = asNumeric(ixk), prec = .getPrec(ixk), as.data.frame(.mpfr2str(ixk), stringsAsFactors = FALSE))[n:1,] attr(dix, "row.names") <- .set_row_names(n) stopifnot(identical(dxi, dix)) ## somewhat (but not so much) revealing : cbind(prec = .getPrec(kxi), kxi = asNumeric(kxi), str = es, fi.10 = ceiling(ef/log2(10)), str.2 = as.integer(es*log2(10)), fi = ef) ## Bug example from RMH 2018-03-16 : (x <- mpfr(c(65, 650, 6500, 65000, 650000), precBits=6)) data.frame(fDec = formatDec(x), f = formatMpfr(x)) x. <- as.numeric(xDec <- formatDec(x)) stopifnot(abs(x - x.) <= c(0, 0, 2, 12, 360)) cat("Checking compatibility .mpfr_formatinfo() <--> .mpfr2str(*, base=2) :\n") for(nm in ls()) if(is(OO <- get(nm), "mpfr")) { cat(nm,": str(*) :\n"); str(OO); cat("compatibility: ") I <- .mpfr_formatinfo(OO) S <- .mpfr2str(OO, base = 2L) if(identical(I, S[-1])) cat("[Ok]\n") else { if(any(B <- !I$finite)) I$exp[B] <- S$exp[B] if(any(B <- I $ is.0)) I$exp[B] <- S$exp[B] if(identical(I, S[-1])) cat(" after fixup [Ok]\n") else stop(".mpfr_formatinfo(*) and .mpfr2str(*, base=2) do not match") } }