R Under development (unstable) (2025-07-07 r88391 ucrt) -- "Unsuffered Consequences" Copyright (C) 2025 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 stirlerr() > #### =================== {previous 2nd part of this, now -->>> ./bd0-tst.R <<<--- > require(DPQ) Loading required package: DPQ > for(pkg in c("Rmpfr", "DPQmpfr")) + if(!requireNamespace(pkg)) { + cat("no CRAN package", sQuote(pkg), " ---> no tests here.\n") + q("no") + } Loading required namespace: Rmpfr Loading required namespace: DPQmpfr > 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 object is masked from 'package:DPQ': log1mexp The following objects are masked from 'package:stats': dbinom, dgamma, dnbinom, dnorm, dpois, dt, pgamma, pnorm The following objects are masked from 'package:base': cbind, pmax, pmin, rbind > > source(system.file(package="DPQ", "test-tools.R", mustWork=TRUE)) > ## => showProc.time(), ... list_() , loadList() , readRDS_() , save2RDS() > ##_ options(conflicts.policy = list(depends.ok=TRUE, error=FALSE, warn=FALSE)) > require(sfsmisc) # masking 'list_' *and* gmp's factorize(), is.whole() Loading required package: sfsmisc Attaching package: 'sfsmisc' The following object is masked _by_ '.GlobalEnv': list_ The following objects are masked from 'package:gmp': factorize, is.whole > ##_ options(conflicts.policy = NULL)o > > ## plot1cuts() , etc: ---> ../inst/extraR/relErr-plots.R <<<<<<< > ## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ > source(system.file(package="DPQ", "extraR", "relErr-plots.R", mustWork=TRUE)) > > do.pdf <- TRUE # (manually) > do.pdf <- !dev.interactive(orNone = TRUE) > do.pdf [1] TRUE > if(do.pdf) { + pdf.options(width = 9, height = 6.5) # for all pdf plots {9/6.5 = 1.38 ; 4/3 < 1.38 < sqrt(2) [A4] + pdf("stirlerr-tst.pdf") + } > > showProc.time() Time (user system elapsed): 0.05 0.02 0.06 > (doExtras <- DPQ:::doExtras()) [1] FALSE > (noLdbl <- (.Machine$sizeof.longdouble <= 8)) ## TRUE when --disable-long-double or 'M1mac' .. [1] FALSE > (M.mac <- grepl("aarch64-apple", R.version$platform)) # Mac with M1, M2, ... proc [1] FALSE > > abs19 <- function(r) pmax(abs(r), 1e-19) # cut |err| to positive {for log-plots} > > options(width = 100, nwarnings = 1e5, warnPartialMatchArgs = FALSE) > > > > ##=== Really, dpois_raw() and dbinom_raw() *both* use stirlerr(x) for "all" 'x > 0' > ## ~~~~~~~~~~~ ~~~~~~~~~~~~ =========== ===== ! > > ## below, 6 "it's okay, but *far* from perfect:" ===> need more terms in stirlerr() [ > ## April 20: MM added more terms up to S10; 2024-01: up to S12 ..helps a little only > x <- lseq(1/16, 6, length=2048) > system.time(stM <- DPQmpfr::stirlerrM(Rmpfr::mpfr(x,2048))) # 1.7 sec elapsed user system elapsed 1.95 0.01 1.97 > plot(x, stirlerr(x, use.halves=FALSE) - stM, type="l", log="x", main="absolute Error") > plot(x, stirlerr(x, use.halves=FALSE) / stM - 1, type="l", log="x", main="relative Error") > plot(x, abs(stirlerr(x, use.halves=FALSE) / stM - 1), type="l", log="xy",main="|relative Error|") > drawEps.h(-(52:50)) > ## lgammacor() does *NOT* help, as it is *designed* for x >= 10 ... (but is interesting there) > ## > ## ==> Need another chebyshev() or rational-approx. for x in [.1, 7] or so !! > > ##=============> see also ../Misc/stirlerr-trms.R <=============== > ## ~~~~~~~~~~~~~~~~~~~~~~~ > > cutoffs <- c(15,35,80,500) # cut points, n=*, in the stirlerr() "algorithm" > ## > n <- c(seq(1,15, by=1/4),seq(16, 25, by=1/2), 26:30, seq(32,50, by=2), seq(55,1000, by=5), + 20*c(51:99), 50*(40:80), 150*(27:48), 500*(15:20)) > st.n <- stirlerr(n, "R3")# rather use.halves=TRUE; but here , use.halves=FALSE > plot(st.n ~ n, log="xy", type="b") ## looks good now (straight line descending {in log-log !} > nM <- mpfr(n, 2048) > st.nM <- stirlerr(nM, use.halves=FALSE) ## << on purpose > all.equal(asNumeric(st.nM), st.n)# TRUE [1] TRUE > all.equal(st.nM, as(st.n,"mpfr"))# .. difference: 3.381400e-14 was 1.05884.........e-15 [1] "Mean relative difference: 3.697822e-14" > all.equal(roundMpfr(st.nM, 64), as(st.n,"mpfr"), tolerance=1e-16)# (ditto) [1] "Mean relative difference: 3.697822e-14" > > > ## --- Look at the direct formula -- why is it not good for n ~= 5 ? > ## > ## Preliminary Conclusions : > ## 1. there is *some* cancellation even for small n (how much?) > ## 2. lgamma1p(n) does really not help much compared to lgamma(n+1) --- but a tiny bit in some cases > > ### 1. Investigating lgamma1p(n) vs lgamma(n+1) for n < 1 ============================================= > > ##' @title Relative Error of lgamma(n+1) vs lgamma1p() vs MM's stirlerrD2(): > ##' @param n numeric, typically n << 1 > ##' @param precBits > ##' @return relative error WRT mpfr(n, precBits) > ##' @author Martin Maechler > relE.lgam1 <- function(n, precBits = if(doExtras) 1024 else 320) { + M_LN2PI <- 1.837877066409345483 # ~ log(2*Const("pi",60)); very slightly more accurate than log(2*pi) + st <- lgamma(n +1) - (n +0.5)*log(n) + n - M_LN2PI/2 + st. <- lgamma1p(n) - (n +0.5)*log(n) + n - M_LN2PI/2 # "lgamma1p" + st2 <- lgamma(n) + n*(1-(l.n <- log(n))) + (l.n - M_LN2PI)/2 # "MM2" + st0 <- -(l.n + M_LN2PI)/2 # "n0" + nM <- mpfr(n, precBits) + stM <- lgamma(nM+1) - (nM+0.5)*log(nM) + nM - log(2*Const("pi", precBits))/2 + ## stM <- roundMpfr(stM, 128) + cbind("R3" = asNumeric(relErrV(stM, st)) + , "lgamma1p"= asNumeric(relErrV(stM, st.)) + , "MM2" = asNumeric(relErrV(stM, st2)) + , "n0" = asNumeric(relErrV(stM, st0)) + ) + } > > n <- 2^-seq.int(1022, 1, by = -(if(doExtras) 1/4 else 9/4)) > relEx <- relE.lgam1(n) > showProc.time() Time (user system elapsed): 2.4 0.01 2.42 > > ## Is *equivalent* to 'new' stirlerr_simpl(n version = *) [not for though, see 'relEmat']: > (simpVer <- eval(formals(stirlerr_simpl)$version)) [1] "R3" "lgamma1p" "MM2" "n0" > if(is.null(simpVer)) { warning("got wrong old version of package 'DPQ':") + print(packageDescription("DPQ")) + stop("invalid outdated version package 'DPQ'") + } > stir.allS <- function(n) sapply(simpVer, function(v) stirlerr_simpl(n, version=v)) > stirS <- stir.allS(n) # matrix > nM <- mpfr(n, 256) # "high" precision = 256 should suffice! > stirM <- stirlerr(nM) > releS <- asNumeric(relErrV(stirM, stirS)) > all.equal(relEx, releS, tolerance = 0) # see TRUE on Linux [1] TRUE > stopifnot(all.equal(relEx, releS, tolerance = if(noLdbl) 2e-15 else 1e-15)) > simpVer3 <- simpVer[simpVer != "lgamma1p"] # have no mpfr-ified lgamma1p()! > ## stirlerr_simpl(, *) : > stirM2 <- sapplyMpfr(simpVer3, function(v) stirlerr_simpl(nM, version=v)) > > ## TODO ?: > ## apply(stirM2, 2, function(v) all.equal(v, stirS, check.class=FALSE)) > ## releS2 <- asNumeric(relErrV(stirM2, stirS)) > ## all.equal(relEx, releS, tolerance = 0) # see TRUE on Linux > ## stopifnot(all.equal(relEx, releS, tolerance = 1e-15)) > relEmat <- matrix(NA, ncol(stirM2), ncol(stirS), + dimnames = list(simpVer3, colnames(stirS))) > for(j in seq_len(ncol(stirM2))) + for(k in seq_len(ncol(stirS))) + relEmat[j,k] <- asNumeric(relErr(stirM2[,j], stirS[,k])) > relEmat R3 lgamma1p MM2 n0 R3 5.948059e-17 5.947554e-17 6.015406e-17 6.661766e-06 MM2 5.948059e-17 5.947554e-17 6.015406e-17 6.661766e-06 n0 6.661811e-06 6.661811e-06 6.661811e-06 5.934193e-17 > round(-log10(relEmat), 2) # well .. {why? / expected ?} R3 lgamma1p MM2 n0 R3 16.23 16.23 16.22 5.18 MM2 16.23 16.23 16.22 5.18 n0 5.18 5.18 5.18 16.23 > > cols <- c("gray30", adjustcolor(c(2,3,4), 1/2)); lwd <- c(1, 3,3,3) > stopifnot((k <- length(cols)) == ncol(relEx), k == length(lwd)) > matplot(n, relEx, type = "l", log="x", col=cols, lwd=lwd, ylim = c(-1,1)*4.5e-16, + main = "relative errors of direct (approx.) formula for stirlerr(n), small n") > mtext("really small errors are dominated by small (< 2^-53) errors of log(n)") > ## very interesting: there are different intervals <---> log(n) Qpattern !! > ## -- but very small difference, only for n >~= 1/1000 but not before > drawEps.h(negative=TRUE) # abline(h= c(-4,-2:2, 4)*2^-53, lty=c(2,2,2, 1, 2,2,2), col="gray") > legend("topleft", legend = colnames(relEx), col=cols, lwd=3) > > ## zoomed in a bit: > n. <- 2^-seq.int(400,2, by = -(if(doExtras) 1/4 else 9/4)) > relEx. <- relE.lgam1(n.) > matplot(n., relEx., type = "l", log="x", col=cols, lwd=lwd, ylim = c(-1,1)*4.5e-16, + main = "relative errors of direct (approx.) formula for stirlerr(n), small n") > drawEps.h(negative=TRUE) > legend("topleft", legend = colnames(relEx.), col=cols, lwd=3) > > ##====> Absolute errors (and look at "n0") -------------------------------------- > matplot(n., abs19(relEx.), type = "l", log="xy", col=cols, lwd=lwd, ylim = c(4e-17, 5e-16), + main = quote(abs(relErr(stirlerr_simpl(n, '*'))))) > drawEps.h(); legend("top", legend = colnames(relEx.), col=cols, lwd=3) > lines(n., abs19(relEx.[,"n0"]), type = "o", cex=1/4, col=cols[4], lwd=2) > > ## more zooom-in > n.2 <- 2^-seq.int(85, 50, by= -(if(doExtras) 1/100 else 1/4)) > stirS.2 <- sapply(c("R3", "lgamma1p", "n0"), function(v) stirlerr_simpl(n.2, version=v)) > releS.2 <- asNumeric(relErrV(stirlerr(mpfr(n.2, 320)), stirS.2)) > > matplot(n.2, abs19(releS.2), type = "l", log="xy", col=cols, lwd=lwd, ylim = c(4e-17, 5e-16), + main = quote(abs(relErr(stirlerr_simpl(n, '*'))))) > drawEps.h(); legend("top", legend = colnames(releS.2), col=cols, lwd=3) > abline(v = 5e-17, col=(cb <- adjustcolor("skyblue4", 1/2)), lwd=2, lty=3) > axis(1, at=5e-17, col.axis=cb, line=-1/4, cex = 3/4) > > matplot(n.2, abs19(releS.2), type = "l", log="xy", col=cols, lwd=lwd, ylim = c(4e-17, 5e-16), + xaxt="n", xlim = c(8e-18, 1e-15), ## <<<<<<<<<<<<<<<<<<< Zoom-in + xlab = quote(n), main = quote(abs(relErr(stirlerr_simpl(n, '*'))))) > eaxis(1); drawEps.h(); legend("top", legend = colnames(releS.2), col=cols, lwd=3) > abline(v = 5e-17, col=(cb <- adjustcolor("skyblue4", 1/2)), lwd=2, lty=3) > mtext('stirlerr_simpl(*, "n0") is as good as others for n <= 5e-17', col=adjustcolor(cols[3], 2)) > ## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ > > ## ===> "all *but* "n0" approximations for "larger" small n: > n2 <- 2^-seq.int(20,0, length.out = if(doExtras) 1000 else 200) > relEx2 <- relE.lgam1(n2)[,c("R3", "lgamma1p", "MM2")] # "n0" is "bad": |relE| >= 2.2e-6 ! > cols <- c("gray30", adjustcolor(c(2,4), 1/2)); lwd <- c(1, 3,3) > stopifnot((k <- length(cols)) == ncol(relEx2), k == length(lwd)) > matplot(n2, relEx2, type = "l", log="x", col=cols, lwd=lwd, ylim = c(-3,3)*1e-15, xaxt="n", + main = "relative errors of direct (approx.) formula for stirlerr(n), small n") > eaxis(1, sub10=c(-3,0)); drawEps.h(negative=TRUE) > legend("topleft", legend = colnames(relEx2), col=cols, lwd=3) > ##==> "MM" is *NOT* good for n < 1 *but* > > ## "the same" -- even larger small n: > n3 <- seq(.01, 5, length = if(doExtras) 1000 else 200) > relEx3 <- relE.lgam1(n3)[,c("R3", "lgamma1p", "MM2")] # "no" is "bad" .. > stopifnot((k <- length(cols)) == ncol(relEx3), k == length(lwd)) > > matplot(n3, relEx3, type = "l", col=cols, lwd=lwd, + main = "relative errors of direct (approx.) formula for stirlerr(n), small n") > legend("topleft", legend = colnames(relEx3), col=cols, lwd=3) > drawEps.h(negative=TRUE) > > matplot(n3, abs19(relEx3), type = "l", col=cols, lwd=lwd, + log="y", ylim = 2^-c(54, 44), yaxt = "n", ylab = quote(abs(relE)), xlab=quote(n), + main = "|relative errors| of direct (approx.) formula for stirlerr(n), small n") > eaxis(2, cex.axis=0.9); legend("topleft", legend = colnames(relEx3), col=cols, lwd=3) > drawEps.h() > lines(n3, smooth.spline(abs(relEx3)[,1], df=12)$y, lwd=3, col=cols[1]) > lines(n3, smooth.spline(abs(relEx3)[,2], df=12)$y, lwd=3, col=adjustcolor(cols[2], 1/2)) > lines(n3, smooth.spline(abs(relEx3)[,3], df=12)$y, lwd=4, col=adjustcolor(cols[3], offset = rep(.2,4))) > ## ===> from n >~= 1, "MM2" is definitely better up to n = 5 !! > > ## Check log() only : > plot(n, asNumeric(relErrV(log(mpfr(n, 256)), log(n))), ylim = c(-1,1)*2^-53, + log="x", type="l", xaxt="n") ## ===> indeed --- log(n) approximation pattern !! > eaxis(1) ; drawEps.h(negative=TRUE) > showProc.time() Time (user system elapsed): 0.63 0 0.63 > > ## =========== "R3" vs "lgamma1p" -------------------------- which is better? > > ## really for the very small n, all is dominated by -(n+0.5)*log(n); and lgamma1p() is unnecessary! > i <- 1:20; ni <- n[i] > lgamma1p(ni) [1] -1.284347e-308 -6.109421e-308 -2.906147e-307 -1.382404e-306 -6.575859e-306 -3.128023e-305 [7] -1.487947e-304 -7.077909e-304 -3.366840e-303 -1.601548e-302 -7.618289e-302 -3.623889e-301 [13] -1.723822e-300 -8.199926e-300 -3.900564e-299 -1.855431e-298 -8.825969e-298 -4.198362e-297 [19] -1.997089e-296 -9.499809e-296 > - (ni +0.5)*log(ni) + ni [1] 354.1982 353.4184 352.6386 351.8588 351.0790 350.2993 349.5195 348.7397 347.9599 347.1801 [11] 346.4003 345.6205 344.8407 344.0609 343.2811 342.5014 341.7216 340.9418 340.1620 339.3822 > > ## much less extreme: > n2 <- lseq(2^-12, 1/2, length = if(doExtras) 1000 else 200) > relE2 <- relE.lgam1(n2)[,-4] > > cols <- c("gray30", adjustcolor(2:3, 1/2)); lwd <- c(1,3,3) > matplot(n2, relE2, type = "l", log="x", col=cols, lwd=lwd) > legend("topleft", legend=colnames(relE2), col=cols, lwd=2, lty=1:3) > drawEps.h(negative=TRUE) > > matplot(n2, abs19(relE2), type = "l", log="xy", col=cols, lwd=lwd, ylim = c(6e-17, 1e-15), + xaxt = "n"); eaxis(1, sub10=c(-2,0)) > legend("topleft", legend=colnames(relE2), col=cols, lwd=2, lty=1:3) > drawEps.h() > ## "MM2" is *worse* here, n < 1/2 > for(j in 1:3) lines(n2, smooth.spline(abs(relE2[,j]), df=10)$y, lwd=3, + col=adjustcolor(cols[j], 1.5, offset = rep(-1/4, 4))) > ## "lgammap very slightly better in [0.002, 0.05] ... > ## "TODO": draw 0.90-quantile curves {--> cobs::cobs() ?} instead of mean-curves? > > ## which is better? ... "random difference" > d.absrelE <- abs(relE2[,"R3"]) - abs(relE2[,"lgamma1p"]) > plot (n2, d.absrelE, type = "l", log="x", # no clear picture ... + main = "|relE_R3| - |relE_lgamma1p|", axes=FALSE, frame.plot=TRUE) > eaxis(1, sub10=c(-2,1)); eaxis(2); axis(3, at=max(n2)); abline(v = max(n2), lty=3, col="gray") > ## 'lgamma1p' very slightly better: > lines(n2, smooth.spline(d.absrelE, df=12)$y, lwd=3, col=2) > > ## not really small n at all == here see, how "bad" the direct formula gets for 1 < n < 10 or so > n3 <- lseq(2^-14, 2^2, length = if(doExtras) 800 else 150) > relE3 <- relE.lgam1(n3)[, -4] > > matplot(n3, relE3, type = "l", log="x", col=cols, lty=1, lwd = c(1,3), + main = quote(rel.lgam1(n)), xlab=quote(n)) > > matplot(n3, abs19(relE3), type = "l", log="xy", col=cols, lwd = c(1,3), xaxt="n", + main = quote(abs(rel.lgam1(n))), xlab=quote(n), ylim = c(2e-17, 4e-14)) > drawEps.h(); eaxis(1, sub10=c(-2,3)) > legend("topleft", legend=colnames(relE3), col=cols, lwd=2) > ## very small difference --- draw the 3 smoothers : > for(j in 1:3) { + ll <- lowess(log(n3), abs19(relE3[, j]), f= 1/12) + with(ll, lines(exp(x), y, col=adjustcolor(cols[j], 1.5), lwd=3)) + } > ## ==> lgamma1p(.) very slightly in n ~ 10^-4 -- 10^-2 --- but not where it matters: n ~ 0.1 -- 1 !! > ## "MM2" gets best from n >~ 1 ! > abline(v=1, lty=3, col = adjustcolor(1, 3/4)) > showProc.time() Time (user system elapsed): 0.12 0 0.12 > > > ### 2. relErr( stirlerr(.) ) ============================================================ > > ##' Very revealing plot showing the *relative* approximation error of stirlerr() > ##' > p.stirlerrDev <- function(n, precBits = if(doExtras) 2048L else 512L, + stnM = stirlerr(mpfr(n, precBits), use.halves=use.halves, verbose=verbose), + abs = FALSE, + ## cut points, n=*, in the stirlerr() algorithm; "FIXME": sync with ../R/dgamma.R <<<< + scheme = c("R3", "R4.4_0"), + cutoffs = switch(match.arg(scheme) + , R3 = c(15, 35, 80, 500) + , R4.4_0 = c(4.9, 5.0, 5.1, 5.2, 5.3, 5.4, 5.7, + 6.1, 6.5, 7, 7.9, 8.75, 10.5, 13, + 20, 26, 60, 200, 3300, 17.4e6) + ## {FIXME: need to sync} <==> ../man/stirlerr.Rd <==> ../R/dgamma.R + ), + use.halves = missing(cutoffs), + direct.ver = c("R3", "lgamma1p", "MM2", "n0"), + verbose = getOption("verbose"), + type = "b", cex = 1, + col = adjustcolor(1, 3/4), colnB = adjustcolor("orange4", 1/3), + log = if(abs) "xy" else "x", + xlim=NULL, ylim = if(abs) c(8e-18, max(abs(N(relE))))) + { + op <- par(las = 1, mgp=c(2, 0.6, 0)) + on.exit(par(op)) + require("Rmpfr"); require("sfsmisc") + st <- stirlerr(n, scheme=scheme, cutoffs=cutoffs, use.halves=use.halves, direct.ver=direct.ver, + verbose=verbose) + relE <- relErrV(stnM, st) # eps0 = .Machine$double.xmin + N <- asNumeric + form <- if(abs) abs(N(relE)) ~ n else N(relE) ~ n + plot(form, log=log, type=type, cex=cex, col=col, xlim=xlim, ylim=ylim, + ylab = quote(relErrV(stM, st)), axes=FALSE, frame.plot=TRUE, + main = sprintf("stirlerr(n, cutoffs) rel.error [wrt stirlerr(Rmpfr::mpfr(n, %d))]", + precBits)) + eaxis(1, sub10=3) + eaxis(2) + mtext(paste("cutoffs =", deparse1(cutoffs))) + ylog <- par("ylog") + ## FIXME: improve this ---> drawEps.h() above + if(ylog) { + epsC <- c(1,2,4,8)*2^-52 + epsCxp <- expression(epsilon[C],2*epsilon[C], 4*epsilon[C], 8*epsilon[C]) + } else { + epsC <- (-2:2)*2^-52 + epsCxp <- expression(-2*epsilon[C],-epsilon[C], 0, +epsilon[C], +2*epsilon[C]) + } + dy <- diff(par("usr")[3:4]) + if(diff(range(if(ylog) log10(epsC) else epsC)) > dy/50) { + lw <- rep(1/2, 5); lw[if(ylog) 1 else 3] <- 2 + abline( h=epsC, lty=3, lwd=lw) + axis(4, at=epsC, epsCxp, las=2, cex.axis = 3/4, mgp=c(3/4, 1/4, 0), tck=0) + } else ## only x-axis + abline(h=if(ylog) epsC else 0, lty=3, lwd=2) + abline(v = cutoffs, col=colnB) + axis(3, at=cutoffs, col=colnB, col.axis=colnB, + labels = formatC(cutoffs, digits=3, width=1)) + invisible(relE) + } ## p.stirlerrDev() > > showProc.time() Time (user system elapsed): 0.02 0 0.02 > > n <- lseq(2^-10, 1e10, length = if(doExtras) 4096 else 512) > n <- lseq(2^-10, 5000, length = if(doExtras) 4096 else 512) > ## store "expensive" stirlerr() result, and re-use many times below: > nM <- mpfr(n, if(doExtras) 2048 else 512) > st.nM <- stirlerr(nM, use.halves=FALSE) ## << on purpose > > p.stirlerrDev(n=n, stnM=st.nM, use.halves = FALSE) # default cutoffs= c(15, 40, 85, 600) > p.stirlerrDev(n=n, stnM=st.nM, use.halves = FALSE, ylim = c(-1,1)*1e-12) # default cutoffs= c(15, 40, 85, 600) > > ## show the zoom-in region in next plot > yl2 <- 3e-14*c(-1,1) > abline(h = yl2, col=adjustcolor("tomato", 1/4), lwd=3, lty=2) > > if(do.pdf) { dev.off() ; pdf("stirlerr-relErr_1.pdf") } > > ## drop n < 7: > p.stirlerrDev(n=n, stnM=st.nM, xlim = c(7, max(n)), use.halves=FALSE) # default cutoffs= c(15, 40, 85, 600) > abline(h = yl2, col=adjustcolor("tomato", 1/4), lwd=3, lty=2) > > ## The first plot clearly shows we should do better: > ## Current code is switching to less terms too early, loosing up to 2 decimals precision > if(FALSE) # no visible difference {use.halves = T / F }: + p.stirlerrDev(n=n, stnM=st.nM, ylim = yl2, use.halves = FALSE) > p.stirlerrDev(n=n, stnM=st.nM, ylim = yl2, use.halves = TRUE)# exact at n/2 (n <= ..) > abline(h = yl2, col=adjustcolor("tomato", 1/4), lwd=3, lty=2) > > showProc.time() Time (user system elapsed): 0.22 0.02 0.23 > > > if(do.pdf) { dev.off(); pdf("stirlerr-relErr_6-fin-1.pdf") } > > ### ~19.April 2021: "This is close to *the* solution" (but see 'cuts' below) > cuts <- c(7, 12, 20, 26, 60, 200, 3300) > ## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ > st. <- stirlerr(n=n , cutoffs = cuts, verbose=TRUE) stirlerr(n, cutoffs = 7,12,20,26,60,200,3300) : case I (n <= 7), using direct formula for n= num [1:294] 0.000977 0.001007 0.001037 0.001069 0.001102 ... case II (n > 7 ), 7 cutoffs: ( 7, 12, 20, 26, 60, 200, 3300 ): n in cutoff intervals: (7,12] (12,20] (20,26] (26,60] (60,200] (200,3.3e+03] (3.3e+03,Inf] 18 17 9 27 40 93 14 > st.nM <- stirlerr(n=nM, cutoffs = cuts, use.halves=FALSE) ## << on purpose > relE <- asNumeric(relErrV(st.nM, st.)) > head(cbind(n, relE), 20) n relE [1,] 0.0009765625 -5.082434e-18 [2,] 0.0010065370 -4.617853e-17 [3,] 0.0010374315 8.200891e-17 [4,] 0.0010692742 1.764253e-16 [5,] 0.0011020944 5.113875e-17 [6,] 0.0011359219 -5.412941e-17 [7,] 0.0011707877 -1.196703e-16 [8,] 0.0012067237 -1.165306e-16 [9,] 0.0012437627 1.116723e-16 [10,] 0.0012819386 -9.937293e-18 [11,] 0.0013212862 -1.852831e-16 [12,] 0.0013618416 -5.152507e-17 [13,] 0.0014036418 -2.835242e-18 [14,] 0.0014467250 3.078313e-17 [15,] 0.0014911305 -2.388434e-16 [16,] 0.0015368991 1.367014e-16 [17,] 0.0015840724 -7.481361e-17 [18,] 0.0016326937 1.315003e-16 [19,] 0.0016828074 9.236353e-17 [20,] 0.0017344592 -2.455704e-17 > ## nice printout : > print(cbind(n = format(n, drop0trailing = TRUE), + stirlerr= format(st.,scientific=FALSE, digits=4), + relErr = signif(relE, 4)) + , quote=FALSE) n stirlerr relErr [1,] 9.765625e-04 2.55398004 -5.082e-18 [2,] 1.006537e-03 2.53905399 -4.618e-17 [3,] 1.037431e-03 2.52413284 8.201e-17 [4,] 1.069274e-03 2.50921671 1.764e-16 [5,] 1.102094e-03 2.49430574 5.114e-17 [6,] 1.135922e-03 2.47940003 -5.413e-17 [7,] 1.170788e-03 2.46449973 -1.197e-16 [8,] 1.206724e-03 2.44960497 -1.165e-16 [9,] 1.243763e-03 2.43471589 1.117e-16 [10,] 1.281939e-03 2.41983263 -9.937e-18 [11,] 1.321286e-03 2.40495532 -1.853e-16 [12,] 1.361842e-03 2.39008412 -5.153e-17 [13,] 1.403642e-03 2.37521918 -2.835e-18 [14,] 1.446725e-03 2.36036065 3.078e-17 [15,] 1.491131e-03 2.34550868 -2.388e-16 [16,] 1.536899e-03 2.33066344 1.367e-16 [17,] 1.584072e-03 2.31582510 -7.481e-17 [18,] 1.632694e-03 2.30099381 1.315e-16 [19,] 1.682807e-03 2.28616975 9.236e-17 [20,] 1.734459e-03 2.27135310 -2.456e-17 [21,] 1.787696e-03 2.25654404 1.333e-16 [22,] 1.842568e-03 2.24174275 9.869e-17 [23,] 1.899123e-03 2.22694942 1.75e-16 [24,] 1.957415e-03 2.21216423 2.298e-17 [25,] 2.017495e-03 2.19738740 1.539e-16 [26,] 2.079420e-03 2.18261911 -5.538e-17 [27,] 2.143245e-03 2.16785957 -8.089e-18 [28,] 2.209030e-03 2.15310900 1.456e-16 [29,] 2.276834e-03 2.13836760 2.27e-16 [30,] 2.346718e-03 2.12363561 7.185e-17 [31,] 2.418748e-03 2.10891323 -7.158e-17 [32,] 2.492989e-03 2.09420070 -7.567e-17 [33,] 2.569508e-03 2.07949825 -2.601e-16 [34,] 2.648376e-03 2.06480613 -7.146e-17 [35,] 2.729665e-03 2.05012458 -1.743e-16 [36,] 2.813449e-03 2.03545384 -5.568e-17 [37,] 2.899805e-03 2.02079417 1.915e-16 [38,] 2.988811e-03 2.00614584 -3.373e-16 [39,] 3.080549e-03 1.99150910 1.422e-16 [40,] 3.175103e-03 1.97688423 1.416e-16 [41,] 3.272559e-03 1.96227151 2.276e-16 [42,] 3.373007e-03 1.94767122 5.332e-17 [43,] 3.476537e-03 1.93308365 5.797e-17 [44,] 3.583246e-03 1.91850910 8.484e-17 [45,] 3.693229e-03 1.90394786 -1.687e-16 [46,] 3.806589e-03 1.88940024 -2.104e-16 [47,] 3.923428e-03 1.87486657 -1.251e-16 [48,] 4.043853e-03 1.86034716 -9.698e-18 [49,] 4.167974e-03 1.84584233 -3.746e-17 [50,] 4.295905e-03 1.83135242 6.331e-17 [51,] 4.427763e-03 1.81687778 2.45e-16 [52,] 4.563668e-03 1.80241875 3.129e-16 [53,] 4.703745e-03 1.78797568 1.701e-16 [54,] 4.848121e-03 1.77354894 2.69e-17 [55,] 4.996928e-03 1.75913889 -1.562e-16 [56,] 5.150303e-03 1.74474592 1.094e-17 [57,] 5.308386e-03 1.73037040 5.266e-17 [58,] 5.471321e-03 1.71601274 -2.156e-16 [59,] 5.639257e-03 1.70167332 -4.997e-17 [60,] 5.812347e-03 1.68735256 -1.317e-16 [61,] 5.990751e-03 1.67305086 3.599e-16 [62,] 6.174630e-03 1.65876866 2.062e-16 [63,] 6.364153e-03 1.64450638 3.104e-16 [64,] 6.559494e-03 1.63026447 2.961e-16 [65,] 6.760830e-03 1.61604336 1.711e-16 [66,] 6.968346e-03 1.60184353 4.538e-17 [67,] 7.182231e-03 1.58766542 -5.506e-18 [68,] 7.402681e-03 1.57350952 3.318e-16 [69,] 7.629898e-03 1.55937630 -1.015e-16 [70,] 7.864089e-03 1.54526626 1.093e-16 [71,] 8.105468e-03 1.53117990 1.838e-16 [72,] 8.354257e-03 1.51711772 2.27e-16 [73,] 8.610681e-03 1.50308025 -1.823e-16 [74,] 8.874976e-03 1.48906800 2.763e-16 [75,] 9.147383e-03 1.47508152 -5.773e-17 [76,] 9.428152e-03 1.46112135 2.12e-16 [77,] 9.717538e-03 1.44718804 1.786e-16 [78,] 1.001581e-02 1.43328216 4.22e-16 [79,] 1.032323e-02 1.41940429 2.179e-17 [80,] 1.064009e-02 1.40555500 6.187e-16 [81,] 1.096668e-02 1.39173490 -1.428e-16 [82,] 1.130328e-02 1.37794458 1.669e-16 [83,] 1.165023e-02 1.36418465 5.435e-17 [84,] 1.200782e-02 1.35045575 -1.93e-16 [85,] 1.237638e-02 1.33675849 -3.173e-16 [86,] 1.275626e-02 1.32309354 -5.002e-17 [87,] 1.314780e-02 1.30946153 -1.318e-17 [88,] 1.355136e-02 1.29586313 1.268e-17 [89,] 1.396730e-02 1.28229901 2.717e-17 [90,] 1.439601e-02 1.26876987 9.635e-17 [91,] 1.483788e-02 1.25527638 -6.562e-17 [92,] 1.529331e-02 1.24181926 5.204e-17 [93,] 1.576272e-02 1.22839922 1.471e-18 [94,] 1.624654e-02 1.21501698 -3.29e-16 [95,] 1.674521e-02 1.20167328 1.943e-16 [96,] 1.725919e-02 1.18836886 1.32e-16 [97,] 1.778894e-02 1.17510448 -4.49e-16 [98,] 1.833495e-02 1.16188089 2.876e-16 [99,] 1.889772e-02 1.14869888 7.296e-16 [100,] 1.947776e-02 1.13555923 3.49e-17 [101,] 2.007561e-02 1.12246274 -2.129e-16 [102,] 2.069181e-02 1.10941020 3.364e-16 [103,] 2.132692e-02 1.09640244 -1.391e-16 [104,] 2.198152e-02 1.08344027 1.911e-16 [105,] 2.265622e-02 1.07052453 7.137e-17 [106,] 2.335163e-02 1.05765606 -9.943e-17 [107,] 2.406838e-02 1.04483572 5.531e-17 [108,] 2.480713e-02 1.03206437 -4.731e-16 [109,] 2.556856e-02 1.01934287 8.037e-17 [110,] 2.635335e-02 1.00667211 -8.774e-17 [111,] 2.716224e-02 0.99405297 -4.024e-17 [112,] 2.799595e-02 0.98148635 3.736e-16 [113,] 2.885526e-02 0.96897315 -9.534e-17 [114,] 2.974094e-02 0.95651430 1.561e-17 [115,] 3.065380e-02 0.94411070 8.567e-17 [116,] 3.159468e-02 0.93176328 1.156e-16 [117,] 3.256445e-02 0.91947298 -3.301e-16 [118,] 3.356398e-02 0.90724075 3.528e-18 [119,] 3.459418e-02 0.89506753 -1.274e-16 [120,] 3.565601e-02 0.88295427 2.864e-17 [121,] 3.675043e-02 0.87090194 -7.201e-17 [122,] 3.787845e-02 0.85891151 1.682e-16 [123,] 3.904108e-02 0.84698394 -3.473e-16 [124,] 4.023940e-02 0.83512021 -2.394e-16 [125,] 4.147450e-02 0.82332130 2.555e-18 [126,] 4.274752e-02 0.81158820 3.882e-16 [127,] 4.405960e-02 0.79992190 6.649e-17 [128,] 4.541196e-02 0.78832338 -2.048e-16 [129,] 4.680583e-02 0.77679364 2.321e-17 [130,] 4.824248e-02 0.76533368 -3.847e-16 [131,] 4.972323e-02 0.75394449 -1.664e-16 [132,] 5.124943e-02 0.74262707 2.521e-17 [133,] 5.282247e-02 0.73138242 -2.333e-16 [134,] 5.444379e-02 0.72021153 -1.738e-17 [135,] 5.611488e-02 0.70911541 1.059e-16 [136,] 5.783727e-02 0.69809505 9.541e-17 [137,] 5.961252e-02 0.68715144 -1.962e-16 [138,] 6.144225e-02 0.67628557 -7.252e-17 [139,] 6.332815e-02 0.66549842 -2.304e-16 [140,] 6.527194e-02 0.65479099 1.205e-16 [141,] 6.727539e-02 0.64416424 9.112e-17 [142,] 6.934033e-02 0.63361914 2.165e-16 [143,] 7.146865e-02 0.62315666 -6.79e-17 [144,] 7.366230e-02 0.61277776 -2.332e-16 [145,] 7.592328e-02 0.60248338 -9.743e-17 [146,] 7.825365e-02 0.59227445 3.514e-16 [147,] 8.065556e-02 0.58215192 1.902e-16 [148,] 8.313119e-02 0.57211669 -5.762e-16 [149,] 8.568281e-02 0.56216966 -4.422e-16 [150,] 8.831274e-02 0.55231174 -2.864e-16 [151,] 9.102340e-02 0.54254380 -2.739e-16 [152,] 9.381726e-02 0.53286670 1.238e-16 [153,] 9.669687e-02 0.52328129 -2.44e-16 [154,] 9.966487e-02 0.51378840 -5.457e-16 [155,] 1.027240e-01 0.50438884 -8.706e-19 [156,] 1.058770e-01 0.49508341 -3.148e-16 [157,] 1.091267e-01 0.48587288 -9.274e-17 [158,] 1.124763e-01 0.47675801 -1.941e-16 [159,] 1.159286e-01 0.46773953 8.576e-17 [160,] 1.194869e-01 0.45881815 2.611e-17 [161,] 1.231544e-01 0.44999455 2.191e-16 [162,] 1.269345e-01 0.44126939 -2.471e-17 [163,] 1.308306e-01 0.43264332 -4.582e-16 [164,] 1.348463e-01 0.42411694 -8.929e-16 [165,] 1.389852e-01 0.41569083 2.523e-16 [166,] 1.432512e-01 0.40736554 -3.684e-16 [167,] 1.476482e-01 0.39914160 2.227e-16 [168,] 1.521801e-01 0.39101950 2.247e-16 [169,] 1.568510e-01 0.38299970 -1.866e-16 [170,] 1.616654e-01 0.37508263 5.162e-16 [171,] 1.666275e-01 0.36726868 -4.9e-16 [172,] 1.717420e-01 0.35955822 4.47e-16 [173,] 1.770134e-01 0.35195157 -4.29e-16 [174,] 1.824466e-01 0.34444902 -2.053e-16 [175,] 1.880466e-01 0.33705084 -6.914e-16 [176,] 1.938185e-01 0.32975722 7.402e-16 [177,] 1.997675e-01 0.32256836 -2.994e-16 [178,] 2.058992e-01 0.31548440 -6.361e-16 [179,] 2.122190e-01 0.30850543 -6.348e-16 [180,] 2.187328e-01 0.30163152 -2.111e-16 [181,] 2.254466e-01 0.29486269 -6.115e-16 [182,] 2.323664e-01 0.28819892 2.833e-16 [183,] 2.394986e-01 0.28164016 -7.802e-17 [184,] 2.468498e-01 0.27518629 2.394e-16 [185,] 2.544265e-01 0.26883719 1.048e-16 [186,] 2.622359e-01 0.26259266 5.306e-16 [187,] 2.702849e-01 0.25645248 -4.927e-16 [188,] 2.785810e-01 0.25041639 -7.581e-16 [189,] 2.871317e-01 0.24448407 1.159e-15 [190,] 2.959449e-01 0.23865516 -7.609e-16 [191,] 3.050286e-01 0.23292929 3.061e-16 [192,] 3.143911e-01 0.22730600 -3.57e-16 [193,] 3.240410e-01 0.22178483 -3.899e-16 [194,] 3.339870e-01 0.21636526 -8.041e-16 [195,] 3.442384e-01 0.21104673 -4.607e-16 [196,] 3.548044e-01 0.20582863 -2.497e-16 [197,] 3.656947e-01 0.20071034 4.085e-16 [198,] 3.769193e-01 0.19569118 -6.634e-16 [199,] 3.884884e-01 0.19077043 -1.255e-16 [200,] 4.004126e-01 0.18594734 6.421e-16 [201,] 4.127028e-01 0.18122112 -1.109e-15 [202,] 4.253702e-01 0.17659094 3.218e-16 [203,] 4.384265e-01 0.17205595 2.142e-16 [204,] 4.518835e-01 0.16761526 -3.358e-16 [205,] 4.657535e-01 0.16326792 1.69e-17 [206,] 4.800493e-01 0.15901300 1.246e-16 [207,] 4.947839e-01 0.15484949 6.442e-16 [208,] 5.099707e-01 0.15077639 4.607e-16 [209,] 5.256237e-01 0.14679264 -8.463e-16 [210,] 5.417571e-01 0.14289718 -3.454e-17 [211,] 5.583857e-01 0.13908891 8.42e-16 [212,] 5.755247e-01 0.13536671 -2.139e-15 [213,] 5.931898e-01 0.13172943 1.03e-15 [214,] 6.113970e-01 0.12817591 -1.075e-15 [215,] 6.301632e-01 0.12470497 -2.034e-16 [216,] 6.495053e-01 0.12131541 2.74e-16 [217,] 6.694411e-01 0.11800601 6.725e-16 [218,] 6.899889e-01 0.11477554 -1.83e-16 [219,] 7.111673e-01 0.11162275 1.027e-16 [220,] 7.329957e-01 0.10854638 1.307e-15 [221,] 7.554942e-01 0.10554517 -3.429e-16 [222,] 7.786832e-01 0.10261785 -4.672e-16 [223,] 8.025840e-01 0.09976312 -1.532e-15 [224,] 8.272184e-01 0.09697970 -1.226e-15 [225,] 8.526090e-01 0.09426630 1.038e-15 [226,] 8.787788e-01 0.09162162 -1.372e-15 [227,] 9.057519e-01 0.08904436 -8.952e-16 [228,] 9.335529e-01 0.08653323 1.457e-15 [229,] 9.622073e-01 0.08408693 -2.146e-16 [230,] 9.917411e-01 0.08170416 5.857e-16 [231,] 1.022181 0.07938364 -7.016e-16 [232,] 1.053556 0.07712409 2.615e-15 [233,] 1.085894 0.07492422 -1.481e-15 [234,] 1.119224 0.07278276 9.102e-17 [235,] 1.153577 0.07069846 6.244e-16 [236,] 1.188985 0.06867005 -1.782e-15 [237,] 1.225480 0.06669631 -1.084e-15 [238,] 1.263094 0.06477599 -1.079e-15 [239,] 1.301864 0.06290789 1.556e-15 [240,] 1.341823 0.06109079 -6.892e-16 [241,] 1.383009 0.05932350 9.288e-16 [242,] 1.425458 0.05760485 2.242e-15 [243,] 1.469211 0.05593369 7.296e-17 [244,] 1.514307 0.05430885 -3.486e-15 [245,] 1.560787 0.05272922 -2.088e-15 [246,] 1.608694 0.05119368 -1.686e-15 [247,] 1.658071 0.04970113 2.864e-15 [248,] 1.708963 0.04825051 -5.645e-15 [249,] 1.761418 0.04684075 3.899e-16 [250,] 1.815482 0.04547082 1.06e-15 [251,] 1.871207 0.04413968 -3.861e-16 [252,] 1.928641 0.04284634 -6.73e-15 [253,] 1.987839 0.04158982 -1.919e-15 [254,] 2.048853 0.04036914 -1.569e-15 [255,] 2.111740 0.03918337 7.659e-15 [256,] 2.176558 0.03803158 -3.946e-15 [257,] 2.243365 0.03691285 8.329e-15 [258,] 2.312222 0.03582631 -1.107e-14 [259,] 2.383193 0.03477108 -1.046e-14 [260,] 2.456343 0.03374632 -1.528e-14 [261,] 2.531737 0.03275119 -5.937e-15 [262,] 2.609446 0.03178488 -1.232e-14 [263,] 2.689540 0.03084660 1.045e-15 [264,] 2.772092 0.02993557 -5.673e-15 [265,] 2.857178 0.02905105 7.394e-15 [266,] 2.944876 0.02819228 5.523e-15 [267,] 3.035266 0.02735856 -8.791e-15 [268,] 3.128430 0.02654918 -6.932e-15 [269,] 3.224453 0.02576345 1.525e-14 [270,] 3.323424 0.02500071 2.953e-14 [271,] 3.425433 0.02426030 7.679e-15 [272,] 3.530573 0.02354160 6.762e-15 [273,] 3.638940 0.02284398 2.706e-15 [274,] 3.750633 0.02216684 -1.86e-14 [275,] 3.865754 0.02150961 3.216e-14 [276,] 3.984409 0.02087170 4.841e-15 [277,] 4.106706 0.02025256 3.919e-14 [278,] 4.232757 0.01965165 -3.739e-14 [279,] 4.362676 0.01906846 2.377e-14 [280,] 4.496584 0.01850245 4.483e-14 [281,] 4.634601 0.01795315 -3.642e-14 [282,] 4.776855 0.01742006 -5.002e-14 [283,] 4.923475 0.01690271 -5.425e-14 [284,] 5.074595 0.01640064 -8.572e-14 [285,] 5.230354 0.01591342 -1.773e-14 [286,] 5.390894 0.01544061 7.777e-14 [287,] 5.556361 0.01498178 -2.167e-14 [288,] 5.726907 0.01453653 4.676e-14 [289,] 5.902688 0.01410446 -7.586e-14 [290,] 6.083864 0.01368519 -2.034e-14 [291,] 6.270602 0.01327834 -4.127e-14 [292,] 6.463071 0.01288355 -4.094e-16 [293,] 6.661447 0.01250046 8.344e-14 [294,] 6.865913 0.01212872 1.011e-14 [295,] 7.076654 0.01176802 -4.737e-14 [296,] 7.293864 0.01141801 -2.954e-14 [297,] 7.517741 0.01107839 -1.833e-14 [298,] 7.748489 0.01074884 -1.131e-14 [299,] 7.986320 0.01042908 -7.051e-15 [300,] 8.231451 0.01011881 -4.541e-15 [301,] 8.484106 0.00981776 -2.781e-15 [302,] 8.744516 0.00952564 -1.699e-15 [303,] 9.012919 0.00924221 -1.042e-15 [304,] 9.289560 0.00896719 -6.288e-16 [305,] 9.574692 0.00870034 -3.56e-16 [306,] 9.868577 0.00844143 -3.624e-16 [307,] 1.017148e+01 0.00819021 -3.341e-16 [308,] 1.048368e+01 0.00794646 -1.412e-16 [309,] 1.080547e+01 0.00770995 -2.076e-16 [310,] 1.113713e+01 0.00748047 -9.087e-17 [311,] 1.147897e+01 0.00725782 -1.937e-16 [312,] 1.183130e+01 0.00704179 -2.157e-17 [313,] 1.219445e+01 0.00683218 2.671e-16 [314,] 1.256875e+01 0.00662881 2.99e-17 [315,] 1.295453e+01 0.00643148 1.162e-16 [316,] 1.335216e+01 0.00624002 6.775e-17 [317,] 1.376198e+01 0.00605426 2.234e-17 [318,] 1.418439e+01 0.00587403 -9.246e-17 [319,] 1.461977e+01 0.00569916 -4.305e-17 [320,] 1.506850e+01 0.00552949 -5.258e-17 [321,] 1.553101e+01 0.00536487 -8.654e-17 [322,] 1.600772e+01 0.00520514 -6.793e-17 [323,] 1.649906e+01 0.00505018 -9.834e-17 [324,] 1.700548e+01 0.00489982 -1.841e-16 [325,] 1.752744e+01 0.00475393 -1.561e-16 [326,] 1.806543e+01 0.00461239 4.877e-17 [327,] 1.861993e+01 0.00447506 -8.2e-18 [328,] 1.919144e+01 0.00434182 -9.85e-17 [329,] 1.978050e+01 0.00421254 -3.257e-17 [330,] 2.038764e+01 0.00408712 -4.636e-17 [331,] 2.101342e+01 0.00396542 4.236e-17 [332,] 2.165840e+01 0.00384735 -1.606e-16 [333,] 2.232318e+01 0.00373279 -9.668e-18 [334,] 2.300836e+01 0.00362164 -8.364e-17 [335,] 2.371458e+01 0.00351380 -1.504e-16 [336,] 2.444247e+01 0.00340918 5.717e-17 [337,] 2.519271e+01 0.00330766 -1.036e-16 [338,] 2.596597e+01 0.00320917 -8.276e-17 [339,] 2.676296e+01 0.00311361 5.781e-17 [340,] 2.758442e+01 0.00302090 7.786e-17 [341,] 2.843109e+01 0.00293094 1.128e-16 [342,] 2.930375e+01 0.00284367 5.834e-17 [343,] 3.020320e+01 0.00275899 1.354e-17 [344,] 3.113025e+01 0.00267683 -2.765e-17 [345,] 3.208576e+01 0.00259712 -7.872e-19 [346,] 3.307059e+01 0.00251979 -5.821e-17 [347,] 3.408566e+01 0.00244475 -6.499e-17 [348,] 3.513188e+01 0.00237195 1.016e-16 [349,] 3.621021e+01 0.00230132 -2.027e-16 [350,] 3.732164e+01 0.00223279 4.664e-17 [351,] 3.846719e+01 0.00216630 -1.874e-17 [352,] 3.964789e+01 0.00210179 -6.279e-17 [353,] 4.086484e+01 0.00203920 -2.186e-18 [354,] 4.211914e+01 0.00197848 -1.886e-16 [355,] 4.341194e+01 0.00191956 1.02e-17 [356,] 4.474442e+01 0.00186240 -6.62e-17 [357,] 4.611780e+01 0.00180694 -3.302e-17 [358,] 4.753333e+01 0.00175313 -1.145e-16 [359,] 4.899231e+01 0.00170092 -5.072e-17 [360,] 5.049607e+01 0.00165027 -2.727e-17 [361,] 5.204599e+01 0.00160113 -6.057e-17 [362,] 5.364348e+01 0.00155345 -1.004e-16 [363,] 5.529001e+01 0.00150719 -8.365e-17 [364,] 5.698707e+01 0.00146230 -2.147e-17 [365,] 5.873623e+01 0.00141876 -9.664e-17 [366,] 6.053907e+01 0.00137651 -9.109e-17 [367,] 6.239725e+01 0.00133552 -1.26e-16 [368,] 6.431246e+01 0.00129575 -9.707e-17 [369,] 6.628645e+01 0.00125716 -2.326e-17 [370,] 6.832104e+01 0.00121972 -7.206e-17 [371,] 7.041808e+01 0.00118340 -1.019e-16 [372,] 7.257948e+01 0.00114816 -1.926e-16 [373,] 7.480722e+01 0.00111397 -1.273e-16 [374,] 7.710335e+01 0.00108079 1.044e-17 [375,] 7.946994e+01 0.00104861 -2.502e-17 [376,] 8.190918e+01 0.00101738 6.244e-17 [377,] 8.442329e+01 0.00098708 -8.512e-17 [378,] 8.701457e+01 0.00095769 1.691e-17 [379,] 8.968538e+01 0.00092917 8.207e-18 [380,] 9.243817e+01 0.00090150 -4.955e-17 [381,] 9.527546e+01 0.00087465 -2.991e-17 [382,] 9.819983e+01 0.00084861 -8.697e-17 [383,] 1.012140e+02 0.00082334 -2.666e-17 [384,] 1.043206e+02 0.00079882 -9.38e-17 [385,] 1.075226e+02 0.00077503 -7.517e-17 [386,] 1.108229e+02 0.00075195 -1.672e-16 [387,] 1.142245e+02 0.00072956 -1.771e-17 [388,] 1.177305e+02 0.00070783 -2.052e-16 [389,] 1.213441e+02 0.00068675 -4.559e-17 [390,] 1.250686e+02 0.00066630 4.914e-17 [391,] 1.289074e+02 0.00064646 -6.877e-17 [392,] 1.328641e+02 0.00062721 -8.279e-17 [393,] 1.369422e+02 0.00060853 -8.847e-17 [394,] 1.411455e+02 0.00059041 -1.288e-16 [395,] 1.454778e+02 0.00057282 -1.225e-16 [396,] 1.499430e+02 0.00055577 -3.718e-17 [397,] 1.545454e+02 0.00053922 -9.812e-18 [398,] 1.592890e+02 0.00052316 -3.229e-17 [399,] 1.641782e+02 0.00050758 -1.263e-16 [400,] 1.692174e+02 0.00049246 5.768e-17 [401,] 1.744114e+02 0.00047780 9.598e-18 [402,] 1.797647e+02 0.00046357 -1.217e-16 [403,] 1.852824e+02 0.00044976 -1.26e-16 [404,] 1.909694e+02 0.00043637 -9.519e-17 [405,] 1.968310e+02 0.00042337 -1.528e-16 [406,] 2.028725e+02 0.00041077 3.91e-17 [407,] 2.090994e+02 0.00039853 1.221e-16 [408,] 2.155175e+02 0.00038667 -4.148e-17 [409,] 2.221326e+02 0.00037515 -9.274e-17 [410,] 2.289507e+02 0.00036398 8.656e-17 [411,] 2.359781e+02 0.00035314 -1.988e-17 [412,] 2.432211e+02 0.00034262 -6.736e-17 [413,] 2.506865e+02 0.00033242 8.559e-17 [414,] 2.583811e+02 0.00032252 -5.875e-17 [415,] 2.663118e+02 0.00031292 4.992e-17 [416,] 2.744859e+02 0.00030360 -6.795e-17 [417,] 2.829110e+02 0.00029456 -1.886e-16 [418,] 2.915946e+02 0.00028578 -4.244e-17 [419,] 3.005447e+02 0.00027727 -7.257e-17 [420,] 3.097696e+02 0.00026902 -1.212e-17 [421,] 3.192776e+02 0.00026101 4.323e-17 [422,] 3.290775e+02 0.00025323 1.222e-18 [423,] 3.391782e+02 0.00024569 1.143e-17 [424,] 3.495888e+02 0.00023838 -3.239e-17 [425,] 3.603191e+02 0.00023128 -2.036e-17 [426,] 3.713787e+02 0.00022439 -1.649e-17 [427,] 3.827777e+02 0.00021771 -1.087e-16 [428,] 3.945266e+02 0.00021122 -3.306e-17 [429,] 4.066362e+02 0.00020493 -6.344e-17 [430,] 4.191174e+02 0.00019883 -9.779e-17 [431,] 4.319817e+02 0.00019291 -1.4e-16 [432,] 4.452409e+02 0.00018716 -4.602e-17 [433,] 4.589071e+02 0.00018159 -2.595e-17 [434,] 4.729927e+02 0.00017618 -6.214e-17 [435,] 4.875107e+02 0.00017094 -2.749e-19 [436,] 5.024742e+02 0.00016585 -1.218e-16 [437,] 5.178971e+02 0.00016091 -5.937e-17 [438,] 5.337934e+02 0.00015612 -1.306e-16 [439,] 5.501776e+02 0.00015147 -1.013e-16 [440,] 5.670646e+02 0.00014696 -1.361e-16 [441,] 5.844700e+02 0.00014258 2.628e-18 [442,] 6.024097e+02 0.00013833 -1.708e-16 [443,] 6.208999e+02 0.00013421 -4.021e-17 [444,] 6.399577e+02 0.00013022 -1.776e-16 [445,] 6.596005e+02 0.00012634 -9.53e-18 [446,] 6.798462e+02 0.00012258 4.098e-17 [447,] 7.007133e+02 0.00011893 -7.074e-17 [448,] 7.222209e+02 0.00011538 -7.154e-17 [449,] 7.443886e+02 0.00011195 5.213e-17 [450,] 7.672368e+02 0.00010861 1.19e-17 [451,] 7.907862e+02 0.00010538 -6.028e-17 [452,] 8.150585e+02 0.00010224 -1.509e-16 [453,] 8.400758e+02 0.00009920 -7.66e-17 [454,] 8.658610e+02 0.00009624 -6.065e-17 [455,] 8.924376e+02 0.00009338 -8.568e-17 [456,] 9.198299e+02 0.00009060 -3.954e-17 [457,] 9.480631e+02 0.00008790 -6.013e-17 [458,] 9.771628e+02 0.00008528 2.89e-17 [459,] 1.007156e+03 0.00008274 -5.789e-17 [460,] 1.038069e+03 0.00008028 -6.901e-17 [461,] 1.069932e+03 0.00007789 -5.364e-17 [462,] 1.102772e+03 0.00007557 -5.602e-17 [463,] 1.136620e+03 0.00007332 -2.722e-17 [464,] 1.171507e+03 0.00007113 5.041e-17 [465,] 1.207465e+03 0.00006902 -9.22e-17 [466,] 1.244527e+03 0.00006696 1.568e-18 [467,] 1.282727e+03 0.00006497 -1.345e-16 [468,] 1.322098e+03 0.00006303 4.978e-17 [469,] 1.362679e+03 0.00006115 -1.015e-16 [470,] 1.404505e+03 0.00005933 -2.064e-17 [471,] 1.447614e+03 0.00005757 -1.672e-17 [472,] 1.492047e+03 0.00005585 -1.072e-17 [473,] 1.537844e+03 0.00005419 4.032e-17 [474,] 1.585046e+03 0.00005257 4.8e-17 [475,] 1.633697e+03 0.00005101 -1.26e-17 [476,] 1.683842e+03 0.00004949 -6.588e-17 [477,] 1.735525e+03 0.00004802 -6.692e-17 [478,] 1.788795e+03 0.00004659 -1.216e-16 [479,] 1.843700e+03 0.00004520 -1.079e-16 [480,] 1.900291e+03 0.00004385 -7.245e-17 [481,] 1.958618e+03 0.00004255 6.521e-18 [482,] 2.018735e+03 0.00004128 -7.591e-17 [483,] 2.080698e+03 0.00004005 -1.044e-16 [484,] 2.144563e+03 0.00003886 -7.502e-17 [485,] 2.210388e+03 0.00003770 -1.287e-16 [486,] 2.278233e+03 0.00003658 1.402e-17 [487,] 2.348161e+03 0.00003549 -1.047e-16 [488,] 2.420235e+03 0.00003443 1.079e-16 [489,] 2.494521e+03 0.00003341 -2.871e-17 [490,] 2.571088e+03 0.00003241 1.029e-17 [491,] 2.650004e+03 0.00003145 -1.527e-16 [492,] 2.731343e+03 0.00003051 -6.406e-17 [493,] 2.815179e+03 0.00002960 -1.602e-18 [494,] 2.901587e+03 0.00002872 -1.671e-16 [495,] 2.990648e+03 0.00002786 -9e-17 [496,] 3.082443e+03 0.00002703 -2.641e-17 [497,] 3.177055e+03 0.00002623 -1.684e-17 [498,] 3.274571e+03 0.00002545 -1.052e-16 [499,] 3.375080e+03 0.00002469 -1.066e-16 [500,] 3.478674e+03 0.00002396 -1.862e-16 [501,] 3.585448e+03 0.00002324 -2.461e-16 [502,] 3.695499e+03 0.00002255 -1.546e-16 [503,] 3.808929e+03 0.00002188 -1.135e-17 [504,] 3.925839e+03 0.00002123 -2.003e-16 [505,] 4.046338e+03 0.00002059 -2.845e-17 [506,] 4.170536e+03 0.00001998 -9.463e-17 [507,] 4.298546e+03 0.00001939 1.223e-17 [508,] 4.430485e+03 0.00001881 -1.434e-16 [509,] 4.566474e+03 0.00001825 4.747e-17 [510,] 4.706636e+03 0.00001771 -7.073e-17 [511,] 4.851101e+03 0.00001718 -7.228e-17 [512,] 5.000000e+03 0.00001667 1.57e-17 > > p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", cutoffs = cuts) > ## and zoom in: > p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", cutoffs = cuts, ylim = yl2) > p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", cutoffs = cuts, ylim = yl2/20) > > if(do.pdf) { dev.off(); pdf("stirlerr-relErr_6-fin-2.pdf") } > > ## zoom in ==> {good for n >= 10} > p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", ylim = 2e-15*c(-1,1), + cutoffs = cuts)## old default cutoffs = c(15,35, 80, 500) > > if(do.pdf) { dev.off(); pdf("stirlerr-relErr_6-fin-3.pdf") } > showProc.time() Time (user system elapsed): 0.22 0.03 0.25 > > > ##-- April 20: have more terms up to S10 in stirlerr() --> can use more cutoffs > n <- n5m <- lseq(1/64, 5000, length = if(doExtras) 4096 else 512) > nM <- mpfr(n, if(doExtras) 2048L # a *lot* accuracy for stirlerr(nM,*) + else 512L) > ct10.1 <- c( 5.4, 7.5, 8.5, 10.625, 12.125, 20, 26, 60, 200, 3300)# till 2024-01-19 > ct10.2 <- c( 5.4, 7.9, 8.75,10.5 , 13, 20, 26, 60, 200, 3300) > cuts <- + ct12.1 <- c(5.22, 6.5, 7.0, 7.9, 8.75,10.5 , 13, 20, 26, 60, 200, 3300) > ## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ > ## 5.25 is "too small" but the direct formula is already really bad there, ... > st.nM <- roundMpfr(stirlerr(nM, use.halves=FALSE, ## << on purpose; + verbose=TRUE), precBits = 128) stirlerr(n): As 'n' is "mpfr", using "mpfr" & stirlerrM(): > ## NB: for x=xM ; `cutoffs` are *not* used. > p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", scheme = "R4.4_0") > p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", scheme = "R4.4_0", ylim = c(-1,1)*3e-15) > p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", scheme = "R4.4_0", ylim = c(-1,1)*1e-15) > > p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", scheme = "R4.4_0", abs=TRUE) > axis(1,at= 2:6, col=NA, col.axis=(cola <- "lightblue"), line=-3/4) > abline(v = 2:6, lty=3, col=cola) > if(FALSE)## using exact values sferr_halves[] *instead* of MPFR ones: ==> confirmation they lay on top + lines((0:30)/2, abs(stirlerr((0:30)/2, cutoffs=cuts, verbose=TRUE)/DPQ:::sferr_halves - 1), type="o", col=2,lwd=2) > > if(FALSE) ## nice (but unneeded) printout : + print(cbind(n = format(n, drop0trailing = TRUE), + stirlerr= format(st.,scientific=FALSE, digits=4), + relErr = signif(relE, 4)) + , quote=FALSE) > > showProc.time() Time (user system elapsed): 0.25 0.02 0.27 > > > ## ========== where should the cutoffs be ? =================================================== > > .stirl.cutoffs <- function(scheme) + eval(do.call(substitute, list(formals(stirlerr)$cutoffs, list(scheme = scheme)))) > drawCuts <- function(scheme, axis=NA, lty = 3, col = "skyblue", ...) { + abline(v = (ct <- .stirl.cutoffs(scheme)), lty=lty, col=col, ...) + if(is.finite(axis)) axisCuts(side = axis, at = ct, col=col, ...) + } > axisCuts <- function(scheme, side = 3, at = .stirl.cutoffs(scheme), col = "skyblue", line = -3/4, ...) + axis(side, at=at, labels=formatC(at), col.axis = col, col=NA, col.ticks=NA, line=line, ...) > mtextCuts <- function(cutoffs, scheme, ...) { + if(!missing(scheme)) cutoffs <- .stirl.cutoffs(scheme) + mtext(paste("cutoffs =", deparse1(cutoffs)), ...) + } > > > if(do.pdf) { dev.off(); pdf("stirlerr-tst_order_k.pdf") } > > mK <- 20L # := max(k) > ## order = k = 1:mK terms in series approx: > k <- 1:mK > n <- 2^seq(1, 28, by=1/16) > nM <- mpfr(n, 1024) > stnM <- stirlerr(nM) # the "true" values > > stirlOrd <- sapply(k, function(k.) stirlerr(n, order = k.)) > stirlO_lgcor <- cbind(stirlOrd, sapply(5:6, function(nal) lgammacor(n, nalgm = nal))) > relE <- asNumeric(stirlO_lgcor/stnM -1) # "true" relative error > > ## use a "smooth" but well visible polette : > palROBG <- colorRampPalette(c("red", "darkorange2", "blue", "seagreen"), space = "Lab") > palette(adjustcolor(palROBG(mK+2), 3/4)) > ## -- 2 lgamcor()'s > > (tit.k <- substitute(list( stirlerr(n, order=k) ~~"error", k == 1:mK), list(mK = mK))) list(stirlerr(n, order = k) ~ ~"error", k == 1:20L) > (tit.kA <- substitute(list(abs(stirlerr(n, order=k) ~~"error"), k == 1:mK), list(mK = mK))) list(abs(stirlerr(n, order = k) ~ ~"error"), k == 1:20L) > lgammacorTit <- function(...) mtext("+ lgammacor(x, 5) [bad] + lgammacor(x, 6) [good]", col=1, ...) > > matplotB(n, relE, cex=2/3, ylim = c(-1,1)*1e-13, col=k, + log = "x", xaxt="n", main = tit.k) > lgammacorTit() > eaxis(1, nintLog = 20) > drawCuts("R4.4_0") > > ## zoom in (ylim) > matplotB(n, relE, cex=2/3, ylim = c(-1,1)*5e-15, col=k, + log = "x", xaxt="n", main = tit.k) > lgammacorTit() > eaxis(1, nintLog = 20); abline(h = (-2:2)*2^-53, lty=3, lwd=1/2) > drawCuts("R4.4_0", axis = 3) > > ## log-log |rel.Err| -- "linear" > matplotB(n, abs19(relE), cex=2/3, col=k, ylim = c(8e-17, 1e-7), log = "xy", main=tit.kA) > mtext(paste("k =", deparse(k))) ; abline(h = 2^-(53:51), lty=3, lwd=1/2) > lgammacorTit(line=-1) > drawCuts("R4.4_0", axis = 3) > > ## zoom in -- still "large" {no longer carry the lgammacor() .. }: > n2c <- 2^seq(2, 8, by = 1/(if(doExtras) 256 else 32)) > nMc <- mpfr(n2c, 1024) > stnMc <- stirlerr(nMc) # the "true" values > stirlOrc <- sapply(k, function(k.) stirlerr(n2c, order = k.)) > relEc <- asNumeric(stirlOrc/stnMc -1) # "true" relative error > > matplotB(n2c, relEc, cex=2/3, ylim = c(-1,1)*1e-13, col=k, + log = "x", xaxt="n", main = tit.k) > eaxis(1, sub10 = 2) > drawCuts("R4.4_0", axis=3) > > ## log-log |rel.Err| -- "linear" > matplotB(n2c, abs19(relEc), cex=2/3, col=k, ylim = c(8e-17, 1e-3), log = "xy", main=tit.kA) > mtext(paste("k =", deparse(k))) ; abline(h = 2^-(53:51), lty=3, lwd=1/2) > drawCuts("R4.4_0", axis = 3) > > > ## zoom into the critical n region > nc <- seq(3.5, 11, by = 1/(if(doExtras) 128 else 16)) > ncM <- mpfr(nc, 256) > stncM <- stirlerr(ncM) # the "true" values > stirlO.c <- sapply(k, function(k) stirlerr(nc, order = k)) > relEc <- asNumeric(stirlO.c/stncM -1) # "true" relative error > > > ## log-log |rel.Err| -- "linear" > matplotB(nc, abs19(relEc), cex=2/3, col=k, ylim = c(2e-17, 1e-8), + log = "xy", xlab = quote(n), main = quote(abs(relErr(stirlerr(n, order==k))))) > mtext(paste("k =", deparse(k))) ; drawEps.h(lwd = 1/2) > lines(nc, abs19(asNumeric(stirlerr_simpl(nc, "R3" )/stncM - 1)), lwd=1.5, col=adjustcolor("thistle", .6)) > lines(nc, abs19(asNumeric(stirlerr_simpl(nc, "MM2")/stncM - 1)), lwd=4, col=adjustcolor(20, .4)) > ## lines(nc, abs19(asNumeric(stirlerr_simpl(nc,"MM2")/stncM - 1)), lwd=3, col=adjustcolor("purple", 2/3)) > legend(10^par("usr")[1], 1e-9, legend=paste0("k=", k), bty="n", lwd=2, + col=k, lty=1:5, pch= c(1L:9L, 0L, letters)[seq_along(k)]) > drawCuts("R4.4_0", axis=3) > > ## Zoom-in [only] > matplotB(nc, abs19(relEc), cex=2/3, col=k, ylim = c(4e-17, 1e-11), xlim = c(4.8, 6.5), + log = "xy", xlab = quote(n), main = quote(abs(relErr(stirlerr(n, order==k))))) > mtext(paste("k =", deparse(k))) ; drawEps.h(lwd = 1/2) > lines(nc, abs19(asNumeric(stirlerr_simpl(nc, "R3" )/stncM - 1)), lwd=1.5, col=adjustcolor("thistle", .6)) > lines(nc, abs19(asNumeric(stirlerr_simpl(nc, "MM2")/stncM - 1)), lwd=4, col=adjustcolor(20, .4)) > > k. <- k[-(1:6)] > legend("bottomleft", legend=paste0("k=", k.), bty="n", lwd=2, + col=k., lty=1:5, pch= c(1L:9L, 0L, letters)[k.]) > drawCuts("R4.4_0", axis=3) > > showProc.time() Time (user system elapsed): 0.62 0.03 0.65 > > ##--- Accuracy of "R4.4_0" ------------------------------------------------------- > > for(nc in list(seq(4.75, 28, by=1/(if(doExtras) 512 else 64)), # for a bigger pix + seq(4.75, 9, by=1/(if(doExtras) 1024 else 128)))) + { + ncM <- mpfr(nc, 1024) + stncM <- stirlerr(ncM) # the "true" values + stirl.440 <- stirlerr(nc, scheme = "R4.4_0") + stirl.3 <- stirlerr(nc, scheme = "R3") + relE440 <- asNumeric(relErrV(stncM, stirl.440)) + relE3 <- asNumeric(relErrV(stncM, stirl.3 )) + ## + plot(nc, abs19(relE440), xlab=quote(n), main = quote(abs(relErr(stirlerr(n, '"R4.4_0"')))), + type = "l", log = "xy", ylim = c(4e-17, 1e-13)) + mtextCuts(scheme="R4.4_0", cex=4/5) + drawCuts("R4.4_0", lty=2, lwd=2, axis=4) + drawEps.h() + if(max(nc) <= 10) abline(v = 5+(0:20)/10, lty=3, col=adjustcolor(4, 1/2)) + if(TRUE) { # but just so ... + c3 <- adjustcolor("royalblue", 1/2) + lines(nc, pmax(abs(relE3), 1e-18), col=c3) + title(quote(abs(relErr(stirlerr(n, '"R3"')))), adj=1, col.main = c3) + drawCuts("R3", lty=4, col=c3); mtextCuts(scheme="R3", adj=1, col=c3) + } + addOrd <- TRUE + addOrd <- dev.interactive(orNone=TRUE) + if(addOrd) { + if(max(nc) >= 100) { + i <- (15 <= nc & nc <= 85) # (no-op in [4.75, 7] !) + ni <- nc[i] + } else { i <- TRUE; ni <- nc } + for(k in 8:17) lines(ni, abs19(asNumeric(relErrV(stncM[i], stirlerr(ni, order=k)))), col=adjustcolor(k, 1/3)) + title(sub = "stirlerr(*, order k = 8:17)") + } + } ## for(nc ..) > > if(FALSE) + lines(nc, abs19(relE440))# *re* draw! > > showProc.time() Time (user system elapsed): 1.05 0 1.05 > > > palette("Tableau") > > ## Focus more: > stirlerrPlot <- function(nc, k, res=NULL, legend.xy = "left", full=TRUE, precB = 1024, + ylim = c(3e-17, 2e-13), cex = 5/4) { + stopifnot(require("Rmpfr"), require("graphics")) + if(is.list(res) && all(c("nc", "k", "relEc","splarelE") %in% names(res))) { ## do *not* recompute + list2env(res, envir = environment()) + } else { ## compute + stopifnot(is.finite(nc), nc > 0, length(nc) >= 100, k == as.integer(k), 0 <= k, k <= 20) + ncM <- mpfr(nc, precB) + stncM <- stirlerr(ncM) # the "true" values + stirlO.c <- sapply(k, function(k) stirlerr(nc, order = k)) + relEc <- asNumeric(stirlO.c/stncM -1) # "true" relative error + ## log |rel.Err| -- "linear" + ## smooth() on log-scale {and transform back}: + splarelE <- apply(log(abs19(relEc)), 2, function(y) exp(smooth.spline(y, df=4)$y)) + ## the direct formulas (default "R3", "MM2"): + arelEs0 <- abs(asNumeric(stirlerr_simpl(nc )/stncM - 1)) + arelEs2 <- abs(asNumeric(stirlerr_simpl(nc, "MM2")/stncM - 1)) + } + pch. <- c(1L:9L, 0L, letters)[k] + if(full) + matplotB(nc, abs19(relEc), col=k, pch = pch., cex=cex, ylim=ylim, log = "y", + xlab = quote(n), main = quote(abs(relErr(stirlerr(n, order==k))))) + else ## smooth only + matplotB(nc, splarelE, col=adjustcolor(k,2/3), pch=pch., lwd=2, cex=cex, ylim=ylim, log = "y", + xlab = quote(n), main = quote(abs(relErr(stirlerr(n, order==k))))) + mtext(paste("k =", deparse(k))) ; abline(h = 2^-(53:51), lty=3, lwd=1/2) + legend(legend.xy, legend=paste0("k=", k), bty="n", lwd=2, col=k, lty=1:5, pch = pch.) + abline(v = 5+(0:20)/10, lty=3, col=adjustcolor(10, 1/2)) + drawCuts("R4.4_0", axis=3) + if(full) { + matlines(nc, splarelE, col=adjustcolor(k,2/3), lwd=4) + lines(nc, pmax(arelEs0, 1e-19), lwd=1.5, col=adjustcolor( 2, 0.2)) + lines(nc, pmax(arelEs2, 1e-19), lwd=1, col=adjustcolor(10, 0.2)) + } + lines(nc, smooth.spline(arelEs0, df=12)$y, lwd=3, col= adjustcolor( 2, 1/2)) + lines(nc, smooth.spline(arelEs2, df=12)$y, lwd=3, col= adjustcolor(10, 1/2)) + invisible(list(nc=nc, k=k, relEc = relEc, splarelE = splarelE, arelEs0=arelEs0, arelEs2=arelEs2)) + } > > rr1 <- stirlerrPlot(nc = seq(4.75, 9.0, by=1/(if(doExtras) 1024 else 256)), + k = 7:20) > stirlerrPlot(res = rr1, full=FALSE, ylim = c(8e-17, 1e-13)) > if(interactive()) + stirlerrPlot(res = rr1) > > rr <- stirlerrPlot(nc = seq(5, 6.25, by=1/(if(doExtras) 2048 else 256)), k = 9:18) > stirlerrPlot(res = rr, full=FALSE, ylim = c(8e-17, 1e-13)) > > showProc.time() Time (user system elapsed): 1.48 0 1.48 > > > palette("default") > > if(do.pdf) { dev.off(); pdf("stirlerr-tst_order_k-vs-k1.pdf") } > > ##' Find 'cuts', i.e., a region c(k) +/- s(k) i.e. intervals [c(k) - s(k), c(k) + s(k)] > ##' where c(k) is such that relE(n=c(k), k) ~= eps) > > ##' 1. Find the c1(k) such that |relE(n, k)| ~= c1(k) * n^{-2k} > findC1 <- function(n, ks, e1 = 1e-15, e2 = 1e-5, res=NULL, + precBits = if(doExtras) 1024 else 256, do.plot = TRUE, ...) + { + if(is.list(res) && all(c("n", "ks", "arelE") %in% names(res))) { ## do *not* recompute, take from 'res': + list2env(res, envir = environment()) + } else { ## compute + stopifnot(require("Rmpfr"), + is.numeric(ks), ks == (k. <- as.integer(ks)), length(ks <- k.) >= 1, + length(e1) == 1L, length(e2) == 1L, is.finite(c(e1,e2)), e1 >= 0, e2 >= e1, + 0 <= ks, ks <= 20, is.numeric(n), n > 0, is.finite(n), length(n) >= 100) + nM <- mpfr(n, precBits) + stirM <- stirlerr(nM) # the "true" values + stirlOrd <- sapply(ks, function(k) stirlerr(n, order = k)) + arelE <- abs(asNumeric(stirlOrd/stirM -1)) # "true" relative error + } + arelE19 <- pmax(arelE, 1e-19) + ## log |rel.Err| -- "linear" + ## on log-scale {and transform back; for linear fit, only use values inside [e1, e2] + if(do.plot) # experi + matplotB(n, arelE19, log="xy", ...) + ## matplot(n, arelE19, type="l", log="xy", xlim = c(min(n), 20)) + + ## re-compute these, as they *also* depend on (e1, e2) + c1 <- vapply(seq_along(ks), function(i) { + k <- ks[i] + y <- arelE19[,i] + iUse <- e1 <= y & y <= e2 + if(sum(iUse) < 10) stop("only", sum(iUse), "values in [e1,e2]") + ## .lm.fit(cbind(1, log(n[iUse])), log(y[iUse]))$coefficients + ## rather, we *know* the error is c* n^{-2k} , i.e., + ## log |relE| = log(c) - 2k * log(n) + ## <==> c = exp( log|relE| + 2k * log(n)) + exp(mean(log(y[iUse]) + 2*k * log(n[iUse]))) + }, numeric(1)) + if(do.plot) { + drawEps.h() + for(i in seq_along(ks)) + lines(n, c1[i] * n^(-2*ks[i]), col=adjustcolor(i, 1/3), lwd = 4, lty = 2) + } + invisible(list(n=n, ks=ks, arelE = arelE, c1 = c1)) + } ## findC1() > > c1.Res <- findC1(n = 2^seq(2, 26, by = 1/(if(doExtras) 128 else 32)), ks = 1:18) > (s.c1.fil <- paste0("stirlerr-c1Res-", myPlatform(), ".rds")) [1] "stirlerr-c1Res-R-_88391__windows_WS2022x64_20_.rds" > saveRDS(c1.Res, file = s.c1.fil) > > > if(!exists("c1.Res")) { + c1.Res <- readRDS(s.c1.fil) + ## re-do the "default" plot of findC1(): + findC1(res = c1.Res, xaxt="n"); eaxis(1, sub10=2) + } > > ## the same, zoomed in: > findC1(res = c1.Res, xlim = c(4, 40), ylim = c(2e-17, 1e-12)) > ks <- c1.Res$ks; pch. <- c(1L:9L, 0L, letters)[ks] > legend("left", legend=paste0("k=", ks), bty="n", lwd=2, col=ks, lty=1:5, pch = pch.) > > ## smaller set : larger e1 : > c1.r2 <- findC1(res = c1.Res, xlim = c(4, 30), ylim = c(4e-17, 1e-13), e1 = 4e-15) > legend("left", legend=paste0("k=", ks), bty="n", lwd=2, col=ks, lty=1:5, pch = pch.) > > print(digits = 4, + cbind(ks, c1. = c1.Res$c1, c1.2 = c1.r2$c1, + relD = round(relErrV(c1.r2$c1, c1.Res$c1), 4))) ks c1. c1.2 relD [1,] 1 3.326e-02 3.331e-02 -0.0017 [2,] 2 9.528e-03 9.511e-03 0.0018 [3,] 3 7.037e-03 7.045e-03 -0.0012 [4,] 4 9.819e-03 9.804e-03 0.0016 [5,] 5 2.170e-02 2.167e-02 0.0014 [6,] 6 7.018e-02 6.957e-02 0.0087 [7,] 7 3.063e-01 3.037e-01 0.0084 [8,] 8 1.766e+00 1.739e+00 0.0156 [9,] 9 1.272e+01 1.255e+01 0.0134 [10,] 10 1.161e+02 1.126e+02 0.0310 [11,] 11 1.241e+03 1.218e+03 0.0187 [12,] 12 1.627e+04 1.584e+04 0.0272 [13,] 13 2.457e+05 2.427e+05 0.0123 [14,] 14 4.425e+06 4.334e+06 0.0210 [15,] 15 8.982e+07 8.822e+07 0.0181 [16,] 16 2.131e+09 2.080e+09 0.0245 [17,] 17 5.563e+10 5.467e+10 0.0175 [18,] 18 1.683e+12 1.639e+12 0.0270 > > c1.. <- c1.r2[c("ks", "c1")] # just the smallest part is needed here: > (s.c1.fil <- paste0("stirlerr-c1r2-", myPlatform(), ".rds")) [1] "stirlerr-c1r2-R-_88391__windows_WS2022x64_20_.rds" > saveRDS(c1.., file = s.c1.fil) # was "stirlerr-c1.rds" > > ## 2. Now, find the n(k) +/- se.n(k) intervals > ## Use these c1 from above > > ##' Given c1-results relErr(n,k); |relE(n,k) ~= c1 * n^{-2k} , find n such that > ##' |relE(n,k)| line {in log-log scale} cuts y = eps, i.e., n* such that |relE(n,k)| <= eps for all n >= n* > n.ep <- function(eps, c1Res, ks = c1Res$ks, c1 = c1Res$c1, ...) { + stopifnot(is.finite(eps), length(eps) == 1L, eps > 0, + length(ks) == length(c1), is.numeric(c1), is.integer(ks), ks >= 1) + ## n: given k, the location where the |relE(n,k)| line {in log-log} cuts y = eps + ## |relE(n,k) ~= c1 * n^{-2k} <==> + ## log|relE(n,k)| ~= log(c1) - 2k* log(n) <==> + ## c := mean{ exp( log|relE(n,k)| + 2k* log(n) ) } ------- see findC1() + ## now, solve for n : + ## c1 * n^{-2k} == eps + ## log(c1) - 2k* log(n) == log(eps) + ## log(n) == (log(eps) - log(c1)) / (-2k) <==> + ## n == exp((log(c1) - log(eps)) / 2k) + exp((log(c1) - log(eps))/(2*ks)) + } > > ## get c1.. > if(!exists("c1..")) c1.. <- readRDS("stirlerr-c1.rds") > > ne2 <- n.ep(2^-51, c1Res = c1..) ## ok > ne1 <- n.ep(2^-52, c1Res = c1..) > ne. <- n.ep(2^-53, c1Res = c1..) > > form <- function(n) format(signif(n, 3), scientific=FALSE) > data.frame(k = ks, ne2 = form(ne2), ne1 = form(ne1), ne. = form(ne.), + cutoffs = form(rev(.stirl.cutoffs("R4.4_0")[-1]))) k ne2 ne1 ne. cutoffs 1 1 8660000.00 12200000.00 17300000.00 17400000.00 2 2 2150.00 2560.00 3040.00 3700.00 3 3 159.00 178.00 200.00 200.00 4 4 46.60 50.80 55.40 81.00 5 5 23.40 25.10 26.90 36.00 6 6 15.20 16.10 17.10 25.00 7 7 11.50 12.10 12.70 19.00 8 8 9.43 9.85 10.30 14.00 9 9 8.20 8.52 8.86 11.00 10 10 7.42 7.68 7.95 9.50 11 11 6.89 7.11 7.34 8.80 12 12 6.53 6.72 6.91 8.25 13 13 6.27 6.44 6.62 7.60 14 14 6.10 6.25 6.41 7.10 15 15 5.98 6.12 6.26 6.50 16 16 5.90 6.03 6.16 6.50 17 17 5.85 5.97 6.10 6.50 18 18 5.83 5.94 6.06 6.50 > > ## ------- Linux F 36/38 x86_64 (nb-mm5|v-lynne) > ## k ne2 ne1 ne. cutoffs > ## 1 8660000.00 12200000.00 17300000.00 17400000.00 > ## 2 2150.00 2560.00 3040.00 3700.00 > ## 3 159.00 178.00 200.00 200.00 > ## 4 46.60 50.80 55.40 81.00 > ## 5 23.40 25.10 26.90 36.00 > ## 6 15.20 16.10 17.10 25.00 > ## 7 11.50 12.10 12.70 19.00 > ## 8 9.43 9.85 10.30 14.00 > ## 9 8.20 8.53 8.86 11.00 > ## 10 7.42 7.68 7.95 9.50 > ## 11 6.89 7.11 7.34 8.80 > ## 12 6.53 6.72 6.92 8.25 > ## 13 6.27 6.44 6.62 7.60 > ## 14 6.10 6.25 6.41 7.10 > ## 15 5.98 6.12 6.26 6.50 * (not used) > ## 16 5.90 6.03 6.16 6.50 * " " > ## 17 5.85 5.97 6.10 6.50 * " " > ## 18 5.83 5.95 6.06 6.50 << used all the way down to 5.25 > > ## ok --- correct order of magnitude ! --- good! > > > ## 2b. find *interval* around the 'n(eps)' values > > ## -- Try simply > d.k <- ne. - ne1 > ## interval > int.k <- cbind(ne1 - d.k, + ne1 + d.k) > ## look at e.g. > data.frame(k=ks, `n(k)` = form(ne1), int = form(int.k)) k n.k. int.1 int.2 1 1 12200000.00 7180000.00 17300000.00 2 2 2560.00 2070.00 3040.00 3 3 178.00 156.00 200.00 4 4 50.80 46.20 55.40 5 5 25.10 23.30 26.90 6 6 16.10 15.20 17.10 7 7 12.10 11.40 12.70 8 8 9.85 9.41 10.30 9 9 8.52 8.19 8.86 10 10 7.68 7.41 7.95 11 11 7.11 6.88 7.34 12 12 6.72 6.52 6.91 13 13 6.44 6.27 6.62 14 14 6.25 6.10 6.41 15 15 6.12 5.98 6.26 16 16 6.03 5.90 6.16 17 17 5.97 5.85 6.10 18 18 5.94 5.83 6.06 > ## k n.k. int.1 int.2 > ## 1 12200000.00 7180000.00 17300000.00 > ## 2 2560.00 2070.00 3040.00 > ## 3 178.00 156.00 200.00 > ## 4 50.80 46.20 55.40 > ## 5 25.10 23.30 26.90 > ## 6 16.10 15.20 17.10 > ## 7 12.10 11.40 12.70 > ## 8 9.85 9.41 10.30 > ## 9 8.53 8.19 8.86 > ## 10 7.68 7.41 7.95 > ## 11 7.11 6.88 7.34 > ## 12 6.72 6.52 6.92 > ## 13 6.44 6.27 6.62 > ## 14 6.25 6.10 6.41 > ## 15 6.12 5.98 6.26 > ## 16 6.03 5.90 6.16 > ## 17 5.97 5.85 6.10 > ## 18 5.95 5.83 6.06 > > ##' as function {well, *not* computing c1.k from scratch > nInt <- function(k, c1.k, ep12 = 2^-(52:53)) { + if(length(k) == 1L) { # special convention to call for *one* k, with c1.k vector + stopifnot(k == (k <- as.integer(k)), k >= 1, length(c1.k) >= k) + c1.k <- c1.k[k] + } + + ## see n.ep() above + n_ <- function(eps, k, c1) { + stopifnot(is.finite(eps), length(eps) == 1L, eps > 0, + length(k) == length(c1), is.numeric(c1), is.integer(k), k >= 1) + exp((log(c1) - log(eps))/(2*k)) + } + + ne1 <- n_(ep12[1], k, c1.k) + ne. <- n_(ep12[2], k, c1.k) + d.k <- ne. - ne1 + stopifnot(d.k > 0) + ## interval: {"fudge" 0.5 / 2.5} from results -- also 'noLdbl' gives *quite* different pic!: + odd <- k %% 2 == 1 + cbind(ne1 - ifelse( odd & noLdbl, 1.5, 0.5) * d.k, + ne1 + ifelse(!odd , 8, 2.5) * d.k) + } > > nInt(k= 1, c1..$c1) [,1] [,2] [1,] 9711841 24932464 > nInt(k= 2, c1..$c1) [,1] [,2] [1,] 2316.234 6430.572 > nInt(k=18, c1..$c1) [,1] [,2] [1,] 5.886718 6.869019 > > nints.k <- nInt(ks, c1..$c1) > ## for printing > form(as.data.frame( nints.k )) V1 V2 1 9710000.00 24900000.00 2 2320.00 6430.00 3 167.00 232.00 4 48.50 87.50 5 24.20 29.60 6 15.70 23.80 7 11.70 13.60 8 9.63 13.30 9 8.36 9.36 10 7.54 9.84 11 6.99 7.68 12 6.62 8.29 13 6.36 6.88 14 6.17 7.51 15 6.05 6.48 16 5.96 7.09 17 5.91 6.28 18 5.89 6.87 > > ## (-.5 , +2.5) ## originally ( -1, +1) > ## 1 9710000.00 24900000.00 # 7180000.00 17300000.00 > ## 2 2320.00 3770.00 # 2070.00 3040.00 > ## 3 167.00 232.00 # 156.00 200.00 > ## 4 48.50 62.30 # 46.20 55.40 > ## 5 24.20 29.60 # 23.30 26.90 > ## 6 15.70 18.50 # 15.20 17.10 > ## 7 11.70 13.60 # 11.40 12.70 > ## 8 9.63 10.90 # 9.41 10.30 > ## 9 8.36 9.36 # 8.19 8.86 > ## 10 7.54 8.36 # 7.41 7.95 > ## 11 7.00 7.68 # 6.88 7.34 > ## 12 6.62 7.21 # 6.52 6.92 > ## 13 6.36 6.88 # 6.27 6.62 > ## 14 6.17 6.64 # 6.10 6.41 > ## 15 6.05 6.48 # 5.98 6.26 > ## 16 5.96 6.36 # 5.90 6.16 > ## 17 5.91 6.28 # 5.85 6.10 > ## 18 5.89 6.24 # 5.83 6.06 > > ## 3. Then for each of the intervals, compare order k vs k+1 > ## -- ==> optimal cutoff { how much platform dependency ?? } > > ### Here, compute only > find1cuts <- function(k, c1, + n = nInt(k, c1), # the *set* of n's or the 'range' + len.n = if(doExtras) 1000 else 250, + precBits = if(doExtras) 1024 else 256, nM = mpfr(n, precBits), + stnM = stirlerr(nM), + stirlOrd = sapply(k+(0:1), function(.k.) stirlerr(n, order = .k.)), + relE = asNumeric(stirlOrd/stnM -1), # "true" relative error for the {k, k+1} + do.spl=TRUE, df.spline = 9, # df = 5 gives 2 cutpoints for k==1 + do.low=TRUE, f.lowess = 0.2, + do.cobs = require("cobs", quietly=TRUE), tau = 0.90) + { + ## check relErrV( stirlerr(n, order=k ) vs + ## stirlerr(n, order=k+1) + if(length(n) == 2L) + if(n[1] < n[2]) n <- seq(n[1], n[2], length.out = len.n) + else stop("'n' must be *increasing") + force(relE) + y <- abs19(relE) + ## NB: all smoothing --- as in stirlerrPlot() above -- should happen in log-space + ## (log(abs19(relEc)), 2, function(y) exp(smooth.spline(y, df=4)$y)) + ly <- log(y) # == log(abs19(relE)) == log(max(|r|, 1e-19)) + if(do.spl) {## + s1 <- exp(smooth.spline(ly[,1], df=df.spline)$y) + s2 <- exp(smooth.spline(ly[,2], df=df.spline)$y) + } + if(do.low) { ## lowess + s1l <- exp(lowess(ly[,1], f=f.lowess)$y) + s2l <- exp(lowess(ly[,2], f=f.lowess)$y) + } + ## also use cobs() splines for the 90% quantile !! + EE <- environment() # so can set do.cobs to FALSE in case of error + if(do.cobs) { ## <==> require("cobs") # yes, this is in tests/ + cobsF <- function(Y) cobs(n, Y, tau=tau, nknots = 6, lambda = -1, + print.warn=FALSE, print.mesg=FALSE) + ## sparseM::chol() now gives error when it gave warning {about singularity} + cobsF <- function(Y) { + r <- tryCatch(cobs(n, Y, tau=tau, nknots = 6, lambda = -1, + print.warn=FALSE, print.mesg=FALSE), + error = identity) + if(inherits(r, "error")) { + assign("do.cobs", FALSE, envir = EE) # and return + list(fitted = FALSE) + } + else + r + } + cs1 <- exp(cobsF(ly[,1])$fitted) + cs2 <- exp(cobsF(ly[,2])$fitted) + } + smooths <- list(spl = if(do.spl ) cbind(s1, s2 ), + low = if(do.low ) cbind(s1l,s2l), + cobs= if(do.cobs) cbind(cs1,cs2)) + ## diffL <- list(spl = if(do.spl ) s2 -s1, + ## low = if(do.low ) s2l-s1l, + ## cobs= if(do.cobs) cs2-cs1) + ### FIXME: simplification does not always *work* -- (i, n.) are not always ok + ### ------ notably within R-devel-no-ldouble {hence probably macOS M1 ... ..} + sapply(smooths, function(s12) { d <- s12[,2] - s12[,1] + ## typically a (almost or completely) montone increasing function, crossing zero *once* + ## compute cutpoint: + ## i := the first n[i] with d(n[i]) >= 0 + i <- which(d >= 0)[1] + if(length(i) == 1L && !is.na(i) && i > 1L) { + i_ <- i - 1L # ==> d(n[i_]) < 0 + ## cutpoint must be in [n[i_], n[i]] --- do linear interpolation + n. <- n[i_] - (n[i] - n[i_])* d[i_] / (d[i] - d[i_]) + } else { + if(length(i) != 1L) i <- -length(i) + n. <- NA_integer_ + } + c(i=i, n.=n.) + }) -> n.L + + list(k=k, n=n, relE = unname(relE), smooths=smooths, i.n = n.L) + } ## find1cuts() > > k. <- 1:15 > system.time( + ## Failed on lynne [2024-06-04, R 4.4.1 beta] with + ## Error in .local(x, ...) : insufficient space ---> SparseM :: chol(<..>) + ## ==> now we catch this inside find1cuts(): + resL <- lapply(setNames(,k.), function(k) find1cuts(k=k, c1=c1..$c1)) + ) ## -- warnings, notably from cobs() not converging user system elapsed 8.08 0.29 8.36 Warning messages: 1: In chol(, *): Replaced 5 tiny diagonal entries by 'Large' 2: In chol(, *): Replaced 7 tiny diagonal entries by 'Large' 3: In chol(, *): Replaced 7 tiny diagonal entries by 'Large' 4: In cobs(n, Y, tau = tau, nknots = 6, lambda = -1, print.warn = FALSE, : drqssbc2(): Not all flags are normal (== 1), ifl : 11111111124182224242424242424242424242424 5: In chol(, *): Replaced 5 tiny diagonal entries by 'Large' 6: In chol(, *): Replaced 7 tiny diagonal entries by 'Large' 7: In chol(, *): Replaced 7 tiny diagonal entries by 'Large' 8: In cobs(n, Y, tau = tau, nknots = 6, lambda = -1, print.warn = FALSE, : drqssbc2(): Not all flags are normal (== 1), ifl : 1111111111182424242424242424242424242424 > ## needs 12 sec (!!) user system elapsed = > (s.find15.fil <- paste0("stirlerr-find1_1-15_", myPlatform(), ".rds")) [1] "stirlerr-find1_1-15_R-_88391__windows_WS2022x64_20_.rds" > ## now we catch cobs() errors {from SparseM::chol}, > ## okCuts <- !inherits(resL, "error") > ## if(okCuts) { > saveRDS(resL, file = s.find15.fil) # was "stirlerr-find1_1-15.rds" > ## } else traceback() > ## 11: stop(mess) > ## 10: .local(x, ...) > ## 9: chol(e, tmpmax = tmpmax, nsubmax = nsubmax, nnzlmax = nnzlmax) > ## 8: chol(e, tmpmax = tmpmax, nsubmax = nsubmax, nnzlmax = nnzlmax) > ## 7: rq.fit.sfnc(Xeq, Yeq, Xieq, Yieq, tau = tau, rhs = rhs, control = rqCtrl) > ## 6: drqssbc2(x, y, w, pw = pw, knots = knots, degree = degree, Tlambda = if (select.lambda) lambdaSet else lambda, > ## constraint = constraint, ptConstr = ptConstr, maxiter = maxiter, > ## trace = trace - 1, nrq, nl1, neqc, niqc, nvar, tau = tau, > ## select.lambda = select.lambda, give.pseudo.x = keep.x.ps, > ## rq.tol = rq.tol, tol.0res = tol.0res, print.warn = print.warn) > ## 5: cobs(n, Y, tau = tau, nknots = 6, lambda = -1, print.warn = FALSE, > ## print.mesg = FALSE) at stirlerr-tst.R!udBzpT#32 > ## 4: cobsF(ly[, 1]) at stirlerr-tst.R!udBzpT#34 > ## 3: find1cuts(k = k, c1 = c1..$c1) at #1 > ## 2: FUN(X[[i]], ...) > ## 1: lapply(setNames(, k.), function(k) find1cuts(k = k, c1 = c1..$c1)) > > ok1cutsLst <- function(res) { + stopifnot(is.list(res), sapply(res, is.list)) # must be list of lists + cobsL <- lapply(lapply(res, `[[`, "smooths"), `[[`, "cobs") + vapply(cobsL, is.array, NA) + } > (resLok <- ok1cutsLst(resL)) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE > ## 1 2 3 4 5 6 ...... 15 > ## FALSE TRUE TRUE TRUE TRUE TRUE ...... TRUE > > > if(FALSE) { + ## e.g. in R-devel-no-ldouble: + (r1 <- find1cuts(k=1, c1=c1..$c1))$i.n # list --- no longer ok {SparseM::chol -> insufficient space} + (r2 <- find1cuts(k=2, c1=c1..$c1))$i.n # "good" + (r3 <- find1cuts(k=3, c1=c1..$c1))$i.n # 3-vector: only 3x 'i' == 1 + } > > ## if(okCuts) { > mult.fig(15, main = "stirlerr(n, order=k) vs order = k+1")$old.par -> opar > invisible(lapply(resL[resLok], plot1cuts)) > ## plus the "last" ones {also showing that k=15 is worse here anyway than k=17} > str(r17 <- find1cuts(k=17, n = seq(5.1, 6.5, length.out = if(doExtras) 1500 else 400), c1=c1..$c1)) List of 5 $ k : num 17 $ n : num [1:400] 5.1 5.1 5.11 5.11 5.11 ... $ relE : num [1:400, 1:2] 5.31e-14 5.18e-14 5.08e-14 4.95e-14 4.84e-14 ... $ smooths:List of 3 ..$ spl : num [1:400, 1:2] 5.29e-14 5.17e-14 5.05e-14 4.94e-14 4.83e-14 ... .. ..- attr(*, "dimnames")=List of 2 .. .. ..$ : NULL .. .. ..$ : chr [1:2] "s1" "s2" ..$ low : num [1:400, 1:2] 5.29e-14 5.17e-14 5.06e-14 4.95e-14 4.84e-14 ... .. ..- attr(*, "dimnames")=List of 2 .. .. ..$ : NULL .. .. ..$ : chr [1:2] "s1l" "s2l" ..$ cobs: num [1:400, 1:2] 5.32e-14 5.20e-14 5.08e-14 4.97e-14 4.86e-14 ... .. ..- attr(*, "dimnames")=List of 2 .. .. ..$ : NULL .. .. ..$ : chr [1:2] "cs1" "cs2" $ i.n : int [1:2, 1:3] 1 NA 1 NA 1 NA ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:2] "i" "n." .. ..$ : chr [1:3] "spl" "low" "cobs" > plot1cuts(r17) # no-ldouble is *very* different than normal: k *much better* than k+1 > > (s.find17.fil <- paste0("stirlerr-find1_17_", myPlatform(), ".rds")) [1] "stirlerr-find1_17_R-_88391__windows_WS2022x64_20_.rds" > saveRDS(r17, file = s.find17.fil) # was "stirlerr-find1_17.rds" > ## } > > ## ditto for tau = 0.8 {quantile for cobs()}: > system.time( + resL.8 <- lapply(setNames(,k.), function(k) find1cuts(k=k, c1=c1..$c1, tau = 0.80)) + ) # warnings from cobs {again 12.1 sec} user system elapsed 6.78 0.24 7.01 Warning messages: 1: In chol(, *): Replaced 5 tiny diagonal entries by 'Large' 2: In chol(, *): Replaced 7 tiny diagonal entries by 'Large' 3: In chol(, *): Replaced 7 tiny diagonal entries by 'Large' 4: In cobs(n, Y, tau = tau, nknots = 6, lambda = -1, print.warn = FALSE, : drqssbc2(): Not all flags are normal (== 1), ifl : 1111111111202324242424242424242424242424 5: In chol(, *): Replaced 5 tiny diagonal entries by 'Large' 6: In chol(, *): Replaced 7 tiny diagonal entries by 'Large' 7: In chol(, *): Replaced 7 tiny diagonal entries by 'Large' 8: In cobs(n, Y, tau = tau, nknots = 6, lambda = -1, print.warn = FALSE, : drqssbc2(): Not all flags are normal (== 1), ifl : 1111111919181192424242424242424242424242424 9: In cobs(n, Y, tau = tau, nknots = 6, lambda = -1, print.warn = FALSE, : drqssbc2(): Not all flags are normal (== 1), ifl : 111111111111111111111919191918 > (resL8ok <- ok1cutsLst(resL.8)) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE > > mult.fig(15, main = "stirlerr(n, order=k) vs order = k+1 -- tau = 0.8") > invisible(lapply(resL.8[resL8ok], plot1cuts)) > ## plus "last" one > str(r17.8 <- find1cuts(k=17, n = seq(5.1, 6.5, length.out = if(doExtras) 1500 else 400), c1=c1..$c1, tau = 0.80)) List of 5 $ k : num 17 $ n : num [1:400] 5.1 5.1 5.11 5.11 5.11 ... $ relE : num [1:400, 1:2] 5.31e-14 5.18e-14 5.08e-14 4.95e-14 4.84e-14 ... $ smooths:List of 3 ..$ spl : num [1:400, 1:2] 5.29e-14 5.17e-14 5.05e-14 4.94e-14 4.83e-14 ... .. ..- attr(*, "dimnames")=List of 2 .. .. ..$ : NULL .. .. ..$ : chr [1:2] "s1" "s2" ..$ low : num [1:400, 1:2] 5.29e-14 5.17e-14 5.06e-14 4.95e-14 4.84e-14 ... .. ..- attr(*, "dimnames")=List of 2 .. .. ..$ : NULL .. .. ..$ : chr [1:2] "s1l" "s2l" ..$ cobs: num [1:400, 1:2] 5.31e-14 5.20e-14 5.08e-14 4.97e-14 4.85e-14 ... .. ..- attr(*, "dimnames")=List of 2 .. .. ..$ : NULL .. .. ..$ : chr [1:2] "cs1" "cs2" $ i.n : int [1:2, 1:3] 1 NA 1 NA 1 NA ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:2] "i" "n." .. ..$ : chr [1:3] "spl" "low" "cobs" > plot1cuts(r17.8) > par(opar) > > n.mn <- sapply(resL, `[[`, "i.n", simplify = "array") > n.mn8 <- sapply(resL.8, `[[`, "i.n", simplify = "array") > ## of course, only the cobs part differs: > ## ... when I later see errors here ... what's going on?? > str(n.mn) # all NULL for "no-double" !! num [1:2, 1:3, 1:15] 91 15211420 113 16523626 98 ... - attr(*, "dimnames")=List of 3 ..$ : chr [1:2] "i" "n." ..$ : chr [1:3] "spl" "low" "cobs" ..$ : chr [1:15] "1" "2" "3" "4" ... > dim(n.mn["n.",,]) [1] 3 15 > str(n.mn8["n.", "cobs",]) Named num [1:15] 1.58e+07 NA 1.98e+02 NA 2.64e+01 ... - attr(*, "names")= chr [1:15] "1" "2" "3" "4" ... > > sessionInfo() R Under development (unstable) (2025-07-07 r88391 ucrt) Platform: x86_64-w64-mingw32/x64 Running under: Windows Server 2022 x64 (build 20348) Matrix products: default LAPACK version 3.12.1 locale: [1] LC_COLLATE=C LC_CTYPE=German_Germany.utf8 LC_MONETARY=C [4] LC_NUMERIC=C 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] cobs_1.3-9-1 sfsmisc_1.1-20 Rmpfr_1.1-0 gmp_0.7-5 DPQ_0.6-0 loaded via a namespace (and not attached): [1] MASS_7.3-65 DPQmpfr_0.3-3 compiler_4.6.0 Matrix_1.7-3 quantreg_6.1 [6] tools_4.6.0 SparseM_1.84-2 survival_3.8-3 MatrixModels_0.5-4 splines_4.6.0 [11] grid_4.6.0 lattice_0.22-7 > > form(data.frame(k = k., cutoffs = rev(.stirl.cutoffs("R4.4_0")[-1])[k.], + t(n.mn ["n.",,]), cobs.80 = n.mn8["n.", "cobs",])) k cutoffs spl low cobs cobs.80 1 1 17400000.00 15200000.00 16500000.00 15600000.00 15800000.00 2 2 3700.00 6200.00 NA NA NA 3 3 200.00 194.00 195.00 198.00 198.00 4 4 81.00 NA NA 68.60 NA 5 5 36.00 26.60 26.90 26.20 26.40 6 6 25.00 23.40 NA 21.30 22.10 7 7 19.00 12.70 13.10 12.80 12.90 8 8 14.00 12.90 NA 12.50 13.10 9 9 11.00 8.93 8.97 8.87 8.80 10 10 9.50 9.56 NA 9.31 9.33 11 11 8.80 7.28 7.32 7.30 7.26 12 12 8.25 NA NA NA NA 13 13 7.60 6.51 6.52 6.50 6.54 14 14 7.10 7.38 NA 7.34 7.50 15 15 6.50 6.09 6.09 6.09 6.08 > ## k cutoffs spl low cobs cobs.80 > ## 1 1 17400000.00 15600000.00 15800000.00 15800000.00 15700000.00 > ## 2 2 3700.00 NA NA NA NA == from tau=0.90 I'd use ~ 5000 > ## 3 3 200.00 200.00 201.00 206.00 205.00 205 > ## 4 4 81.00 NA NA 86.60 86.20 86 > ## 5 5 36.00 27.20 27.10 27.00 27.10 27 > ## 6 6 25.00 23.50 NA NA NA 23.5 > ## 7 7 19.00 12.80 12.80 12.80 12.80 12.8 > ## 8 8 14.00 NA NA 12.30 12.80 12.3 > ## 9 9 11.00 8.91 8.95 8.89 8.86 8.9 > ## 10 10 9.50 9.69 NA 9.72 NA --skip-- > ## 11 11 8.80 7.26 7.31 7.32 7.29 7.3 > ## 12 12 8.25 8.16 NA 8.10 7.76 --skip-- > ## 13 13 7.60 6.51 6.53 6.51 6.52 6.52 > ## 14 14 7.10 7.43 NA 7.47 7.43 --skip-- > ## 15 15 6.50 6.10 6.10 6.08 6.09 6.10 > ## 16 --skip-- > ## 17 5.25 or something all the way dow > > ##---- In the no-ldouble case --- this is very different: > ## k cutoffs spl low cobs cobs.80 > ## 1 17400000.00 8580000.00 8580000.00 8370000.00 8370000.00 > ## 2 3700.00 NA NA 5890.00 6180.00 > ## 3 200.00 159.00 159.00 160.00 159.00 > ## 4 81.00 87.10 87.10 84.80 86.00 > ## 5 36.00 23.50 23.50 23.50 23.50 > ## 6 25.00 23.70 23.70 NA NA > ## 7 19.00 11.50 11.50 11.50 11.50 > ## 8 14.00 NA NA 12.90 13.00 > ## 9 11.00 8.20 8.20 8.21 8.21 > ## 10 9.50 NA NA NA 9.69 > ## 11 8.80 6.84 6.84 6.85 6.83 > ## 12 8.25 NA NA 8.27 7.95 > ## 13 7.60 NA NA NA NA > ## 14 7.10 NA NA 7.43 7.40 > ## 15 6.50 NA NA NA NA > > > ## M1 macOS is similar to no-ldouble : "==" if equaln > ## M1 mac | aarch64-apple-darwin20 | macOS Ventura 13.3.1 | R-devel (2024-01-29 r85841) ; LAPACK version 3.12.0 > ## k cutoffs spl low cobs cobs.80 > ## 1 1 17400000.00 8580000.00 8580000.00 8370000.00 8370000.00 == > ## 2 2 3700.00 NA NA 5900.00 NA ~= > ## 3 3 200.00 159.00 159.00 160.00 159.00 == > ## 4 4 81.00 87.10 87.10 84.80 86.00 == > ## 5 5 36.00 23.50 23.50 23.50 23.50 > ## 6 6 25.00 23.70 23.70 NA NA == > ## 7 7 19.00 11.50 11.50 11.50 11.50 > ## 8 8 14.00 NA NA 12.90 13.00 == > ## 9 9 11.00 8.20 8.20 8.21 8.21 > ## 10 10 9.50 NA NA NA 9.69 == > ## 11 11 8.80 6.84 6.84 6.85 6.83 > ## 12 12 8.25 NA NA 8.27 7.95 == > ## 13 13 7.60 NA NA NA NA == > ## 14 14 7.10 NA NA 7.43 7.40 == > ## 15 15 6.50 NA NA NA NA == > > > > > ## ========== Try a slightly better direct formula ====================================================== > > ## after some trial error: gamma(n+1) = n*gamma(n) ==> lgamma(n+1) = lgamma(n) + log(n) > ## stirlerrD2 <- function(n) lgamma(n) + n*(1-(l.n <- log(n))) + (l.n - log(2*pi))/2 > > > palette("default") > > if(do.pdf) { dev.off(); pdf("stirlerr-tst_others.pdf") } > > n <- n5m # (1/64 ... 5000) -- goes with st.nM > p.stirlerrDev(n=n, stnM=st.nM, cex=1/4, type="o", cutoffs = cuts, abs=TRUE) > axis(1,at= 2:6, col=NA, col.axis=(cola <- "lightblue"), line=-3/4) > abline(v = 2:6, lty=3, col=cola) > i.n <- 1 <= n & n <= 15 ; cr2 <- adjustcolor(2, 0.6) > lines(n[i.n], abs(relErrV(st.nM[i.n], stirlerr_simpl(n[i.n], "MM2" ))), col=cr2, lwd=2) > legend(20, 1e-13, legend = quote( relErr(stirlerr_simpl(n, '"MM2"'))), col=cr2, lwd=3, bty="n") > > > n <- seq(1,6, by= 1/(if(doExtras) 200 else 64)) > stM <- stirlerr(mpfr(n, if(doExtras) 512 else 256), use.halves = FALSE) > relE <- asNumeric(relErrV(stM, cbind(stD = stirlerr_simpl(n), stD2 = stirlerr_simpl(n, "MM2")))) > signif(apply(abs(relE), 2, summary), 4) stD stD2 Min. 7.568e-18 6.131e-17 1st Qu. 1.955e-15 1.501e-15 Median 7.465e-15 4.893e-15 Mean 1.598e-14 9.633e-15 3rd Qu. 2.311e-14 1.346e-14 Max. 9.051e-14 7.273e-14 > ## stD stD2 > ## Min. 3.093e-17 7.081e-18 > ## 1st Qu. 2.777e-15 1.936e-15 > ## Median 8.735e-15 5.238e-15 > ## Mean 1.670e-14 1.017e-14 .. well "67% better" > ## 3rd Qu. 2.380e-14 1.408e-14 > ## Max. 1.284e-13 9.382e-14 > > c2 <- adjustcolor(1:2, 0.6) > matplotB(n, abs19(relE), type="o", cex=3/4, log="y", ylim = c(8e-17, 1.3e-13), yaxt="n", col=c2) > eaxis(2) > abline(h = 2^-53, lty=1, col="gray") > smrelE <- apply(abs(relE), 2, \(y) lowess(n, y, f = 0.1)$y) > matlines(n, smrelE, lwd=3, lty=1) > legend("topleft", legend = expression(stirlerr_simpl(n), stirlerr_simple(n, "MM2")), + bty='n', col=c2, lwd=3, lty=1) > drawEps.h(-(53:48)) > > > ## should we e.g., use interpolation spline through sfserr_halves[] for n <= 7.5 > ## -- doing the interpolation on the log(1 - 12*x*stirlerr(x)) vs log2(x) scale -- maybe ? > stirL <- curve(1-12*x*stirlerr(x, cutoffs = cuts, verbose=TRUE), 1/64, 8, log="xy", n=if(doExtras) 2048 else 512) stirlerr(n, cutoffs = 5.22,6.5,7,7.9,8.75,10.5,13,20,26,60,200,3300) : case I (n <= 5.22), using direct formula for n= num [1:477] 0.0156 0.0158 0.016 0.0162 0.0164 ... case II (n > 5.22 ), 12 cutoffs: ( 5.22, 6.5, 7, 7.9, 8.75, 10.5, 13, 20, 26, 60, 200, 3300 ): n in cutoff intervals: (5.22,6.5] (6.5,7] (7,7.9] (7.9,8.75] (8.75,10.5] (10.5,13] (13,20] 17 7 9 2 0 0 0 (20,26] (26,60] (60,200] (200,3.3e+03] (3.3e+03,Inf] 0 0 0 0 0 > ## just need "true" values for x = 2^-(6,5,4,3,2) in addition to those we already have at x = 1/2, 1.5, 2, 2.5, ..., 7.5, 8 > xM <- mpfr(stirL$x, if(doExtras) 2048 else 256) > stilEM <- roundMpfr(1 - 12*xM*stirlerr(xM, verbose=TRUE), 128) stirlerr(n): As 'n' is "mpfr", using "mpfr" & stirlerrM(): > relEsml <- relErrV(stilEM, stirL$y) > plot(stirL$x, relEsml) # again: stirlerr() is "limited.." for ~ [2, 6] > > ## The function to interpolate: > plot(asNumeric( (stilEM)) ~ x, data = stirL, type = "l", log="x") > plot(asNumeric(sqrt(stilEM)) ~ x, data = stirL, type = "l", log="x") > plot(asNumeric( log(stilEM)) ~ x, data = stirL, type = "l", log="x") > > y <- asNumeric( log(stilEM)) > x <- stirL$x > spl <- splinefun(log(x), y) > > plot(asNumeric(log(stilEM)) ~ x, data = stirL, type = "l", log="x") > lines(x, spl(log(x)), col=2) > summary(rE <- relErrV(target = y, current = spl(log(x)))) # all 0 {of course: interpolation at *these* x} Min. 1st Qu. Median Mean 3rd Qu. Max. 0 0 0 0 0 0 > > ssmpl <- smooth.spline(log(x), y) > str(ssmpl$fit)# only 174 knots, 170 coefficients List of 5 $ knot : num [1:125] 0 0 0 0 0.00783 ... $ nk : num 121 $ min : num -4.16 $ range: num 6.24 $ coef : num [1:121] -0.263 -0.266 -0.272 -0.282 -0.293 ... - attr(*, "class")= chr "smooth.spline.fit" > methods(class="smooth.spline.fit") #-> .. [1] predict see '?methods' for accessing help and source code > str(yp <- predict(ssmpl$fit, x=log(x))) List of 2 $ x: num [1:512] -4.16 -4.15 -4.13 -4.12 -4.11 ... $ y: num [1:512] -0.263 -0.265 -0.267 -0.27 -0.272 ... > lines(x, yp$y, col=adjustcolor(3, 1/3), lwd=3) # looks fine > ## but > summary(reSpl <- relErrV(target = y, current = yp$y)) Min. 1st Qu. Median Mean 3rd Qu. Max. -3.452e-09 -3.272e-10 -3.239e-12 -1.000e-16 3.819e-10 3.236e-09 > ## Min. 1st Qu. Median Mean 3rd Qu. Max. > ## -6.35e-10 -7.55e-11 1.10e-12 0.00e+00 7.32e-11 5.53e-10 > ## which is *NOT* good enough, of course .... > > > > showProc.time() Time (user system elapsed): 17.36 0.62 17.99 > > proc.time() user system elapsed 24.90 0.79 25.67