R Under development (unstable) (2024-01-20 r85814 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 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. > #### Low level stuff - debugging etc > #### ========= ========= > > 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 > 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 4 'mpfr' numbers of precision 32 bits [1] NaN NaN -Inf Inf > 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) [1] TRUE > > 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 1 'mpfr' number of precision 100 bits .mpfr_debug[1]: mpfr2str(*, digits=0, maybeF=False, base=10): .mpfr_debug[1]: [i=0]: prec=100, exp2=-2807 -> (nchar_i,dig.n)=(33,33) .. max_nchar=33 [1] 7.9792116643192417444751621015671e-846 > ## (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") .mpfr_debug[1]: R_mpfr_get_erange(1): -1073741823 .mpfr_debug[1]: R_mpfr_get_erange(2): 1073741823 .mpfr_debug[1]: R_mpfr_get_erange(3): -1073741823 .mpfr_debug[1]: R_mpfr_get_erange(4): 1073741823 .mpfr_debug[1]: R_mpfr_get_erange(5): -1073741823 .mpfr_debug[1]: R_mpfr_get_erange(6): 1073741823 Emin Emax min.emin max.emin min.emax "-1073741823" "1073741823" "-1073741823" "1073741823" "-1073741823" max.emax "1073741823" > stopifnot(identical(.eranges(), .erange2())) .mpfr_debug[1]: R_mpfr_get_erange(1): -1073741823 .mpfr_debug[1]: R_mpfr_get_erange(2): 1073741823 .mpfr_debug[1]: R_mpfr_get_erange(3): -1073741823 .mpfr_debug[1]: R_mpfr_get_erange(4): 1073741823 .mpfr_debug[1]: R_mpfr_get_erange(5): -1073741823 .mpfr_debug[1]: R_mpfr_get_erange(6): 1073741823 .mpfr_debug[1]: R_mpfr_get_erange(1): -1073741823 .mpfr_debug[1]: R_mpfr_get_erange(2): 1073741823 .mpfr_debug[1]: R_mpfr_get_erange(3): -1073741823 .mpfr_debug[1]: R_mpfr_get_erange(4): 1073741823 .mpfr_debug[1]: R_mpfr_get_erange(5): -1073741823 .mpfr_debug[1]: R_mpfr_get_erange(6): 1073741823 > > .mpfr_minPrec() .mpfr_debug[1]: R_mpfr_prec_range(): 1 [1] 1 > .mpfr_maxPrec()# debug printing shows the long integer (on 64 bit) .mpfr_debug[1]: R_mpfr_prec_range(): 2147483391 [1] 2147483391 > > ## Now, level 2 : > stopifnot(Rmpfr:::.mpfr_debug(2) == 1) > r 1 'mpfr' number of precision 100 bits .mpfr_debug[2]: mpfr2str(*, digits=0, maybeF=False, base=10): .mpfr_debug[2]: ex[0:1]= (4294964489,4294967295) -> _exp = 0xfffff509 .mpfr_debug[2]: dd[0:1]= (2684354560, 135379182) -> r..d[i=0]= 0xa0000000 .mpfr_debug[2]: dd[2:3]= (1395976304,3358285975) -> r..d[i=1]= 0x5334e870 .mpfr_debug[2]: [i=0]: prec=100, exp2=-2807 -> (nchar_i,dig.n)=(33,33) .. max_nchar=33 .mpfr_debug[2]: ex[0:1]= (4294964489,4294967295) -> _exp = 0xfffff509 .mpfr_debug[2]: dd[0:1]= (2684354560, 135379182) -> r..d[i=0]= 0xa0000000 .mpfr_debug[2]: dd[2:3]= (1395976304,3358285975) -> r..d[i=1]= 0x5334e870 [1] 7.9792116643192417444751621015671e-846 > ## 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 .mpfr_debug[2]: ex[0:1]= (4294964489,4294967295) -> _exp = 0xfffff509 .mpfr_debug[2]: dd[0:1]= (2684354560, 135379182) -> r..d[i=0]= 0xa0000000 .mpfr_debug[2]: dd[2:3]= (1395976304,3358285975) -> r..d[i=1]= 0x5334e870 .mpfr_debug[2]: _exp = 0xfffbb761 .mpfr_debug[2]: r..d[i=0] = 0x70000000 .mpfr_debug[2]: r..d[i=1] = 0x1cfbf9bc > r2 1 'mpfr' number of precision 100 bits .mpfr_debug[2]: mpfr2str(*, digits=0, maybeF=False, base=10): .mpfr_debug[2]: ex[0:1]= (4294686561,4294967295) -> _exp = 0xfffbb761 .mpfr_debug[2]: dd[0:1]= (1879048192,3308289136) -> r..d[i=0]= 0x70000000 .mpfr_debug[2]: dd[2:3]= ( 486275516,3053452464) -> r..d[i=1]= 0x1cfbf9bc .mpfr_debug[2]: [i=0]: prec=100, exp2=-280735 -> (nchar_i,dig.n)=(33,33) .. max_nchar=33 .mpfr_debug[2]: ex[0:1]= (4294686561,4294967295) -> _exp = 0xfffbb761 .mpfr_debug[2]: dd[0:1]= (1879048192,3308289136) -> r..d[i=0]= 0x70000000 .mpfr_debug[2]: dd[2:3]= ( 486275516,3053452464) -> r..d[i=1]= 0x1cfbf9bc [1] 1.5703576492231738528268016707234e-84510 > L <- r^-100000 .mpfr_debug[2]: ex[0:1]= (4294964489,4294967295) -> _exp = 0xfffff509 .mpfr_debug[2]: dd[0:1]= (2684354560, 135379182) -> r..d[i=0]= 0xa0000000 .mpfr_debug[2]: dd[2:3]= (1395976304,3358285975) -> r..d[i=1]= 0x5334e870 .mpfr_debug[2]: _exp = 0x10bbaf05 .mpfr_debug[2]: r..d[i=0] = 0x20000000 .mpfr_debug[2]: r..d[i=1] = 0x1dc591a5 > L3 <- L^3 .mpfr_debug[2]: ex[0:1]= ( 280735493, 0) -> _exp = 0x10bbaf05 .mpfr_debug[2]: dd[0:1]= ( 536870912,3888314300) -> r..d[i=0]= 0x20000000 .mpfr_debug[2]: dd[2:3]= ( 499487141,2476680137) -> r..d[i=1]= 0x1dc591a5 .mpfr_debug[2]: _exp = 0x32330d0d .mpfr_debug[2]: r..d[i=0] = 0x10000000 .mpfr_debug[2]: r..d[i=1] = 0x58f361ce > 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 > ## 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) Class 'mpfr' [package "Rmpfr"] of length 1 and precision 100 .mpfr_debug[2]: mpfr2str(*, digits=12, maybeF=TRUE, base=10): .mpfr_debug[2]: ex[0:1]= ( 842206477, 0) -> _exp = 0x32330d0d .mpfr_debug[2]: dd[0:1]= ( 268435456, 761715680) -> r..d[i=0]= 0x10000000 .mpfr_debug[2]: dd[2:3]= (1492345294,3294200526) -> r..d[i=1]= 0x58f361ce .mpfr_debug[2]: N_digits: [i=0]: ... -> dig.n = 12 .. max_nchar=12 .mpfr_debug[2]: ex[0:1]= ( 842206477, 0) -> _exp = 0x32330d0d .mpfr_debug[2]: dd[0:1]= ( 268435456, 761715680) -> r..d[i=0]= 0x10000000 .mpfr_debug[2]: dd[2:3]= (1492345294,3294200526) -> r..d[i=1]= 0x58f361ce .mpfr_debug[2]: mpfr2str(*, digits=12, maybeF=TRUE, base=10): .mpfr_debug[2]: ex[0:1]= ( 842206477, 0) -> _exp = 0x32330d0d .mpfr_debug[2]: dd[0:1]= ( 268435456, 761715680) -> r..d[i=0]= 0x10000000 .mpfr_debug[2]: dd[2:3]= (1492345294,3294200526) -> r..d[i=1]= 0x58f361ce .mpfr_debug[2]: N_digits: [i=0]: ... -> dig.n = 12 .. max_nchar=12 .mpfr_debug[2]: ex[0:1]= ( 842206477, 0) -> _exp = 0x32330d0d .mpfr_debug[2]: dd[0:1]= ( 268435456, 761715680) -> r..d[i=0]= 0x10000000 .mpfr_debug[2]: dd[2:3]= (1492345294,3294200526) -> r..d[i=1]= 0x58f361ce 1.00989692356e+253529412 > ## 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) 1 'mpfr' number of precision 100 bits .mpfr_debug[2]: mpfr2str(*, digits=0, maybeF=False, base=10): .mpfr_debug[2]: ex[0:1]= ( 842206477, 0) -> _exp = 0x32330d0d .mpfr_debug[2]: dd[0:1]= ( 268435456, 761715680) -> r..d[i=0]= 0x10000000 .mpfr_debug[2]: dd[2:3]= (1492345294,3294200526) -> r..d[i=1]= 0x58f361ce .mpfr_debug[2]: [i=0]: prec=100, exp2=842206477 -> (nchar_i,dig.n)=(33,33) .. max_nchar=33 .mpfr_debug[2]: ex[0:1]= ( 842206477, 0) -> _exp = 0x32330d0d .mpfr_debug[2]: dd[0:1]= ( 268435456, 761715680) -> r..d[i=0]= 0x10000000 .mpfr_debug[2]: dd[2:3]= (1492345294,3294200526) -> r..d[i=1]= 0x58f361ce [1] 1.0098969235574373617534053306205e+253529412 > ## 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)) List of 4 $ :List of 4 ..$ prec: int 32 ..$ exp : int [1:2] -2147483646 -1 ..$ sign: int 1 ..$ d : int(0) $ :List of 4 ..$ prec: int 32 ..$ exp : int [1:2] -2147483646 -1 ..$ sign: int 1 ..$ d : int(0) $ :List of 4 ..$ prec: int 32 ..$ exp : int [1:2] -2147483645 -1 ..$ sign: int -1 ..$ d : int(0) $ :List of 4 ..$ prec: int 32 ..$ exp : int [1:2] -2147483645 -1 ..$ sign: int 1 ..$ d : int(0) > ## slightly nicer ["uniformly not worse"] (still very similar) : > str(x4, internal=TRUE) Class 'mpfr' [package "Rmpfr"] of length 4 and precision 32 internally @.Data: List of 4 $ :Formal class 'mpfr1' [package "Rmpfr"] with 4 slots .. ..@ prec: int 32 .. ..@ exp : int [1:2] -2147483646 -1 .. ..@ sign: int 1 .. ..@ d : int(0) $ :Formal class 'mpfr1' [package "Rmpfr"] with 4 slots .. ..@ prec: int 32 .. ..@ exp : int [1:2] -2147483646 -1 .. ..@ sign: int 1 .. ..@ d : int(0) $ :Formal class 'mpfr1' [package "Rmpfr"] with 4 slots .. ..@ prec: int 32 .. ..@ exp : int [1:2] -2147483645 -1 .. ..@ sign: int -1 .. ..@ d : int(0) $ :Formal class 'mpfr1' [package "Rmpfr"] with 4 slots .. ..@ prec: int 32 .. ..@ exp : int [1:2] -2147483645 -1 .. ..@ sign: int 1 .. ..@ d : int(0) > x4 ## "similar info" as .mpfr2list(.) 4 'mpfr' numbers of precision 32 bits .mpfr_debug[2]: mpfr2str(*, digits=0, maybeF=False, base=10): .mpfr_debug[2]: ex[0:1]= (2147483650,4294967295) -> _exp = 0x80000002 .. max_nchar=5 .mpfr_debug[2]: ex[0:1]= (2147483650,4294967295) -> _exp = 0x80000002 .. max_nchar=5 .mpfr_debug[2]: ex[0:1]= (2147483651,4294967295) -> _exp = 0x80000003 .. max_nchar=5 .mpfr_debug[2]: ex[0:1]= (2147483651,4294967295) -> _exp = 0x80000003 .. max_nchar=5 .mpfr_debug[2]: ex[0:1]= (2147483650,4294967295) -> _exp = 0x80000002 .mpfr_debug[2]: ex[0:1]= (2147483650,4294967295) -> _exp = 0x80000002 .mpfr_debug[2]: ex[0:1]= (2147483651,4294967295) -> _exp = 0x80000003 .mpfr_debug[2]: ex[0:1]= (2147483651,4294967295) -> _exp = 0x80000003 [1] NaN NaN -Inf Inf > > ## Increase maximal exponent: > > tools:::assertWarning( + .mpfr_erange_set("Emax", 5e18)) # too large {FIXME why only warning and not error ??} > .mpfr_erange("Emax") # is unchanged .mpfr_debug[2]: R_mpfr_get_erange(2): 1073741823 Emax 1073741823 > if(4e18 < .mpfr_erange("max.emax")) { + .mpfr_erange_set("Emax", 4e18) # now ok: + stopifnot(.mpfr_erange("Emax") == 4e18) + } .mpfr_debug[2]: R_mpfr_get_erange(6): 1073741823 > > > ## revert to no debugging: > stopifnot(Rmpfr:::.mpfr_debug(0) == 2) > .mpfr_maxPrec() [1] 2147483391 > > L / (r2^-1000)# 1.00000....448 (could be more accurate?) 1 'mpfr' number of precision 100 bits [1] 1.0000000000000000000000000004481 > > 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)) List of 5 $ gmp.numb.bits: int 64 $ mpfr.version : chr "4.2.1" $ Machine :List of 5 ..$ sizeof.long : int 4 ..$ sizeof.longlong : int 8 ..$ sizeof.longdouble: int 16 ..$ sizeof.pointer : int 8 ..$ sizeof.time_t : int 8 $ Sys.info : Named chr [1:2] "Windows" "x86-64" ..- attr(*, "names")= chr [1:2] "sysname" "machine" $ mpfr :List of 3 ..$ :List of 4 .. ..$ prec: int 100 .. ..$ exp : int [1:2] -2807 -1 .. ..$ sign: int 1 .. ..$ d : int [1:4] -1610612736 135379182 1395976304 -936681321 ..$ :List of 4 .. ..$ prec: int 100 .. ..$ exp : int [1:2] -280735 -1 .. ..$ sign: int 1 .. ..$ d : int [1:4] 1879048192 -986678160 486275516 -1241514832 ..$ :List of 4 .. ..$ prec: int 100 .. ..$ exp : int [1:2] 280735493 0 .. ..$ sign: int 1 .. ..$ d : int [1:4] 536870912 -406652996 499487141 -1818287159 - attr(*, "class")= chr "mpfrXport" > str(mpfrXport(mpfr(2, 64)^(-3:3))) List of 5 $ gmp.numb.bits: int 64 $ mpfr.version : chr "4.2.1" $ Machine :List of 5 ..$ sizeof.long : int 4 ..$ sizeof.longlong : int 8 ..$ sizeof.longdouble: int 16 ..$ sizeof.pointer : int 8 ..$ sizeof.time_t : int 8 $ Sys.info : Named chr [1:2] "Windows" "x86-64" ..- attr(*, "names")= chr [1:2] "sysname" "machine" $ mpfr :List of 7 ..$ :List of 4 .. ..$ prec: int 64 .. ..$ exp : int [1:2] -2 -1 .. ..$ sign: int 1 .. ..$ d : int [1:2] 0 NA ..$ :List of 4 .. ..$ prec: int 64 .. ..$ exp : int [1:2] -1 -1 .. ..$ sign: int 1 .. ..$ d : int [1:2] 0 NA ..$ :List of 4 .. ..$ prec: int 64 .. ..$ exp : int [1:2] 0 0 .. ..$ sign: int 1 .. ..$ d : int [1:2] 0 NA ..$ :List of 4 .. ..$ prec: int 64 .. ..$ exp : int [1:2] 1 0 .. ..$ sign: int 1 .. ..$ d : int [1:2] 0 NA ..$ :List of 4 .. ..$ prec: int 64 .. ..$ exp : int [1:2] 2 0 .. ..$ sign: int 1 .. ..$ d : int [1:2] 0 NA ..$ :List of 4 .. ..$ prec: int 64 .. ..$ exp : int [1:2] 3 0 .. ..$ sign: int 1 .. ..$ d : int [1:2] 0 NA ..$ :List of 4 .. ..$ prec: int 64 .. ..$ exp : int [1:2] 4 0 .. ..$ sign: int 1 .. ..$ d : int [1:2] 0 NA - attr(*, "class")= chr "mpfrXport" > str(mpfrXport(Const("pi")* 2^(-3:3))) List of 5 $ gmp.numb.bits: int 64 $ mpfr.version : chr "4.2.1" $ Machine :List of 5 ..$ sizeof.long : int 4 ..$ sizeof.longlong : int 8 ..$ sizeof.longdouble: int 16 ..$ sizeof.pointer : int 8 ..$ sizeof.time_t : int 8 $ Sys.info : Named chr [1:2] "Windows" "x86-64" ..- attr(*, "names")= chr [1:2] "sysname" "machine" $ mpfr :List of 7 ..$ :List of 4 .. ..$ prec: int 120 .. ..$ exp : int [1:2] -1 -1 .. ..$ sign: int 1 .. ..$ d : int [1:4] -2133058304 -993631605 560513588 -921707870 ..$ :List of 4 .. ..$ prec: int 120 .. ..$ exp : int [1:2] 0 0 .. ..$ sign: int 1 .. ..$ d : int [1:4] -2133058304 -993631605 560513588 -921707870 ..$ :List of 4 .. ..$ prec: int 120 .. ..$ exp : int [1:2] 1 0 .. ..$ sign: int 1 .. ..$ d : int [1:4] -2133058304 -993631605 560513588 -921707870 ..$ :List of 4 .. ..$ prec: int 120 .. ..$ exp : int [1:2] 2 0 .. ..$ sign: int 1 .. ..$ d : int [1:4] -2133058304 -993631605 560513588 -921707870 ..$ :List of 4 .. ..$ prec: int 120 .. ..$ exp : int [1:2] 3 0 .. ..$ sign: int 1 .. ..$ d : int [1:4] -2133058304 -993631605 560513588 -921707870 ..$ :List of 4 .. ..$ prec: int 120 .. ..$ exp : int [1:2] 4 0 .. ..$ sign: int 1 .. ..$ d : int [1:4] -2133058304 -993631605 560513588 -921707870 ..$ :List of 4 .. ..$ prec: int 120 .. ..$ exp : int [1:2] 5 0 .. ..$ sign: int 1 .. ..$ d : int [1:4] -2133058304 -993631605 560513588 -921707870 - attr(*, "class")= chr "mpfrXport" > > ## and a very large one > mil <- mpfr(1025, 111) > str(mm <- mpfrXport(xx <- mil^(2^25))) List of 5 $ gmp.numb.bits: int 64 $ mpfr.version : chr "4.2.1" $ Machine :List of 5 ..$ sizeof.long : int 4 ..$ sizeof.longlong : int 8 ..$ sizeof.longdouble: int 16 ..$ sizeof.pointer : int 8 ..$ sizeof.time_t : int 8 $ Sys.info : Named chr [1:2] "Windows" "x86-64" ..- attr(*, "names")= chr [1:2] "sysname" "machine" $ mpfr :List of 1 ..$ :List of 4 .. ..$ prec: int 111 .. ..$ exp : int [1:2] 335591572 0 .. ..$ sign: int 1 .. ..$ d : int [1:4] -474349568 -2023461216 -832488434 -1890623142 - attr(*, "class")= chr "mpfrXport" > 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"))) + } Cannot set 'Emin' to -2^40 (= -1.09951e+12), as .mpfr_erange("min.emin") is larger, namely -1.07374e+09. > > ## 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)) [,1] [1,] "0.05273" [2,] "0.052631578947" [3,] "0.0526315793" > 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 [1] 128 1024 12288 131072 > str(sk1 <- .mpfr2str(k1)) List of 4 $ str : chr [1:4] "13" "10" "12" "13" $ exp : int [1:4] 3 4 5 6 $ finite: logi [1:4] TRUE TRUE TRUE TRUE $ is.0 : logi [1:4] FALSE FALSE FALSE FALSE > str(sk1. <- .mpfr2str(k1, maybe.full=TRUE)) List of 4 $ str : chr [1:4] "128" "1024" "12288" "131072" $ exp : int [1:4] 3 4 5 6 $ finite: logi [1:4] TRUE TRUE TRUE TRUE $ is.0 : logi [1:4] FALSE FALSE FALSE FALSE > str(sk1.2 <- .mpfr2str(k1, digits=2, base=2)) List of 4 $ str : chr [1:4] "10" "10" "11" "10" $ exp : int [1:4] 8 11 14 18 $ finite: logi [1:4] TRUE TRUE TRUE TRUE $ is.0 : logi [1:4] FALSE FALSE FALSE FALSE > str(sk1.2F <- .mpfr2str(k1, maybe.full=TRUE, base=2)) List of 4 $ str : chr [1:4] "10000000" "10000000000" "11000000000000" "100000000000000000" $ exp : int [1:4] 8 11 14 18 $ finite: logi [1:4] TRUE TRUE TRUE TRUE $ is.0 : logi [1:4] FALSE FALSE FALSE FALSE > 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 ???? [1] "128." "1024." "12288." "131072." > ## was "128." "1024." "12288." "131072." , but now obeying internal precision gives > ## "1.e2" "1.e3" "1.e4" "1.e5" > (fk1 <- formatMpfr(k1, digits=6)) [1] "128.000" "1024.00" "12288.0" "131072." > 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)) Call: lm(formula = log10(colMeans(abs(rE)))[i] ~ i) Residuals: 1 2 3 4 5 -0.12076 0.06070 0.10167 0.09761 -0.13922 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.49710 0.14522 3.423 0.041747 * i -1.15528 0.04378 -26.386 0.000119 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.1385 on 3 degrees of freedom Multiple R-squared: 0.9957, Adjusted R-squared: 0.9943 F-statistic: 696.2 on 1 and 3 DF, p-value: 0.0001194 > 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)) 64 'mpfr' numbers of precision 2 .. 64 bits [1] NaN NaN [3] -Inf -2 [5] -1.2000000000000002 -1.1000000000000001 [7] -1 -1 [9] -0.89999999999999991 -0.80000000000000004 [11] -0.69999999999999996 -0.60000000000000009 [13] -0.5 -0.40000000000000002 [15] -0.30000000000000004 -0.20000000000000001 [17] -0.109999999999999999999 -0.10000000000000001 [19] -0.100000000000000000001 -0.0900000000000000000033 [21] -0.0799999999999999999984 -0.0700000000000000000003 [23] -0.0599999999999999999988 -0.0500000000000000000007 [25] -0.0399999999999999999992 -0.0299999999999999999994 [27] -0.0199999999999999999996 -0.0099999999999999999998 [29] 0 0 [31] 0 0.0099999999999999999998 [33] 0.0199999999999999999996 0.0299999999999999999994 [35] 0.0399999999999999999992 0.0500000000000000000007 [37] 0.0599999999999999999988 0.0700000000000000000003 [39] 0.0799999999999999999984 0.0900000000000000000033 [41] 0.100000000000000000001 0.10000000000000001 [43] 0.109999999999999999999 0.20000000000000001 [45] 0.30000000000000004 0.40000000000000002 [47] 0.5 0.60000000000000009 [49] 0.69999999999999996 0.80000000000000004 [51] 0.89999999999999991 1 [53] 1 1.1000000000000001 [55] 1.2000000000000002 2 [57] 3 4 [59] 5 130 [61] 1e+3 1.2e+4 [63] 1.3e+5 Inf > 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 !! x prec str exp finite is.0 1 NaN 32 @NaN@ 0 FALSE FALSE 2 NaN 32 @NaN@ 0 FALSE FALSE 3 -Inf 32 -@Inf@ 0 FALSE FALSE 4 -2.00 32 -20000000000 1 TRUE FALSE 5 -1.20 52 -12000000000000002 1 TRUE FALSE 6 -1.10 52 -11000000000000001 1 TRUE FALSE 7 -1.00 32 -10000000000 1 TRUE FALSE 8 -1.00 52 -10000000000000000 1 TRUE FALSE 9 -0.90 52 -89999999999999991 0 TRUE FALSE 10 -0.80 52 -80000000000000004 0 TRUE FALSE 11 -0.70 52 -69999999999999996 0 TRUE FALSE 12 -0.60 52 -60000000000000009 0 TRUE FALSE 13 -0.50 52 -50000000000000000 0 TRUE FALSE 14 -0.40 52 -40000000000000002 0 TRUE FALSE 15 -0.30 52 -30000000000000004 0 TRUE FALSE 16 -0.20 52 -20000000000000001 0 TRUE FALSE 17 -0.11 64 -109999999999999999999 0 TRUE FALSE 18 -0.10 52 -10000000000000001 0 TRUE FALSE 19 -0.10 64 -100000000000000000001 0 TRUE FALSE 20 -0.09 64 -900000000000000000033 -1 TRUE FALSE 21 -0.08 64 -799999999999999999984 -1 TRUE FALSE 22 -0.07 64 -700000000000000000003 -1 TRUE FALSE 23 -0.06 64 -599999999999999999988 -1 TRUE FALSE 24 -0.05 64 -500000000000000000007 -1 TRUE FALSE 25 -0.04 64 -399999999999999999992 -1 TRUE FALSE 26 -0.03 64 -299999999999999999994 -1 TRUE FALSE 27 -0.02 64 -199999999999999999996 -1 TRUE FALSE 28 -0.01 64 -999999999999999999980 -2 TRUE FALSE 29 0.00 32 0 0 TRUE TRUE 30 0.00 52 0 0 TRUE TRUE 31 0.00 64 0 0 TRUE TRUE 32 0.01 64 999999999999999999980 -2 TRUE FALSE 33 0.02 64 199999999999999999996 -1 TRUE FALSE 34 0.03 64 299999999999999999994 -1 TRUE FALSE 35 0.04 64 399999999999999999992 -1 TRUE FALSE 36 0.05 64 500000000000000000007 -1 TRUE FALSE 37 0.06 64 599999999999999999988 -1 TRUE FALSE 38 0.07 64 700000000000000000003 -1 TRUE FALSE 39 0.08 64 799999999999999999984 -1 TRUE FALSE 40 0.09 64 900000000000000000033 -1 TRUE FALSE 41 0.10 64 100000000000000000001 0 TRUE FALSE 42 0.10 52 10000000000000001 0 TRUE FALSE 43 0.11 64 109999999999999999999 0 TRUE FALSE 44 0.20 52 20000000000000001 0 TRUE FALSE 45 0.30 52 30000000000000004 0 TRUE FALSE 46 0.40 52 40000000000000002 0 TRUE FALSE 47 0.50 52 50000000000000000 0 TRUE FALSE 48 0.60 52 60000000000000009 0 TRUE FALSE 49 0.70 52 69999999999999996 0 TRUE FALSE 50 0.80 52 80000000000000004 0 TRUE FALSE 51 0.90 52 89999999999999991 0 TRUE FALSE 52 1.00 32 10000000000 1 TRUE FALSE 53 1.00 52 10000000000000000 1 TRUE FALSE 54 1.10 52 11000000000000001 1 TRUE FALSE 55 1.20 52 12000000000000002 1 TRUE FALSE 56 2.00 32 20000000000 1 TRUE FALSE 57 3.00 32 30000000000 1 TRUE FALSE 58 4.00 32 40000000000 1 TRUE FALSE 59 5.00 32 50000000000 1 TRUE FALSE 60 128.00 2 13 3 TRUE FALSE 61 1024.00 2 10 4 TRUE FALSE 62 12288.00 2 12 5 TRUE FALSE 63 131072.00 2 13 6 TRUE FALSE 64 Inf 32 @Inf@ 0 FALSE FALSE > ## 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) prec kxi str fi.10 str.2 fi [1,] 32 NaN 0 -646456992 0 -2147483646 [2,] 32 NaN 0 -646456992 0 -2147483646 [3,] 32 -Inf 0 -646456992 0 -2147483645 [4,] 32 -2.00 1 1 3 2 [5,] 52 -1.20 1 1 3 1 [6,] 52 -1.10 1 1 3 1 [7,] 32 -1.00 1 1 3 1 [8,] 52 -1.00 1 1 3 1 [9,] 52 -0.90 0 0 0 0 [10,] 52 -0.80 0 0 0 0 [11,] 52 -0.70 0 0 0 0 [12,] 52 -0.60 0 0 0 0 [13,] 52 -0.50 0 0 0 0 [14,] 52 -0.40 0 0 0 -1 [15,] 52 -0.30 0 0 0 -1 [16,] 52 -0.20 0 0 0 -2 [17,] 64 -0.11 0 0 0 -3 [18,] 52 -0.10 0 0 0 -3 [19,] 64 -0.10 0 0 0 -3 [20,] 64 -0.09 -1 0 -3 -3 [21,] 64 -0.08 -1 0 -3 -3 [22,] 64 -0.07 -1 0 -3 -3 [23,] 64 -0.06 -1 -1 -3 -4 [24,] 64 -0.05 -1 -1 -3 -4 [25,] 64 -0.04 -1 -1 -3 -4 [26,] 64 -0.03 -1 -1 -3 -5 [27,] 64 -0.02 -1 -1 -3 -5 [28,] 64 -0.01 -2 -1 -6 -6 [29,] 32 0.00 0 -646456992 0 -2147483647 [30,] 52 0.00 0 -646456992 0 -2147483647 [31,] 64 0.00 0 -646456992 0 -2147483647 [32,] 64 0.01 -2 -1 -6 -6 [33,] 64 0.02 -1 -1 -3 -5 [34,] 64 0.03 -1 -1 -3 -5 [35,] 64 0.04 -1 -1 -3 -4 [36,] 64 0.05 -1 -1 -3 -4 [37,] 64 0.06 -1 -1 -3 -4 [38,] 64 0.07 -1 0 -3 -3 [39,] 64 0.08 -1 0 -3 -3 [40,] 64 0.09 -1 0 -3 -3 [41,] 64 0.10 0 0 0 -3 [42,] 52 0.10 0 0 0 -3 [43,] 64 0.11 0 0 0 -3 [44,] 52 0.20 0 0 0 -2 [45,] 52 0.30 0 0 0 -1 [46,] 52 0.40 0 0 0 -1 [47,] 52 0.50 0 0 0 0 [48,] 52 0.60 0 0 0 0 [49,] 52 0.70 0 0 0 0 [50,] 52 0.80 0 0 0 0 [51,] 52 0.90 0 0 0 0 [52,] 32 1.00 1 1 3 1 [53,] 52 1.00 1 1 3 1 [54,] 52 1.10 1 1 3 1 [55,] 52 1.20 1 1 3 1 [56,] 32 2.00 1 1 3 2 [57,] 32 3.00 1 1 3 2 [58,] 32 4.00 1 1 3 3 [59,] 32 5.00 1 1 3 3 [60,] 2 128.00 3 3 9 8 [61,] 2 1024.00 4 4 13 11 [62,] 2 12288.00 5 5 16 14 [63,] 2 131072.00 6 6 19 18 [64,] 32 Inf 0 -646456992 0 -2147483645 > > > > ## Bug example from RMH 2018-03-16 : > (x <- mpfr(c(65, 650, 6500, 65000, 650000), precBits=6)) 5 'mpfr' numbers of precision 6 bits [1] 64 656 6530 6.45e+4 6.55e+5 > data.frame(fDec = formatDec(x), f = formatMpfr(x)) fDec f 1 64.0 64.0 2 656. 656. 3 6530. 6530. 4 6.45e+4 6.45e+4 5 6.55e+5 6.55e+5 > 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") Checking compatibility .mpfr_formatinfo() <--> .mpfr2str(*, base=2) : > 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") + } + } L : str(*) : Class 'mpfr' [package "Rmpfr"] of length 1 and precision 100 1.00328815073e+84509804 compatibility: [Ok] L3 : str(*) : Class 'mpfr' [package "Rmpfr"] of length 1 and precision 100 1.00989692356e+253529412 compatibility: [Ok] cents : str(*) : Class 'mpfr' [package "Rmpfr"] of length 23 and precision 64 -0.11 -0.1 -0.09 -0.08 -0.07 -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0 ... compatibility: after fixup [Ok] i8 : str(*) : Class 'mpfr' [package "Rmpfr"] of length 8 and precision 32 -2 -1 0 1 2 3 4 5 compatibility: after fixup [Ok] ixk : str(*) : Class 'mpfr' [package "Rmpfr"] of length 64 and precisions 2..64 Inf 131072 12288 ... compatibility: after fixup [Ok] k1 : str(*) : Class 'mpfr' [package "Rmpfr"] of length 4 and precision 2 1.e+2 1.e+3 1.e+4 1.e+5 compatibility: [Ok] kxi : str(*) : Class 'mpfr' [package "Rmpfr"] of length 64 and precisions 2..64 NaN NaN -Inf -2 -1.2 ... compatibility: after fixup [Ok] mil : str(*) : Class 'mpfr' [package "Rmpfr"] of length 1 and precision 111 1025 compatibility: [Ok] r : str(*) : Class 'mpfr' [package "Rmpfr"] of length 1 and precision 100 7.97921166432e-846 compatibility: [Ok] r2 : str(*) : Class 'mpfr' [package "Rmpfr"] of length 1 and precision 100 1.57035764922e-84510 compatibility: [Ok] rv : str(*) : Class 'mpfr' [package "Rmpfr"] of length 3 and precision 100 7.97921166432e-846 1.57035764922e-84510 1.00328815073e+84509804 compatibility: [Ok] tenth : str(*) : Class 'mpfr' [package "Rmpfr"] of length 25 and precision 52 -1.2 -1.1 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 ... compatibility: after fixup [Ok] x : str(*) : Class 'mpfr' [package "Rmpfr"] of length 5 and precision 6 64 660 6.5e+3 6.5e+4 6.6e+5 compatibility: [Ok] x4 : str(*) : Class 'mpfr' [package "Rmpfr"] of length 4 and precision 32 NaN NaN -Inf Inf compatibility: after fixup [Ok] xx : str(*) : Class 'mpfr' [package "Rmpfr"] of length 1 and precision 111 1.62954142022e+101023129 compatibility: [Ok] z : str(*) : Class 'mpfr' [package "Rmpfr"] of length 3 and precisions 8..32 0.052734375 0.05263157895 0.05263157934 compatibility: [Ok] > > proc.time() user system elapsed 0.68 0.12 0.79