R Under development (unstable) (2023-12-04 r85659 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(round) > > ## Utilities (sync with ../man/roundX.Rd ! ) > > ##' maybe add to package -- it's really useful here : > formatF <- function(x, ...) noquote(format(x, scientific=FALSE, drop0trailing=TRUE, ...)) > > same_cols <- function(m) all(m == m[,1]) > matPm <- function(y) { + matplot(y=y, type = "l", lwd = 2, xlab = "i", ylab = deparse(substitute(y))) + abline(h = 0, lty=2, col="gray") + legend("topleft", legend = setdiff(roundVersions, "r3"), + col = 1:6, lty = 1:5, lwd = 2, bty = "n") + } > > str(.M <- .Machine, digits = 8) List of 29 $ double.eps : num 2.220446e-16 $ double.neg.eps : num 1.110223e-16 $ double.xmin : num 2.2250739e-308 $ double.xmax : num 1.7976931e+308 $ double.base : int 2 $ double.digits : int 53 $ double.rounding : int 5 $ double.guard : int 0 $ double.ulp.digits : int -52 $ double.neg.ulp.digits : int -53 $ double.exponent : int 11 $ double.min.exp : int -1022 $ double.max.exp : int 1024 $ integer.max : int 2147483647 $ sizeof.long : int 4 $ sizeof.longlong : int 8 $ sizeof.longdouble : int 16 $ sizeof.pointer : int 8 $ sizeof.time_t : int 8 $ longdouble.eps : num 1.0842022e-19 $ longdouble.neg.eps : num 5.4210109e-20 $ longdouble.digits : int 64 $ longdouble.rounding : int 5 $ longdouble.guard : int 0 $ longdouble.ulp.digits : int -63 $ longdouble.neg.ulp.digits: int -64 $ longdouble.exponent : int 15 $ longdouble.min.exp : int -16382 $ longdouble.max.exp : int 16384 > capabilities() jpeg png tiff tcltk X11 aqua TRUE TRUE TRUE TRUE FALSE FALSE http/ftp sockets libxml fifo cledit iconv TRUE TRUE FALSE TRUE FALSE TRUE NLS Rprof profmem cairo ICU long.double TRUE TRUE TRUE TRUE TRUE TRUE libcurl TRUE > str(.M[grep("^sizeof", names(.M))]) ## also differentiate long-double.. List of 5 $ sizeof.long : int 4 $ sizeof.longlong : int 8 $ sizeof.longdouble: int 16 $ sizeof.pointer : int 8 $ sizeof.time_t : int 8 > (b64 <- .M$sizeof.pointer == 8) [1] TRUE > (arch <- Sys.info()[["machine"]]) [1] "x86-64" > (onWindows <- .Platform$OS.type == "windows") [1] TRUE > (hasLD <- unname(capabilities("long.double"))) [1] TRUE > win32 <- onWindows && !b64 > > if(!dev.interactive(orNone=TRUE)) pdf("round-extreme.pdf") > > ## the extreme case > M <- .Machine$double.xmax; dig <- -(1:314) > resM <- roundAll(M, dig) There were 50 or more warnings (use warnings() to see the first 50) > if(length(warnings() -> wrngs)) print(summary(wrngs)) 50 identical warnings: In FUN(X[[i]], ...) : probable complete loss of accuracy in modulus > rownames(resM) <- paste0("dig=", dig) > op <- options(digits = 9, width = 149) > print(head(resM, 9), digits = 7) sprintf r0.C r1.C r1a.C r2.C r2a.C r3.C r3d.C r3 dig=-1 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 dig=-2 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 dig=-3 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 dig=-4 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 dig=-5 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 dig=-6 1.797693e+308 Inf Inf Inf Inf 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 dig=-7 1.797693e+308 Inf Inf Inf Inf 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 dig=-8 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 dig=-9 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 > print(tail(resM, 20)) # "r3*" are best sprintf r0.C r1.C r1a.C r2.C r2a.C r3.C r3d.C r3 dig=-295 Inf 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 dig=-296 Inf 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 dig=-297 Inf 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 dig=-298 Inf Inf Inf Inf Inf Inf 1.79769313e+308 1.79769313e+308 1.79769313e+308 dig=-299 Inf Inf Inf Inf Inf Inf 1.79769313e+308 1.79769313e+308 1.79769313e+308 dig=-300 Inf 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 1.79769313e+308 dig=-301 Inf 1.79769310e+308 1.79769310e+308 1.79769310e+308 1.79769310e+308 1.79769310e+308 1.79769310e+308 1.79769310e+308 1.79769310e+308 dig=-302 Inf 1.79769300e+308 1.79769300e+308 1.79769300e+308 1.79769300e+308 1.79769300e+308 1.79769300e+308 1.79769300e+308 1.79769300e+308 dig=-303 Inf 1.79769000e+308 1.79769000e+308 1.79769000e+308 1.79769000e+308 1.79769000e+308 1.79769000e+308 1.79769000e+308 1.79769000e+308 dig=-304 Inf Inf Inf Inf Inf Inf 1.79760000e+308 1.79760000e+308 1.79760000e+308 dig=-305 Inf Inf Inf Inf Inf Inf 1.79700000e+308 1.79700000e+308 1.79700000e+308 dig=-306 Inf Inf Inf Inf Inf Inf 1.79000000e+308 1.79000000e+308 1.79000000e+308 dig=-307 Inf Inf Inf Inf Inf Inf 1.70000000e+308 1.70000000e+308 1.70000000e+308 dig=-308 Inf Inf Inf Inf Inf Inf 1.00000000e+308 1.00000000e+308 1.00000000e+308 dig=-309 Inf NaN NaN NaN NaN 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 dig=-310 Inf NaN NaN NaN NaN 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 dig=-311 Inf NaN NaN NaN NaN 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 dig=-312 Inf NaN NaN NaN NaN 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 dig=-313 Inf NaN NaN NaN NaN 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 dig=-314 Inf NaN NaN NaN NaN 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 > options(op) > iO <- 1:308 > all.eq.M <- function(x,y, ...) all.equal(x, y, countEQ=TRUE, scale=M, ...) > all.eq.M(resM[iO, "r3.C"], resM[iO,"r3"], tol = 0) # "Mean scaled diff.: 4.58e-17" [1] TRUE > all.eq.M(resM[iO,"r3d.C"], resM[iO,"r3"], tol = 0) # " " " 7.17e-17" [1] TRUE > > ## Windows 32 bit (on Winbuilder, gcc 4.9.x, -O3) is "exception" > stopifnot(exprs = { + is.finite(resM[, c("r3.C", if(!win32) "r3d.C", "r3")]) + all.eq.M(resM[iO, "r3.C"], resM[iO,"r3"], tol = 1e-15) + if(win32) TRUE else all.eq.M(resM[iO,"r3d.C"], resM[iO,"r3"], tol = 1e-15) + }) > ## but there are quite a few "small noise" differences: > table(resM[iO, "r3.C"] == resM[iO,"r3"]) FALSE TRUE 125 183 > ## FALSE TRUE > ## 124 184 (F30 Linux gcc; 64-bit) > table(resM[iO,"r3d.C"] == resM[iO,"r3"]) FALSE TRUE 136 172 > ## 136 172 (F30 Linux gcc; 64-bit) > > > ## Unbelievably, this gives *NO* 'Inf' values (on 32-bit Windows), > ## but the default trace=0 *does* give 'Inf' (on Winbuilder, see above): > roundX(.Machine$double.xmax, -(291:310), "r3d.C", trace = 2) fround3d(1.79769313486232e+308, digits=-291): + l10x + dig > DBL_DIG ( = 15): returning x fround3d(1.79769313486232e+308, digits=-292): + l10x + dig > DBL_DIG ( = 15): returning x fround3d(1.79769313486232e+308, digits=-293): + l10x + dig > DBL_DIG ( = 15): returning x fround3d(1.79769313486232e+308, digits=-294): + l10x=308.104, dig=-294, sign=1 + dig <= 308: pow10=1e-294, x10=1.79769e+14, i10=1.79769e+14 ==> (xd, xu) = (1.79769313486231e+308,inf) on WIN32_ fake action on take_xu ==> 0 + d{u,d} = inf,4.79001674288333e+293; i10 is odd ==> choosing xd fround3d(1.79769313486232e+308, digits=-295): + l10x=308.104, dig=-295, sign=1 + dig <= 308: pow10=1e-295, x10=1.79769e+13, i10=1.79769e+13 ==> (xd, xu) = (1.7976931348623e+308,inf) on WIN32_ fake action on take_xu ==> 0 + d{u,d} = inf,1.45696342596035e+294; i10 is odd ==> choosing xd fround3d(1.79769313486232e+308, digits=-296): + l10x=308.104, dig=-296, sign=1 + dig <= 308: pow10=1e-296, x10=1.79769e+12, i10=1.79769e+12 ==> (xd, xu) = (1.797693134862e+308,inf) on WIN32_ fake action on take_xu ==> 0 + d{u,d} = inf,3.14744016813625e+295; i10 is even ==> choosing xd fround3d(1.79769313486232e+308, digits=-297): + l10x=308.104, dig=-297, sign=1 + dig <= 308: pow10=1e-297, x10=1.79769e+11, i10=1.79769e+11 ==> (xd, xu) = (1.79769313486e+308,inf) on WIN32_ fake action on take_xu ==> 0 + d{u,d} = inf,2.31457600696741e+296; i10 is even ==> choosing xd fround3d(1.79769313486232e+308, digits=-298): + l10x=308.104, dig=-298, sign=1 + dig <= 308: pow10=1e-298, x10=1.79769e+10, i10=1.79769e+10 ==> (xd, xu) = (1.7976931348e+308,inf) on WIN32_ fake action on take_xu ==> 0 + d{u,d} = inf,6.23147248963859e+297; i10 is even ==> choosing xd fround3d(1.79769313486232e+308, digits=-299): + l10x=308.104, dig=-299, sign=1 + dig <= 308: pow10=1e-299, x10=1.79769e+09, i10=1.79769e+09 ==> (xd, xu) = (1.797693134e+308,inf) on WIN32_ fake action on take_xu ==> 0 + d{u,d} = inf,8.62314581192302e+298; i10 is even ==> choosing xd fround3d(1.79769313486232e+308, digits=-300): + l10x=308.104, dig=-300, sign=1 + dig <= 308: pow10=1e-300, x10=1.79769e+08, i10=1.79769e+08 ==> (xd, xu) = (1.79769313e+308,inf) on WIN32_ fake action on take_xu ==> 0 + d{u,d} = inf,4.86231446142398e+299; i10 is odd ==> choosing xd fround3d(1.79769313486232e+308, digits=-301): + l10x=308.104, dig=-301, sign=1 + dig <= 308: pow10=1e-301, x10=1.79769e+07, i10=1.79769e+07 ==> (xd, xu) = (1.7976931e+308,inf) on WIN32_ fake action on take_xu ==> 0 + d{u,d} = inf,3.48623146608737e+300; i10 is odd ==> choosing xd fround3d(1.79769313486232e+308, digits=-302): + l10x=308.104, dig=-302, sign=1 + dig <= 308: pow10=1e-302, x10=1.79769e+06, i10=1.79769e+06 ==> (xd, xu) = (1.797693e+308,inf) on WIN32_ fake action on take_xu ==> 0 + d{u,d} = inf,1.3486231486001e+301; i10 is odd ==> choosing xd fround3d(1.79769313486232e+308, digits=-303): + l10x=308.104, dig=-303, sign=1 + dig <= 308: pow10=1e-303, x10=179769, i10=179769 ==> (xd, xu) = (1.79769e+308,inf) on WIN32_ fake action on take_xu ==> 0 + d{u,d} = inf,3.13486231464699e+302; i10 is odd ==> choosing xd fround3d(1.79769313486232e+308, digits=-304): + l10x=308.104, dig=-304, sign=1 + dig <= 308: pow10=1e-304, x10=17976.9, i10=17976 ==> (xd, xu) = (1.7976e+308,inf) on WIN32_ fake action on take_xu ==> 0 + d{u,d} = inf,9.31348623146432e+303; i10 is even ==> choosing xd fround3d(1.79769313486232e+308, digits=-305): + l10x=308.104, dig=-305, sign=1 + dig <= 308: pow10=1e-305, x10=1797.69, i10=1797 ==> (xd, xu) = (1.797e+308,inf) on WIN32_ fake action on take_xu ==> 0 + d{u,d} = inf,6.93134862314751e+304; i10 is odd ==> choosing xd fround3d(1.79769313486232e+308, digits=-306): + l10x=308.104, dig=-306, sign=1 + dig <= 308: pow10=1e-306, x10=179.769, i10=179 ==> (xd, xu) = (1.79e+308,inf) on WIN32_ fake action on take_xu ==> 0 + d{u,d} = inf,7.69313486231455e+305; i10 is odd ==> choosing xd fround3d(1.79769313486232e+308, digits=-307): + l10x=308.104, dig=-307, sign=1 + dig <= 308: pow10=1e-307, x10=17.9769, i10=17 ==> (xd, xu) = (1.7e+308,inf) on WIN32_ fake action on take_xu ==> 0 + d{u,d} = inf,9.7693134862315e+306; i10 is odd ==> choosing xd fround3d(1.79769313486232e+308, digits=-308): + l10x=308.104, dig=-308, sign=1 + dig <= 308: pow10=1e-308, x10=1.79769, i10=1 ==> (xd, xu) = (1e+308,inf) on WIN32_ fake action on take_xu ==> 0 + d{u,d} = inf,7.97693134862315e+307; i10 is odd ==> choosing xd fround3d(1.79769313486232e+308, digits=-309): fround3d(1.79769313486232e+308, digits=-310): [1] 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 [6] 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 [11] 1.797693e+308 1.797693e+308 1.797690e+308 1.797600e+308 1.797000e+308 [16] 1.790000e+308 1.700000e+308 1.000000e+308 0.000000e+00 0.000000e+00 > roundX(.Machine$double.xmax, -(291:310), "r3d.C")# trace = 0 --> contains 7 x Inf [1] 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 [6] 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 1.797693e+308 [11] 1.797693e+308 1.797693e+308 1.797690e+308 1.797600e+308 1.797000e+308 [16] 1.790000e+308 1.700000e+308 1.000000e+308 0.000000e+00 0.000000e+00 > > proc.time() user system elapsed 0.31 0.03 0.29