R Under development (unstable) (2024-08-21 r87038 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. > ### ----------------- Testing expm1x() { := e^x - 1 - x , numerically stable} ---- > require(DPQ) Loading required package: DPQ > require(sfsmisc) Loading required package: sfsmisc > > > ##------------------------ Testing *accuracy* via Rmpfr {R <--> GNU MPFR C-library} ---- > for(pkg in c("Rmpfr")) + if(!requireNamespace(pkg)) { + cat("no CRAN package", sQuote(pkg), " ---> no tests here.\n") + q("no") + } Loading required namespace: Rmpfr > require("Rmpfr") Loading required package: Rmpfr Loading required package: gmp Attaching package: 'gmp' The following objects are masked from 'package:sfsmisc': factorize, is.whole 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 object is masked from 'package:DPQ': log1mexp 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 > > #------------------------------------------------------- > > do.pdf <- TRUE > do.pdf <- !dev.interactive(orNone = TRUE) > do.pdf [1] TRUE > if(do.pdf) endPdf <- if(interactive()) sfsmisc::end.pdf else dev.off > > ## relErrV <- sfsmisc::relErrV > ## eaxis <- sfsmisc::eaxis > > > ## direct formula - may be really "bad" : > expm1x.0 <- function(x) exp(x) -1 - x > ## less direct formula - improved (but still not universally ok): > expm1x.1 <- function(x) expm1(x) - x > > ## a symmetric set of negative and positive > x <- unique(c(2^-seq(-3/8, 54, by = 1/8), seq(7/8, 3, by = 1/128))) > x <- x0 <- sort(c(-x, 0, x)) # negative *and* positive > > ## Mathematically, expm1x() = exp(x) - 1 - x >= 0 (and == 0 only at x=0): > em1x <- expm1x(x) > stopifnot(em1x >= 0, identical(x == 0, em1x == 0)) > > xxTrue <- expm1x.1(mpfr(x, 1024)) > relE <- asNumeric(relErrV(xxTrue, em1x)) > relE1 <- asNumeric(relErrV(xxTrue, expm1(x)-x)) > > if(do.pdf) pdf("expm1x_relE-1.pdf") > plot(x, abs(relE), log="y", type = "b", cex=1/2, ylim = c(2e-17, max(abs(relE)))) Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > lines(x, abs(relE1), col = adjustcolor(2, 1/2), lwd=3) > abline(h = 2^-(52:51), lty=3, col=paste("gray",c(20,50))) > axis(4, at=2^-52, labels = quote(2^{-52}), col.ticks=NA, las=2, cex.axis=3/4, hadj = 1/2) > > plot (abs(x), abs(relE), log="xy", type = "b", cex=1/2, ylim = c(2e-17, max(abs(relE)))) Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log) : 1 x value <= 0 omitted from logarithmic plot 2: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > lines(abs(x), abs(relE1), col = adjustcolor(2, 1/2), lwd=3) > abline(h = 2^-(52:51), lty=3, col=paste("gray",c(20,50))) > axis(4, at=2^-52, labels = quote(2^{-52}), col.ticks=NA, las=2, cex.axis=3/4, hadj = 1/2) > > iL <- which(abs(relE) > 2^-52) > print(cbind(x = x[iL], relE = relE[iL] , relE1 = relE1[iL]), digits = 4) x relE relE1 [1,] 9.686e-16 -2.255e-16 -1.591e-01 [2,] 5.611e-12 2.321e-16 1.572e-06 [3,] 3.250e-08 -2.773e-16 -5.176e-09 [4,] 3.328e-05 -2.512e-16 2.289e-12 [5,] 1.704e-02 -2.318e-16 3.110e-15 [6,] 1.406e+00 -2.580e-16 -2.580e-16 [7,] 1.438e+00 -2.371e-16 -2.371e-16 > ## x relE relE1 > ## 9.686e-16 -2.255e-16 -1.591e-01 > ## 5.611e-12 2.321e-16 1.572e-06 > ## 3.250e-08 -2.773e-16 -5.176e-09 > ## 3.328e-05 -2.512e-16 2.289e-12 > ## 1.704e-02 -2.318e-16 3.110e-15 > ## 1.211e+00 2.221e-16 2.221e-16 > stopifnot(abs(relE) < 5e-16, sum(abs(relE) > 2^-52) <= 6 + 10)# the above + "some slack" > > > ### -------------------- Relative error of Taylor series approximations ----------- > ## > if(do.pdf && .Device == "pdf") { dev.off(); pdf("expm1x_relE_Taylor.pdf") } > > twoP <- seq(-0.75, 54, by = 1/8) > x <- 2^-twoP > x <- sort(c(-x,x)) # negative *and* positive > e1xAll <- cbind(expm1x.0 = expm1x.0(x), + expm1x.1 = expm1x.1(x), + vapply(1:15, \(k) expm1xTser(x, k=k), x)) > colnames(e1xAll)[-(1:2)] <- paste0("k=",1:15) > > xM <- mpfr(x, 1024) # high accuracy (to push cancellation out of dbl.prec. range) > expm1xM <- expm1x.1(xM) > expm1x.relE <- asNumeric(relErrV(expm1xM, e1xAll)) > > pl.relEexpm1x <- function(x, relE, ind = TRUE, type = "l", las = 2, ..., + leg.x = "top", leg.y = NULL, leg.ncol = 4, + leg.cex = if(.Device == "pdf") 4/5 else 3/4) + { + stopifnot(is.numeric(N <- ncol(relE)), N >= 2) + legs <- c("exp(x) - 1 - x", "expm1(x) - x", + paste0("expm1xTser(x, ", colnames(relE)[-(1:2)], ")")) + stopifnot(N == length(legs)) + ## matplot() default has col = 1:6, lty = 1:5 + col <- rep_len(1:6, N)[ind] + lty <- rep_len(1:5, N)[ind] + matplot(x, relE[, ind], type=type, las=las, xaxt = "n", ..., + col = col, lty = lty, ylab = deparse1(substitute(relE)), + main = expression("Accuracy (Relative Error) of" ~~~~~ e^x - 1 - x ~~ "Computations")) + eaxis(1, sub10 = c(-2, 2)) + legend(leg.x, leg.y, legend = legs[ind], ncol=leg.ncol, + cex = leg.cex, bty = "n", col = col, lty = lty) + } > > ## no x-log scale here at first: > pl.relEexpm1x(x, expm1x.relE, ylim = c(-1,1)*1000, leg.x = "topright", leg.ncol = 3) # non sense > rug(x) > pl.relEexpm1x(x, expm1x.relE, xlim = c(-1,1)*0.5, ylim = c(-1,1)*1000) > pl.relEexpm1x(x, expm1x.relE, xlim = c(-1,1)*0.1, ylim = c(-1,1)*1000) # still "non sense" > > I <- x > 0 # positive > pl.relEexpm1x(x[I], expm1x.relE[I,], log="x", ylim = c(-1,1)*100, leg.x = "topright", leg.ncol = 3) > N <- x < 0 # negative > matlines(-x[N], expm1x.relE[N,]) # not much changes in picture -- good ?! > > ## not showing the worst (= the first): > pl.relEexpm1x(x[I], expm1x.relE[I,], log="x", ylim = c(-1,1) / 2 , ind = -1) > # see how Taylor k=1 and k=2 diverge > ## zoom-in (ylim): > pl.relEexpm1x(x[I], expm1x.relE[I,], log="x", ylim = c(-1,1) /100 , ind = -1) > # see how Taylor k=1,2,3,4 diverge > matlines(-x[N], expm1x.relE[N,-1], col = c(2:6,1), lty = c(2:5,1)) > ## with correct colors -- shows different error *sign* of k= ! > > ## Much more relevant: *relative* error (and log-scale) : > pl.relEexpm1x(x[I], abs(expm1x.relE[I,]), log="xy", yaxt="n", leg.cex = .9) ; eaxis(2, nintLog = 20) > abline(h = 2^-(52:51), lty=3, col=paste("gray",c(20,50))) > axis(4, at=2^-52, labels = quote(2^{-52}), col.ticks=NA, las=2, cex.axis=3/4, hadj = 1/2) > if("nintLog" %in% names(formals(grid))) grid(nintLog = 20) else grid() # in the future R >= 4.5 > matlines(-x[N], expm1x.relE[N,]) # onle little changes in picture -- but visible for "large x" (x ~= 1) good! > ## draw some vertical lines "at" cutpoints x[k]: > abl <- function(v, lty = 3, col = adjustcolor(1, 1/2), lwd = 2, ...) { + ## abline(v = v, lty=lty, col=col, lwd=lwd, ...) + uy <- par("usr")[3:4] + segments(x0 = v, y0 = 10^uy[1], y1 = 2^-52, lty=lty, col=col, lwd=lwd, ...) + if(getOption("scipen") > -2) { op <- options(scipen = -2); on.exit(op) } + form <- function(v) sub("e-0", "e-", format(v, digits=3)) + axis(1, at=v, labels=form(v), padj = -3, cex = 0.75) + } > abl(v = 5.00e-16)# k = 1 <==> ord=2 > abl(v = 4.40e-8) > abl(v = 2.02e-5) # k = 3 <==> ord=4 > abl(v = 3.95e-4) > abl(v = 3.00e-3) # k = 5 <==> ord=6 > ## abl(v = 0.0111) > > if(do.pdf && .Device == "pdf") dev.off() null device 1 > > > ###--------------------- Older Experiments for finding cutoffs and visualization --- originally in ../R/expm1x.R > > (doExtras <- DPQ:::doExtras() && !grepl("valgrind", R.home())) # TRUE, typically when interactive [1] FALSE > if(!doExtras) quit("no") > proc.time() user system elapsed 2.26 0.18 2.42