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 1) dbinom_raw(), dnbinomR() and dnbinom.mu() > #### 2) log1pmx(), logcf() etc > require(DPQ) Loading required package: DPQ > > source(system.file(package="DPQ", "test-tools.R", + mustWork=TRUE))# ../inst/test-tools.R > ## => showProc.time(), ... list_() , loadList() , readRDS_() , save2RDS() > options(warnPartialMatchArgs = FALSE) > > (doExtras <- DPQ:::doExtras() && !grepl("valgrind", R.home())) [1] FALSE > > do.pdf <- TRUE > do.pdf <- !dev.interactive(orNone = TRUE) > > ### 1. Testing dbinom_raw(), dnbinomR() and dnbinom.mu() >>> ../R/dbinom-nbinom.R <<< > ### ---------- ../man/dbinom_raw.Rd & ../man/dnbinomR.Rd > > eaxis <- sfsmisc::eaxis > relErrV <- sfsmisc::relErrV > ## "FIXME:" use relErrV() already here > > ### dbinom() vs dbinom.raw() : > > for(n in 1:20) { + cat("n=",n," ") + for(x in 0:n) + cat(".") + for(p in c(0, .1, .5, .8, 1)) { + stopifnot(all.equal(dbinom_raw(x, n, p, q=1-p, log=FALSE), + dbinom (x, n, p, log=FALSE)), + all.equal(dbinom_raw(x, n, p, q=1-p, log =TRUE), + dbinom (x, n, p, log =TRUE))) + } + cat("\n") + } n= 1 .. n= 2 ... n= 3 .... n= 4 ..... n= 5 ...... n= 6 ....... n= 7 ........ n= 8 ......... n= 9 .......... n= 10 ........... n= 11 ............ n= 12 ............. n= 13 .............. n= 14 ............... n= 15 ................ n= 16 ................. n= 17 .................. n= 18 ................... n= 19 .................... n= 20 ..................... > showProc.time() Time (user system elapsed): 0.03 0.02 0.05 > > ### dnbinom*() : > stopifnot(exprs = { + dnbinomR(0, 1, 1) == 1 + }) > > ### exploring 'eps' == "true" tests must be done with Rmpfr !! > > ### 2. Testing log1pmx(), logcf() etc > ### ---------- > > do.pdf [1] TRUE > if(do.pdf) { pdf.options(width = 9, height = 6.5)# for all {9/6.5 = 1.38 ; 4/3 < 1.38 < sqrt(2) [A4] + pdf("dnbinom-logcf.pdf") + } > > ### 2a: logcf() > ## == ======= > x <- c((-20:3)/4, (25:31)/32) # close (but not too close) to upper bound 1 > > (lC <- logcf (x, i=2, d=3, eps=1e-9)) [1] 0.2225364 0.2273165 0.2323964 0.2378089 0.2435921 0.2497910 0.2564582 [8] 0.2636564 0.2714610 0.2799634 0.2892758 0.2995378 0.3109259 0.3236668 [15] 0.3380584 0.3545014 0.3735507 0.3960035 0.4230578 0.4566237 0.5000000 [22] 0.5595850 0.6503112 0.8221318 0.8574109 0.8989257 0.9490572 1.0118259 [29] 1.0948318 1.2152882 1.4286636 > lCt <- logcf (x, i=2, d=3, eps=1e-9, trace=TRUE) ; stopifnot(identical(lCt, lC)) it= 0: ==> |b2|=162720|-> a2/b2=0.2282448377581121 it= 1: ==> |b2|=1.68458e+08|-> a2/b2=0.2235169038765654 it= 2: ==> |b2|=3.02689e+11|-> a2/b2=0.2227073061010817 it= 3: ==> |b2|=8.40216e+14|-> a2/b2=0.2225663427223296 it= 4: ==> |b2|=3.33607e+18|-> a2/b2=0.2225416527107165 it= 5: ==> |b2|=1.79478e+22|-> a2/b2=0.2225373164616655 it= 6: ==> |b2|=1.25703e+26|-> a2/b2=0.2225365537561161 it= 7: ==> |b2|=1.11146e+30|-> a2/b2=0.2225364194783633 it= 8: ==> |b2|=1.21086e+34|-> a2/b2=0.2225363958233605 it= 9: ==> |b2|=1.5936e+38|-> a2/b2=0.2225363916543374 it=10: ==> |b2|=2.49268e+42|-> a2/b2=0.2225363909193379 logcf(*) used 11 iterations. it= 0: ==> |b2|=151400|-> a2/b2=0.2325379788639366 it= 1: ==> |b2|=1.519e+08|-> a2/b2=0.2281774349397349 it= 2: ==> |b2|=2.64707e+11|-> a2/b2=0.2274604630468608 it= 3: ==> |b2|=7.12814e+14|-> a2/b2=0.2273407160828378 it= 4: ==> |b2|=2.74588e+18|-> a2/b2=0.22732060432467 it= 5: ==> |b2|=1.4333e+22|-> a2/b2=0.2273172178197734 it= 6: ==> |b2|=9.73998e+25|-> a2/b2=0.2273166467704412 it= 7: ==> |b2|=8.35605e+29|-> a2/b2=0.2273165503912295 it= 8: ==> |b2|=8.83286e+33|-> a2/b2=0.2273165341149803 it= 9: ==> |b2|=1.12795e+38|-> a2/b2=0.2273165313651214 it=10: ==> |b2|=1.71192e+42|-> a2/b2=0.2273165309003892 logcf(*) used 11 iterations. it= 0: ==> |b2|=140480|-> a2/b2=0.2371369589977221 it= 1: ==> |b2|=1.36437e+08|-> a2/b2=0.233144219132411 it= 2: ==> |b2|=2.30332e+11|-> a2/b2=0.2325159146623991 it= 3: ==> |b2|=6.0102e+14|-> a2/b2=0.2324155839716448 it= 4: ==> |b2|=2.24367e+18|-> a2/b2=0.232399478669337 it= 5: ==> |b2|=1.135e+22|-> a2/b2=0.2323968871196107 it= 6: ==> |b2|=7.47503e+25|-> a2/b2=0.2323964695381487 it= 7: ==> |b2|=6.21522e+29|-> a2/b2=0.2323964021949098 it= 8: ==> |b2|=6.3674e+33|-> a2/b2=0.2323963913282185 it= 9: ==> |b2|=7.88061e+37|-> a2/b2=0.2323963895740212 it=10: ==> |b2|=1.15921e+42|-> a2/b2=0.2323963892907578 logcf(*) used 11 iterations. it= 0: ==> |b2|=129960|-> a2/b2=0.2420764081255771 it= 1: ==> |b2|=1.22034e+08|-> a2/b2=0.2384505742154687 it= 2: ==> |b2|=1.99336e+11|-> a2/b2=0.2379065422904285 it= 3: ==> |b2|=5.03394e+14|-> a2/b2=0.2378237915834095 it= 4: ==> |b2|=1.81889e+18|-> a2/b2=0.2378111429951061 it= 5: ==> |b2|=8.90621e+21|-> a2/b2=0.2378092051966509 it= 6: ==> |b2|=5.67763e+25|-> a2/b2=0.2378089079362223 it= 7: ==> |b2|=4.56957e+29|-> a2/b2=0.237808862298923 it= 8: ==> |b2|=4.53158e+33|-> a2/b2=0.237808855288528 it= 9: ==> |b2|=5.429e+37|-> a2/b2=0.2378088542112312 logcf(*) used 10 iterations. it= 0: ==> |b2|=119840|-> a2/b2=0.2473965287049399 it= 1: ==> |b2|=1.08655e+08|-> a2/b2=0.2441351443683956 it= 2: ==> |b2|=1.71497e+11|-> a2/b2=0.2436705415694199 it= 3: ==> |b2|=4.18587e+14|-> a2/b2=0.243603512590089 it= 4: ==> |b2|=1.46194e+18|-> a2/b2=0.2435937980676084 it= 5: ==> |b2|=6.91963e+21|-> a2/b2=0.2435923870961549 it= 6: ==> |b2|=4.26415e+25|-> a2/b2=0.2435921819100178 it= 7: ==> |b2|=3.31759e+29|-> a2/b2=0.2435921520482697 it= 8: ==> |b2|=3.18042e+33|-> a2/b2=0.243592147700053 it= 9: ==> |b2|=3.68336e+37|-> a2/b2=0.2435921470666651 logcf(*) used 10 iterations. it= 0: ==> |b2|=110120|-> a2/b2=0.2531442971304032 it= 1: ==> |b2|=9.62637e+07|-> a2/b2=0.250243065925647 it= 2: ==> |b2|=1.46601e+11|-> a2/b2=0.2498525958064797 it= 3: ==> |b2|=3.45334e+14|-> a2/b2=0.2497994254950086 it= 4: ==> |b2|=1.16411e+18|-> a2/b2=0.2497921546097616 it= 5: ==> |b2|=5.31835e+21|-> a2/b2=0.249791158319238 it= 6: ==> |b2|=3.16349e+25|-> a2/b2=0.2497910216443144 it= 7: ==> |b2|=2.37577e+29|-> a2/b2=0.2497910028808943 it= 8: ==> |b2|=2.19845e+33|-> a2/b2=0.249791000303667 it= 9: ==> |b2|=2.45771e+37|-> a2/b2=0.2497909999495487 logcf(*) used 10 iterations. it= 0: ==> |b2|=100800|-> a2/b2=0.259375 it= 1: ==> |b2|=8.48232e+07|-> a2/b2=0.2568274658348188 it= 2: ==> |b2|=1.24442e+11|-> a2/b2=0.25650538512364 it= 3: ==> |b2|=2.82452e+14|-> a2/b2=0.2564642259697995 it= 4: ==> |b2|=9.17519e+17|-> a2/b2=0.2564589455665124 it= 5: ==> |b2|=4.03952e+21|-> a2/b2=0.2564582668445072 it= 6: ==> |b2|=2.3156e+25|-> a2/b2=0.2564581795087773 it= 7: ==> |b2|=1.67591e+29|-> a2/b2=0.256458168262858 it= 8: ==> |b2|=1.49457e+33|-> a2/b2=0.2564581668140725 logcf(*) used 9 iterations. it= 0: ==> |b2|=91880|-> a2/b2=0.266154222899434 it= 1: ==> |b2|=7.42974e+07|-> a2/b2=0.2639514124718479 it= 2: ==> |b2|=1.04819e+11|-> a2/b2=0.2636915518134123 it= 3: ==> |b2|=2.28837e+14|-> a2/b2=0.2636605954874903 it= 4: ==> |b2|=7.15064e+17|-> a2/b2=0.263656894453428 it= 5: ==> |b2|=3.02848e+21|-> a2/b2=0.2636564511866551 it= 6: ==> |b2|=1.67007e+25|-> a2/b2=0.263656398042764 it= 7: ==> |b2|=1.1628e+29|-> a2/b2=0.2636563916670788 it= 8: ==> |b2|=9.97611e+32|-> a2/b2=0.2636563909018438 logcf(*) used 9 iterations. it= 0: ==> |b2|=83360|-> a2/b2=0.2735604606525912 it= 1: ==> |b2|=6.46501e+07|-> a2/b2=0.2716904913342721 it= 2: ==> |b2|=8.75389e+10|-> a2/b2=0.2714862953286458 it= 3: ==> |b2|=1.83464e+14|-> a2/b2=0.271463799471145 it= 4: ==> |b2|=5.50387e+17|-> a2/b2=0.2714613129507076 it= 5: ==> |b2|=2.23803e+21|-> a2/b2=0.2714610376580362 it= 6: ==> |b2|=1.18496e+25|-> a2/b2=0.2714610071500224 it= 7: ==> |b2|=7.92152e+28|-> a2/b2=0.2714610037670337 it= 8: ==> |b2|=6.52535e+32|-> a2/b2=0.2714610033917415 logcf(*) used 9 iterations. it= 0: ==> |b2|=75240|-> a2/b2=0.2816885964912281 it= 1: ==> |b2|=5.58449e+07|-> a2/b2=0.2801362592848659 it= 2: ==> |b2|=7.24171e+10|-> a2/b2=0.2799808515150535 it= 3: ==> |b2|=1.45381e+14|-> a2/b2=0.2799651692452077 it= 4: ==> |b2|=4.17809e+17|-> a2/b2=0.2799635819710054 it= 5: ==> |b2|=1.6276e+21|-> a2/b2=0.2799634210721147 it= 6: ==> |b2|=8.25594e+24|-> a2/b2=0.2799634047475746 it= 7: ==> |b2|=5.28764e+28|-> a2/b2=0.279963403090364 logcf(*) used 8 iterations. it= 0: ==> |b2|=67520|-> a2/b2=0.2906546208530806 it= 1: ==> |b2|=4.78456e+07|-> a2/b2=0.2894009584998412 it= 2: ==> |b2|=5.92745e+10|-> a2/b2=0.2892872380428356 it= 3: ==> |b2|=1.13708e+14|-> a2/b2=0.2892768492714667 it= 4: ==> |b2|=3.12287e+17|-> a2/b2=0.2892758976324906 it= 5: ==> |b2|=1.16261e+21|-> a2/b2=0.2892758103387463 it= 6: ==> |b2|=5.6361e+24|-> a2/b2=0.2892758023247274 it= 7: ==> |b2|=3.44989e+28|-> a2/b2=0.2892758015886094 logcf(*) used 8 iterations. it= 0: ==> |b2|=60200|-> a2/b2=0.3006021594684385 it= 1: ==> |b2|=4.06158e+07|-> a2/b2=0.2996240768690056 it= 2: ==> |b2|=4.79397e+10|-> a2/b2=0.2995448531150074 it= 3: ==> |b2|=8.76351e+13|-> a2/b2=0.2995383960829794 it= 4: ==> |b2|=2.2937e+17|-> a2/b2=0.2995378685241684 it= 5: ==> |b2|=8.13827e+20|-> a2/b2=0.2995378253666954 it= 6: ==> |b2|=3.76013e+24|-> a2/b2=0.299537821833512 it= 7: ==> |b2|=2.19363e+28|-> a2/b2=0.29953782154412 logcf(*) used 8 iterations. it= 0: ==> |b2|=53280|-> a2/b2=0.3117117117117117 it= 1: ==> |b2|=3.41194e+07|-> a2/b2=0.3109816831265299 it= 2: ==> |b2|=3.82483e+10|-> a2/b2=0.3109298549498925 it= 3: ==> |b2|=6.64186e+13|-> a2/b2=0.3109261555416012 it= 4: ==> |b2|=1.6515e+17|-> a2/b2=0.3109258909118994 it= 5: ==> |b2|=5.56707e+20|-> a2/b2=0.3109258719607546 it= 6: ==> |b2|=2.44378e+24|-> a2/b2=0.3109258706026756 logcf(*) used 7 iterations. it= 0: ==> |b2|=46760|-> a2/b2=0.3242140718562874 it= 1: ==> |b2|=2.83198e+07|-> a2/b2=0.3237000517835029 it= 2: ==> |b2|=3.0043e+10|-> a2/b2=0.3236688340788258 it= 3: ==> |b2|=4.93794e+13|-> a2/b2=0.32366692941967 it= 4: ==> |b2|=1.16224e+17|-> a2/b2=0.3236668129926885 it= 5: ==> |b2|=3.70875e+20|-> a2/b2=0.3236668058687138 it= 6: ==> |b2|=1.54119e+24|-> a2/b2=0.3236668054325496 logcf(*) used 7 iterations. it= 0: ==> |b2|=40640|-> a2/b2=0.3384104330708662 it= 1: ==> |b2|=2.3181e+07|-> a2/b2=0.3380761409363547 it= 2: ==> |b2|=2.31738e+10|-> a2/b2=0.3380593343638449 it= 3: ==> |b2|=3.59e+13|-> a2/b2=0.3380584861718837 it= 4: ==> |b2|=7.96488e+16|-> a2/b2=0.3380584432964664 it= 5: ==> |b2|=2.39588e+20|-> a2/b2=0.3380584411272919 logcf(*) used 6 iterations. it= 0: ==> |b2|=34920|-> a2/b2=0.3547036082474227 it= 1: ==> |b2|=1.86664e+07|-> a2/b2=0.3545094440828331 it= 2: ==> |b2|=1.74976e+10|-> a2/b2=0.3545017297423328 it= 3: ==> |b2|=2.5422e+13|-> a2/b2=0.3545014222983448 it= 4: ==> |b2|=5.29017e+16|-> a2/b2=0.3545014100293588 it= 5: ==> |b2|=1.49263e+20|-> a2/b2=0.3545014095393998 logcf(*) used 6 iterations. it= 0: ==> |b2|=29600|-> a2/b2=0.3736486486486487 it= 1: ==> |b2|=1.474e+07|-> a2/b2=0.3735535956580733 it= 2: ==> |b2|=1.28785e+10|-> a2/b2=0.3735508120722061 it= 3: ==> |b2|=1.74436e+13|-> a2/b2=0.3735507303649829 it= 4: ==> |b2|=3.38438e+16|-> a2/b2=0.3735507279641467 logcf(*) used 5 iterations. it= 0: ==> |b2|=24680|-> a2/b2=0.3960393030794165 it= 1: ==> |b2|=1.13653e+07|-> a2/b2=0.3960042169989204 it= 2: ==> |b2|=9.18785e+09|-> a2/b2=0.3960035415346633 it= 3: ==> |b2|=1.1517e+13|-> a2/b2=0.3960035285101282 it= 4: ==> |b2|=2.06815e+16|-> a2/b2=0.3960035282588064 logcf(*) used 5 iterations. it= 0: ==> |b2|=20160|-> a2/b2=0.4230654761904762 it= 1: ==> |b2|=8.50608e+06|-> a2/b2=0.4230579185711867 it= 2: ==> |b2|=6.30386e+09|-> a2/b2=0.4230578415832492 it= 3: ==> |b2|=7.24564e+12|-> a2/b2=0.423057840798356 logcf(*) used 4 iterations. it= 0: ==> |b2|=16040|-> a2/b2=0.4566240648379052 it= 1: ==> |b2|=6.12601e+06|-> a2/b2=0.4566236526711513 it= 2: ==> |b2|=4.11202e+09|-> a2/b2=0.45662365139084 logcf(*) used 3 iterations. logcf(*) used 0 iterations. it= 0: ==> |b2|=9000|-> a2/b2=0.5595833333333333 it= 1: ==> |b2|=2.65815e+06|-> a2/b2=0.5595850262400541 it= 2: ==> |b2|=1.38218e+09|-> a2/b2=0.5595850350182284 logcf(*) used 3 iterations. it= 0: ==> |b2|=6080|-> a2/b2=0.6501644736842105 it= 1: ==> |b2|=1.49776e+06|-> a2/b2=0.6503067914752697 it= 2: ==> |b2|=6.50656e+08|-> a2/b2=0.6503110275633139 it= 3: ==> |b2|=4.39124e+11|-> a2/b2=0.65031115303776 it= 4: ==> |b2|=4.24985e+14|-> a2/b2=0.6503111567461622 logcf(*) used 5 iterations. it= 0: ==> |b2|=3560|-> a2/b2=0.8174859550561798 it= 1: ==> |b2|=671330|-> a2/b2=0.8216038498205056 it= 2: ==> |b2|=2.24237e+08|-> a2/b2=0.8220723977911059 it= 3: ==> |b2|=1.16565e+11|-> a2/b2=0.8221251256996892 it= 4: ==> |b2|=8.69636e+13|-> a2/b2=0.8221310310708708 it= 5: ==> |b2|=8.80714e+16|-> a2/b2=0.8221316908161992 it= 6: ==> |b2|=1.16246e+20|-> a2/b2=0.8221317644149767 it= 7: ==> |b2|=1.93847e+23|-> a2/b2=0.8221317726176806 it= 8: ==> |b2|=3.98491e+26|-> a2/b2=0.8221317735312994 logcf(*) used 9 iterations. it= 0: ==> |b2|=3273.12|-> a2/b2=0.8501825949971358 it= 1: ==> |b2|=589700|-> a2/b2=0.8564359269653014 it= 2: ==> |b2|=1.88377e+08|-> a2/b2=0.8572808958391158 it= 3: ==> |b2|=9.36959e+10|-> a2/b2=0.8573936483488037 it= 4: ==> |b2|=6.68994e+13|-> a2/b2=0.8574086115763199 it= 5: ==> |b2|=6.48488e+16|-> a2/b2=0.8574105917266195 it= 6: ==> |b2|=8.19327e+19|-> a2/b2=0.8574108533369799 it= 7: ==> |b2|=1.30789e+23|-> a2/b2=0.8574108878636231 it= 8: ==> |b2|=2.57381e+26|-> a2/b2=0.8574108924170954 it= 9: ==> |b2|=6.12129e+29|-> a2/b2=0.857410893017314 logcf(*) used 10 iterations. it= 0: ==> |b2|=2992.5|-> a2/b2=0.887515664160401 it= 1: ==> |b2|=512650|-> a2/b2=0.8970879760307886 it= 2: ==> |b2|=1.55894e+08|-> a2/b2=0.8986335219604555 it= 3: ==> |b2|=7.3859e+10|-> a2/b2=0.8988795255023756 it= 4: ==> |b2|=5.02475e+13|-> a2/b2=0.8989184315726793 it= 5: ==> |b2|=4.64164e+16|-> a2/b2=0.8989245645338572 it= 6: ==> |b2|=5.58911e+19|-> a2/b2=0.8989255294740479 it= 7: ==> |b2|=8.50347e+22|-> a2/b2=0.8989256811124874 it= 8: ==> |b2|=1.595e+26|-> a2/b2=0.8989257049228231 it= 9: ==> |b2|=3.61574e+29|-> a2/b2=0.8989257086593775 it=10: ==> |b2|=9.74479e+32|-> a2/b2=0.8989257092455057 logcf(*) used 11 iterations. it= 0: ==> |b2|=2718.12|-> a2/b2=0.9306636583122557 it= 1: ==> |b2|=440109|-> a2/b2=0.9454886448454234 it= 2: ==> |b2|=1.26644e+08|-> a2/b2=0.9483747936978182 it= 3: ==> |b2|=5.68225e+10|-> a2/b2=0.9489276102289144 it= 4: ==> |b2|=3.66244e+13|-> a2/b2=0.9490326945801036 it= 5: ==> |b2|=3.20598e+16|-> a2/b2=0.9490525921046072 it= 6: ==> |b2|=3.65864e+19|-> a2/b2=0.9490563512518312 it= 7: ==> |b2|=5.27587e+22|-> a2/b2=0.9490570604554582 it= 8: ==> |b2|=9.37997e+25|-> a2/b2=0.949057194128513 it= 9: ==> |b2|=2.01557e+29|-> a2/b2=0.9490572193069129 it=10: ==> |b2|=5.14924e+32|-> a2/b2=0.9490572240471546 it=11: ==> |b2|=1.54257e+36|-> a2/b2=0.9490572249392525 logcf(*) used 12 iterations. it= 0: ==> |b2|=2450|-> a2/b2=0.98125 it= 1: ==> |b2|=372006|-> a2/b2=1.004588821592379 it= 2: ==> |b2|=1.00485e+08|-> a2/b2=1.010139849506342 it= 3: ==> |b2|=4.23633e+10|-> a2/b2=1.011436268076868 it= 4: ==> |b2|=2.56713e+13|-> a2/b2=1.011736272952714 it= 5: ==> |b2|=2.11343e+16|-> a2/b2=1.011805363744281 it= 6: ==> |b2|=2.26869e+19|-> a2/b2=1.011821231674799 it= 7: ==> |b2|=3.07772e+22|-> a2/b2=1.011824869836296 it= 8: ==> |b2|=5.14811e+25|-> a2/b2=1.011825703044044 it= 9: ==> |b2|=1.04082e+29|-> a2/b2=1.011825893713502 it=10: ==> |b2|=2.50192e+32|-> a2/b2=1.011825937320711 it=11: ==> |b2|=7.05238e+35|-> a2/b2=1.011825947289596 it=12: ==> |b2|=2.30384e+39|-> a2/b2=1.01182594956778 logcf(*) used 13 iterations. it= 0: ==> |b2|=2188.12|-> a2/b2=1.041575621251071 it= 1: ==> |b2|=308271|-> a2/b2=1.079158530385997 it= 2: ==> |b2|=7.72745e+07|-> a2/b2=1.090294707272438 it= 3: ==> |b2|=3.02658e+10|-> a2/b2=1.093530698375819 it= 4: ==> |b2|=1.70531e+13|-> a2/b2=1.094460628360649 it= 5: ==> |b2|=1.30605e+16|-> a2/b2=1.09472621968612 it= 6: ==> |b2|=1.30466e+19|-> a2/b2=1.094801803325279 it= 7: ==> |b2|=1.64734e+22|-> a2/b2=1.094823266186136 it= 8: ==> |b2|=2.56499e+25|-> a2/b2=1.09482935204936 it= 9: ==> |b2|=4.82765e+28|-> a2/b2=1.094831075998267 it=10: ==> |b2|=1.08039e+32|-> a2/b2=1.094831563992047 it=11: ==> |b2|=2.83535e+35|-> a2/b2=1.094831702052934 it=12: ==> |b2|=8.62389e+38|-> a2/b2=1.094831741096325 it=13: ==> |b2|=3.00926e+42|-> a2/b2=1.094831752134144 it=14: ==> |b2|=1.19409e+46|-> a2/b2=1.094831755253794 it=15: ==> |b2|=5.34632e+49|-> a2/b2=1.094831756135323 logcf(*) used 16 iterations. it= 0: ==> |b2|=1932.5|-> a2/b2=1.115014553686934 it= 1: ==> |b2|=248832|-> a2/b2=1.177472194314061 it= 2: ==> |b2|=5.68734e+07|-> a2/b2=1.201237997708952 it= 3: ==> |b2|=2.03226e+10|-> a2/b2=1.210122180173716 it= 4: ==> |b2|=1.04577e+13|-> a2/b2=1.213401220746793 it= 5: ==> |b2|=7.32086e+15|-> a2/b2=1.214601770463929 it= 6: ==> |b2|=6.68834e+18|-> a2/b2=1.215039143200027 it= 7: ==> |b2|=7.72653e+21|-> a2/b2=1.215197983065948 it= 8: ==> |b2|=1.10096e+25|-> a2/b2=1.21525555017665 it= 9: ==> |b2|=1.89662e+28|-> a2/b2=1.215276384510363 it=10: ==> |b2|=3.88536e+31|-> a2/b2=1.215283917220463 it=11: ==> |b2|=9.33474e+34|-> a2/b2=1.215286638689748 it=12: ==> |b2|=2.59938e+38|-> a2/b2=1.215287621371808 it=13: ==> |b2|=8.30457e+41|-> a2/b2=1.21528797604937 it=14: ==> |b2|=3.01718e+45|-> a2/b2=1.215288104018202 it=15: ==> |b2|=1.23692e+49|-> a2/b2=1.215288150176859 it=16: ==> |b2|=5.68258e+52|-> a2/b2=1.21528816682257 it=17: ==> |b2|=2.90768e+56|-> a2/b2=1.21528817282419 it=18: ==> |b2|=1.64796e+60|-> a2/b2=1.21528817498773 logcf(*) used 19 iterations. it= 0: ==> |b2|=1683.12|-> a2/b2=1.206723449684367 it= 1: ==> |b2|=193619|-> a2/b2=1.315303466880831 it= 2: ==> |b2|=3.91439e+07|-> a2/b2=1.371129848795821 it= 3: ==> |b2|=1.23338e+10|-> a2/b2=1.399699792974485 it= 4: ==> |b2|=5.59551e+12|-> a2/b2=1.414183204867562 it= 5: ==> |b2|=3.4562e+15|-> a2/b2=1.421462092550596 it= 6: ==> |b2|=2.78868e+18|-> a2/b2=1.425095727190543 it= 7: ==> |b2|=2.84748e+21|-> a2/b2=1.426900785436191 it= 8: ==> |b2|=3.58854e+24|-> a2/b2=1.427794344145924 it= 9: ==> |b2|=5.4701e+27|-> a2/b2=1.42823558059338 it=10: ==> |b2|=9.91885e+30|-> a2/b2=1.428453069803035 it=11: ==> |b2|=2.10987e+34|-> a2/b2=1.428560130489993 it=12: ==> |b2|=5.20269e+37|-> a2/b2=1.428612779766881 it=13: ==> |b2|=1.47211e+41|-> a2/b2=1.428638651529776 it=14: ==> |b2|=4.73732e+44|-> a2/b2=1.428651357361294 it=15: ==> |b2|=1.72036e+48|-> a2/b2=1.428657594367876 it=16: ==> |b2|=7.00164e+51|-> a2/b2=1.428660654811599 it=17: ==> |b2|=3.17394e+55|-> a2/b2=1.428662156076175 it=18: ==> |b2|=1.59374e+59|-> a2/b2=1.428662892312951 it=19: ==> |b2|=8.8205e+62|-> a2/b2=1.428663253292892 it=20: ==> |b2|=5.35623e+66|-> a2/b2=1.428663430250051 it=21: ==> |b2|=3.5541e+70|-> a2/b2=1.428663516983012 it=22: ==> |b2|=2.56724e+74|-> a2/b2=1.428663559488062 it=23: ==> |b2|=2.01172e+78 Lrg |b2||-> a2/b2=1.428663580315935 it=24: ==> |b2|=147221|-> a2/b2=1.428663590520718 it=25: ==> |b2|=1.34508e+09|-> a2/b2=1.428663595520171 it=26: ==> |b2|=1.32142e+13|-> a2/b2=1.428663597969267 it=27: ==> |b2|=1.39232e+17|-> a2/b2=1.428663599168924 logcf(*) used 28 iterations. > (lR <- logcfR(x, i=2, d=3, eps=1e-9)) [1] 0.2225364 0.2273165 0.2323964 0.2378089 0.2435921 0.2497910 0.2564582 [8] 0.2636564 0.2714610 0.2799634 0.2892758 0.2995378 0.3109259 0.3236668 [15] 0.3380584 0.3545014 0.3735507 0.3960035 0.4230578 0.4566237 0.5000000 [22] 0.5595850 0.6503112 0.8221318 0.8574109 0.8989257 0.9490572 1.0118259 [29] 1.0948318 1.2152882 1.4286636 > all.equal(lC, lR, tol = 0) # x86_64 F40: 6.54e-16 [1] "Mean relative difference: 6.541027e-16" > stopifnot(all.equal(lC, lR, tol = 4e-15)) # 1.08295e-15 (Apple clang 14.0.3) > lRt <- logcfR(x, i=2, d=3, eps=1e-9, trace=TRUE) ; stopifnot(identical(lRt, lR)) logcf(x[], i=2, d=3, eps=1e-09) iterations: it= 1: needIt: 30 TRUE, and 1 F.; length(x[])=30, m.B2= 2.33e+04 it= 2: needIt: 30 TRUE; length(x[])=30, m.B2= 1.01e+07 it= 3: needIt: 30 TRUE; length(x[])=30, m.B2= 7.71e+09 it= 4: needIt: 28 TRUE, and 2 F.; length(x[])=28, m.B2= 1.01e+13 it= 5: needIt: 27 TRUE, and 1 F.; length(x[])=27, m.B2= 1.76e+16 it= 6: needIt: 24 TRUE, and 3 F.; length(x[])=24, m.B2= 4.75e+19 it= 7: needIt: 22 TRUE, and 2 F.; length(x[])=22, m.B2= 1.28e+23 it= 8: needIt: 20 TRUE, and 2 F.; length(x[])=20, m.B2= 3.64e+26 it= 9: needIt: 17 TRUE, and 3 F.; length(x[])=17, m.B2= 6.87e+29 it=10: needIt: 13 TRUE, and 4 F.; length(x[])=13, m.B2= 1.04e+33 it=11: needIt: 9 TRUE, and 4 F.; length(x[])= 9, m.B2= 3.09e+35 it=12: needIt: 5 TRUE, and 4 F.; length(x[])= 5, m.B2= 2.27e+35 it=13: needIt: 4 TRUE, and 1 F.; length(x[])= 4, m.B2= 4.05e+38 it=14: needIt: 3 TRUE, and 1 F.; todo: x[29,30,31]= c(0.906, 0.938, 0.969), m.B2= 7.17e+41 it=15: needIt: 3 TRUE; todo: x[29,30,31]= c(0.906, 0.938, 0.969), m.B2= 2.57e+45 it=16: needIt: 3 TRUE; todo: x[29,30,31]= c(0.906, 0.938, 0.969), m.B2= 1.04e+49 it=17: needIt: 2 TRUE, and 1 F.; todo: x[30,31]= c(0.938, 0.969), m.B2= 1.99e+52 it=18: needIt: 2 TRUE; todo: x[30,31]= c(0.938, 0.969), m.B2= 9.61e+55 it=19: needIt: 2 TRUE; todo: x[30,31]= c(0.938, 0.969), m.B2= 5.12e+59 it=20: needIt: 1 TRUE, and 1 F.; todo: x[31]= c(0.969), m.B2= 8.82e+62 it=21: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 5.36e+66 it=22: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 3.55e+70 it=23: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 2.57e+74 it=24: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 2.01e+78 Lrg it=25: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.47e+05 it=26: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.35e+09 it=27: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.32e+13 it=28: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.39e+17 logcf(*) end: after 28 iterations. > lRt2 <- logcfR(x, i=2, d=3, eps=1e-9, trace= 2) ; stopifnot(identical(lRt2,lR)) logcf(x[], i=2, d=3, eps=1e-09) iterations: it= 1: needIt: 30 TRUE, and 1 F.; length(x[])=30, m.B2= 2.33e+04 it= 2: needIt: 30 TRUE; length(x[])=30, m.B2= 1.01e+07 it= 3: needIt: 30 TRUE; length(x[])=30, m.B2= 7.71e+09 it= 4: needIt: 28 TRUE, and 2 F.; length(x[])=28, m.B2= 1.01e+13 it= 5: needIt: 27 TRUE, and 1 F.; length(x[])=27, m.B2= 1.76e+16 it= 6: needIt: 24 TRUE, and 3 F.; length(x[])=24, m.B2= 4.75e+19 it= 7: needIt: 22 TRUE, and 2 F.; length(x[])=22, m.B2= 1.28e+23 it= 8: needIt: 20 TRUE, and 2 F.; length(x[])=20, m.B2= 3.64e+26 it= 9: needIt: 17 TRUE, and 3 F.; length(x[])=17, m.B2= 6.87e+29 it=10: needIt: 13 TRUE, and 4 F.; length(x[])=13, m.B2= 1.04e+33 it=11: needIt: 9 TRUE, and 4 F.; length(x[])= 9, m.B2= 3.09e+35 it=12: needIt: 5 TRUE, and 4 F.; length(x[])= 5, m.B2= 2.27e+35 it=13: needIt: 4 TRUE, and 1 F.; length(x[])= 4, m.B2= 4.05e+38 it=14: needIt: 3 TRUE, and 1 F.; todo: x[29,30,31]= c(0.906, 0.938, 0.969), m.B2= 7.17e+41 it=15: needIt: 3 TRUE; todo: x[29,30,31]= c(0.906, 0.938, 0.969), m.B2= 2.57e+45 it=16: needIt: 3 TRUE; todo: x[29,30,31]= c(0.906, 0.938, 0.969), m.B2= 1.04e+49 it=17: needIt: 2 TRUE, and 1 F.; todo: x[30,31]= c(0.938, 0.969), m.B2= 1.99e+52 it=18: needIt: 2 TRUE; todo: x[30,31]= c(0.938, 0.969), m.B2= 9.61e+55 it=19: needIt: 2 TRUE; todo: x[30,31]= c(0.938, 0.969), m.B2= 5.12e+59 it=20: needIt: 1 TRUE, and 1 F.; todo: x[31]= c(0.969), m.B2= 8.82e+62 it=21: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 5.36e+66 it=22: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 3.55e+70 it=23: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 2.57e+74 it=24: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 2.01e+78 Lrg it=25: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.47e+05 it=26: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.35e+09 it=27: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.32e+13 it=28: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.39e+17 logcf(*) end: after 28 iterations. > > lR. <- logcfR(x, i=2, d=3, eps=1e-9) > lR.t <- logcfR(x, i=2, d=3, eps=1e-9, trace=TRUE) ; stopifnot(identical(lR.t, lR.)) logcf(x[], i=2, d=3, eps=1e-09) iterations: it= 1: needIt: 30 TRUE, and 1 F.; length(x[])=30, m.B2= 2.33e+04 it= 2: needIt: 30 TRUE; length(x[])=30, m.B2= 1.01e+07 it= 3: needIt: 30 TRUE; length(x[])=30, m.B2= 7.71e+09 it= 4: needIt: 28 TRUE, and 2 F.; length(x[])=28, m.B2= 1.01e+13 it= 5: needIt: 27 TRUE, and 1 F.; length(x[])=27, m.B2= 1.76e+16 it= 6: needIt: 24 TRUE, and 3 F.; length(x[])=24, m.B2= 4.75e+19 it= 7: needIt: 22 TRUE, and 2 F.; length(x[])=22, m.B2= 1.28e+23 it= 8: needIt: 20 TRUE, and 2 F.; length(x[])=20, m.B2= 3.64e+26 it= 9: needIt: 17 TRUE, and 3 F.; length(x[])=17, m.B2= 6.87e+29 it=10: needIt: 13 TRUE, and 4 F.; length(x[])=13, m.B2= 1.04e+33 it=11: needIt: 9 TRUE, and 4 F.; length(x[])= 9, m.B2= 3.09e+35 it=12: needIt: 5 TRUE, and 4 F.; length(x[])= 5, m.B2= 2.27e+35 it=13: needIt: 4 TRUE, and 1 F.; length(x[])= 4, m.B2= 4.05e+38 it=14: needIt: 3 TRUE, and 1 F.; todo: x[29,30,31]= c(0.906, 0.938, 0.969), m.B2= 7.17e+41 it=15: needIt: 3 TRUE; todo: x[29,30,31]= c(0.906, 0.938, 0.969), m.B2= 2.57e+45 it=16: needIt: 3 TRUE; todo: x[29,30,31]= c(0.906, 0.938, 0.969), m.B2= 1.04e+49 it=17: needIt: 2 TRUE, and 1 F.; todo: x[30,31]= c(0.938, 0.969), m.B2= 1.99e+52 it=18: needIt: 2 TRUE; todo: x[30,31]= c(0.938, 0.969), m.B2= 9.61e+55 it=19: needIt: 2 TRUE; todo: x[30,31]= c(0.938, 0.969), m.B2= 5.12e+59 it=20: needIt: 1 TRUE, and 1 F.; todo: x[31]= c(0.969), m.B2= 8.82e+62 it=21: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 5.36e+66 it=22: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 3.55e+70 it=23: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 2.57e+74 it=24: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 2.01e+78 Lrg it=25: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.47e+05 it=26: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.35e+09 it=27: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.32e+13 it=28: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.39e+17 logcf(*) end: after 28 iterations. > > all.equal(lC, lR., tol = 0) # no longer TRUE, but really small [1] "Mean relative difference: 6.541027e-16" > all.equal(lR, lR., tol = 0) # TRUE !! " " [1] TRUE > stopifnot(all.equal(lC, lR., tol = 1e-14)) > ## (even though they used eps=1e-9 .. i.e., are not *so* accurate) > lR.14 <- logcf(x, i=2, d=3, eps=1e-14)# default when used in R, incl. from log1pmx(): > lR.18 <- logcf(x, i=2, d=3, eps=1e-18) > > showProc.time() Time (user system elapsed): 0.08 0 0.08 > > > ##require(Rmpfr) may be not, see if NS loading (via "::") is sufficient: > requireNamespace("Rmpfr") || quit("no") Loading required namespace: Rmpfr [1] TRUE > ## ----- ---------- > asNumeric <- Rmpfr::asNumeric > all.equal <- Rmpfr::all.equal > mpfr <- Rmpfr::mpfr > getPrec <- Rmpfr::getPrec > .getPrec <- Rmpfr::.getPrec > > xM <- mpfr(x, 512) > (ct.24 <- system.time(lRM24 <- logcfR (xM, i=2, d=3, eps=1e-24, trace=TRUE))) # it=76, 0.73 sec logcf(x[], i=2, d=3, eps=1e-24) iterations: it= 1: needIt: 30 TRUE, and 1 F.; length(x[])=30, m.B2= 2.33e+04 it= 2: needIt: 30 TRUE; length(x[])=30, m.B2= 1.01e+07 it= 3: needIt: 30 TRUE; length(x[])=30, m.B2= 7.71e+09 it= 4: needIt: 30 TRUE; length(x[])=30, m.B2= 9.11e+12 it= 5: needIt: 30 TRUE; length(x[])=30, m.B2= 1.54e+16 it= 6: needIt: 30 TRUE; length(x[])=30, m.B2= 3.55e+19 it= 7: needIt: 30 TRUE; length(x[])=30, m.B2= 1.06e+23 it= 8: needIt: 30 TRUE; length(x[])=30, m.B2= 4.01e+26 it= 9: needIt: 30 TRUE; length(x[])=30, m.B2= 1.87e+30 it=10: needIt: 29 TRUE, and 1 F.; length(x[])=29, m.B2= 1.11e+34 it=11: needIt: 28 TRUE, and 1 F.; length(x[])=28, m.B2= 8.86e+37 it=12: needIt: 27 TRUE, and 1 F.; length(x[])=27, m.B2= 7.27e+41 it=13: needIt: 27 TRUE; length(x[])=27, m.B2= 6.75e+45 it=14: needIt: 26 TRUE, and 1 F.; length(x[])=26, m.B2= 7.04e+49 it=15: needIt: 26 TRUE; length(x[])=26, m.B2= 8.51e+53 it=16: needIt: 24 TRUE, and 2 F.; length(x[])=24, m.B2= 1.63e+58 it=17: needIt: 24 TRUE; length(x[])=24, m.B2= 2.54e+62 it=18: needIt: 23 TRUE, and 1 F.; length(x[])=23, m.B2= 3.83e+66 it=19: needIt: 22 TRUE, and 1 F.; length(x[])=22, m.B2= 5.89e+70 it=20: needIt: 21 TRUE, and 1 F.; length(x[])=21, m.B2= 8.99e+74 it=21: needIt: 20 TRUE, and 1 F.; length(x[])=20, m.B2= 1.33e+79 Lrg it=22: needIt: 19 TRUE, and 1 F.; length(x[])=19, m.B2= 1.61e+06 it=23: needIt: 18 TRUE, and 1 F.; length(x[])=18, m.B2= 2.05e+10 it=24: needIt: 17 TRUE, and 1 F.; length(x[])=17, m.B2= 2.3e+14 it=25: needIt: 16 TRUE, and 1 F.; length(x[])=16, m.B2= 2.16e+18 it=26: needIt: 14 TRUE, and 2 F.; length(x[])=14, m.B2= 5.47e+22 it=27: needIt: 13 TRUE, and 1 F.; length(x[])=13, m.B2= 3.02e+26 it=28: needIt: 11 TRUE, and 2 F.; length(x[])=11, m.B2= 5.02e+30 it=29: needIt: 10 TRUE, and 1 F.; length(x[])=10, m.B2= 9.27e+33 it=30: needIt: 8 TRUE, and 2 F.; length(x[])= 8, m.B2= 4.27e+37 it=31: needIt: 6 TRUE, and 2 F.; length(x[])= 6, m.B2= 3.71e+36 it=32: needIt: 5 TRUE, and 1 F.; length(x[])= 5, m.B2= 1.6e+36 it=33: needIt: 4 TRUE, and 1 F.; length(x[])= 4, m.B2= 7.87e+39 it=34: needIt: 4 TRUE; length(x[])= 4, m.B2= 1.41e+44 it=35: needIt: 4 TRUE; length(x[])= 4, m.B2= 2.66e+48 it=36: needIt: 4 TRUE; length(x[])= 4, m.B2= 5.3e+52 it=37: needIt: 4 TRUE; length(x[])= 4, m.B2= 1.12e+57 it=38: needIt: 3 TRUE, and 1 F.; todo: x[29,30,31]= c(0.906, 0.938, 0.969), m.B2= 4.78e+60 it=39: needIt: 3 TRUE; todo: x[29,30,31]= c(0.906, 0.938, 0.969), m.B2= 1.07e+65 it=40: needIt: 3 TRUE; todo: x[29,30,31]= c(0.906, 0.938, 0.969), m.B2= 2.51e+69 it=41: needIt: 3 TRUE; todo: x[29,30,31]= c(0.906, 0.938, 0.969), m.B2= 6.17e+73 it=42: needIt: 3 TRUE; todo: x[29,30,31]= c(0.906, 0.938, 0.969), m.B2= 1.59e+78 Lrg it=43: needIt: 3 TRUE; todo: x[29,30,31]= c(0.906, 0.938, 0.969), m.B2= 3.72e+05 it=44: needIt: 2 TRUE, and 1 F.; todo: x[30,31]= c(0.938, 0.969), m.B2= 1.16e+09 it=45: needIt: 2 TRUE; todo: x[30,31]= c(0.938, 0.969), m.B2= 3.27e+13 it=46: needIt: 2 TRUE; todo: x[30,31]= c(0.938, 0.969), m.B2= 9.59e+17 it=47: needIt: 2 TRUE; todo: x[30,31]= c(0.938, 0.969), m.B2= 2.93e+22 it=48: needIt: 2 TRUE; todo: x[30,31]= c(0.938, 0.969), m.B2= 9.36e+26 it=49: needIt: 2 TRUE; todo: x[30,31]= c(0.938, 0.969), m.B2= 3.11e+31 it=50: needIt: 2 TRUE; todo: x[30,31]= c(0.938, 0.969), m.B2= 1.07e+36 it=51: needIt: 2 TRUE; todo: x[30,31]= c(0.938, 0.969), m.B2= 3.85e+40 it=52: needIt: 2 TRUE; todo: x[30,31]= c(0.938, 0.969), m.B2= 1.44e+45 it=53: needIt: 2 TRUE; todo: x[30,31]= c(0.938, 0.969), m.B2= 5.56e+49 it=54: needIt: 1 TRUE, and 1 F.; todo: x[31]= c(0.969), m.B2= 8.39e+52 it=55: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 3.29e+57 it=56: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.33e+62 it=57: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 5.61e+66 it=58: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 2.44e+71 it=59: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.1e+76 it=60: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 5.11e+80 Lrg it=61: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 2.12e+08 it=62: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.05e+13 it=63: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 5.38e+17 it=64: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 2.84e+22 it=65: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.54e+27 it=66: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 8.66e+31 it=67: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 5e+36 it=68: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 2.98e+41 it=69: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.82e+46 it=70: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.15e+51 it=71: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 7.43e+55 it=72: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 4.94e+60 it=73: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 3.38e+65 it=74: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 2.38e+70 it=75: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.71e+75 it=76: needIt: 1 TRUE; todo: x[31]= c(0.969), m.B2= 1.27e+80 Lrg logcf(*) end: after 76 iterations. user system elapsed 0.72 0.03 0.75 > if(doExtras) withAutoprint({ # ----------------------------------- + ct24 <- system.time(lR24 <- logcfR_vec(xM, i=2, d=3, eps=1e-24, trace=TRUE)); ct24 # 5.7 sec + all.equal(lRM24, lR24, tol=0) # TRUE + identical(lRM24, lR24) # TRUE !! (not sure if on all platforms!) + stopifnot(all.equal(lRM24, lR24, tol = 5e-16)) + }) > > ## ===> use logcfR() {the internally vectorized version from now on) > > SS <- function(ch, digits = 4) sub(paste0("([0-9]{1,",digits,"})[0-9]*e"), "\\1e", ch) > ## double prec <--> MPFR: vvvv (same eps) > lRM9 <- logcfR(xM, 2,3, eps=1e-9) > ## show: > SS(all.equal(Rmpfr::roundMpfr(lRM9, 64), lR, tol=0))# .. 5.1138e-16 [1] "Mean relative difference: 5.1138e-16" > stopifnot(all.equal(lRM9, lR , tol=1e-15)) > SS(all.equal(lRM24, lRM9, tol=0))# .. 3.701..e-10 [1] "Mean relative difference: 3.7013e-10" > stopifnot(all.equal(lRM24, lRM9, tol=1e-9)) > showProc.time() Time (user system elapsed): 1.3 0.06 1.36 > ## now see if small eps makes a relevant difference: > relE14 <- asNumeric(relErrV(lRM24, lR.14)) > relE18 <- asNumeric(relErrV(lRM24, lR.18)) > summary(cbind(relE14, relE18)) relE14 relE18 Min. :-8.766e-15 Min. :-9.169e-17 1st Qu.:-9.971e-16 1st Qu.:-3.435e-17 Median : 1.268e-15 Median :-7.028e-20 Mean : 1.345e-15 Mean :-3.644e-18 3rd Qu.: 3.843e-15 3rd Qu.: 3.514e-17 Max. : 1.090e-14 Max. : 7.480e-17 > ## relE14 relE18 > ## Min. :-8.766e-15 Min. :-9.169e-17 > ## 1st Qu.:-9.971e-16 1st Qu.:-3.435e-17 > ## Median : 1.268e-15 Median :-7.028e-20 > ## Mean : 1.345e-15 Mean :-3.644e-18 > ## 3rd Qu.: 3.843e-15 3rd Qu.: 3.514e-17 > ## Max. : 1.090e-14 Max. : 7.480e-17 > matplot(x, pmax(abs(cbind(relE14, relE18)), 2^-55), log = "y", type = "l", ylim = c(3e-17, 1e-14), + panel.first = abline(h= 2^-(54:51), lty=3, lwd = c(1,1,3,1), col="lightgray")) > title("|relative Errors| (wrt mpfr precBits = 512) of logcf(x, 2,3, eps = *)") > legend("left", paste("eps = ", format(c(1e-14, 1e-18))), col=1:2, lty=1:2) > > ## zoom into [0, 1) -- which is the domain of y = (x/(x+2)^2 for log1pmx(x) : > N <- if(doExtras) 1024 else 128 > x <- (0:(N-1))/N # close (but not too close) to upper bound 1 > xM <- mpfr(x, 512) > system.time(lgRM <- logcfR(xM, i=2, d=3, eps=1e-24, trace=doExtras))# it=152; 1.5 s {doEx.. it=425, 5.6 s} user system elapsed 1.45 0.02 1.47 > relE14 <- asNumeric(relErrV(lgRM, logcf(x, i=2, d=3, eps=1e-14))) > relE18 <- asNumeric(relErrV(lgRM, logcf(x, i=2, d=3, eps=1e-18))) > > summary(cbind(relE14, relE18)) relE14 relE18 Min. :-2.335e-14 Min. :-9.503e-17 1st Qu.:-2.334e-15 1st Qu.:-4.085e-17 Median :-7.081e-16 Median :-1.371e-17 Mean :-1.984e-15 Mean :-6.597e-18 3rd Qu.:-8.660e-17 3rd Qu.: 3.087e-17 Max. : 7.918e-17 Max. : 1.003e-16 > ## relE14 relE18 > ## Min. :-5.935e-14 Min. :-1.087e-16 > ## 1st Qu.:-2.565e-15 1st Qu.:-3.831e-17 > ## Median :-6.919e-16 Median :-1.261e-18 > ## Mean :-2.274e-15 Mean : 9.056e-19 > ## 3rd Qu.:-1.144e-16 3rd Qu.: 4.275e-17 > ## Max. : 9.903e-17 Max. : 1.155e-16 > matplot(x, pmax(abs(cbind(relE14, relE18)), 2^-55), log = "y", type = "l", ylim = c(3e-17, 4e-14), + ylab = "", main = "logcf(x, 2,3, eps = *): |relative Errors| wrt mpfr precBits = 512", + panel.first = abline(h= 2^-(54:51), lty=3, lwd = c(1,1,3,1), col="lightgray"), yaxt="n") > eaxis(2) > legend("topleft", paste("eps = ", format(c(1e-14, 1e-18))), col=1:2, lty=1:2, bty = "n") > ## Wow !!! ===> definitely should use a smaller eps or tol_logcf !! > > if(do.pdf) { dev.off(); pdf("dnbinom-log1pmx.pdf") } > > ### 2b: log1pmx(x) -- calls logcf(y, 2,3), for y = t^2 = (x/(x+2))^2 in [0,1) for x > -1 > ## == ========= > > ## The first part is from original ../man/log1pmx.Rd > > e <- if(doExtras) 2^-12 else 2^-8; by.p <- 1/(if(doExtras) 256 else 64) > xd <- sort(c(seq(-1+e, 0+100*e, by=e), seq(by.p, 5, by=by.p))) # length 676 or 5476 if do.X. > plot(xd, log1pmx(xd), type="l", col=2, lwd=2, main = "log1pmx(x)") > abline(h=0, v=-1:0, lty=3) > > xM <- mpfr(xd, 512) > ## for MPFR numbers, really better than tol_logcf = eps = 1e-14 (default): > if(doExtras) { + print( system.time( + lg1pM <- log1pmx(xM, tol_logcf = 1e-25, eps2 = 1e-4) + )) # 2.3 sec [new logcfR()] + lines(xd, asNumeric(lg1pM), col=adjustcolor(4, 1/4), lwd=3) + } > lg1pdR <- log1pmx(xd, trace=TRUE, logCF=logcfR) logcf(x[], i=3, d=2, eps=1e-14) iterations: it= 1: needIt: 362 TRUE; length(x[])=362, m.B2= 9.41e+03 it= 2: needIt: 321 TRUE, and 41 F.; length(x[])=321, m.B2= 1.73e+06 it= 3: needIt: 191 TRUE, and 130 F.; length(x[])=191, m.B2= 4.72e+08 it= 4: needIt: 106 TRUE, and 85 F.; length(x[])=106, m.B2= 1.7e+11 it= 5: needIt: 56 TRUE, and 50 F.; length(x[])=56, m.B2= 7.52e+13 it= 6: needIt: 31 TRUE, and 25 F.; length(x[])=31, m.B2= 4.35e+16 it= 7: needIt: 13 TRUE, and 18 F.; length(x[])=13, m.B2= 3.06e+19 logcf(*) end: after 7 iterations. > lines(xd, lg1pdR, col=adjustcolor(5, 1/4), lwd=5) > ## > lg1pM. <- log1p(xM) - xM # good enough if(precBits are high enough) > if(doExtras) { + print(summary(relEM <- asNumeric(relErrV(lg1pM, lg1pM.)))) + ## Min. 1st Qu. Median Mean 3rd Qu. Max. + ## -7.910e-28 0.000e+00 5.600e-33 2.282e-28 2.732e-29 1.104e-26 + stopifnot(abs(relEM) < 1e-25) + ## the error of the log1pmx() "algorithm" even for "perfect" mpfr-accuracy: + xM2k <- mpfr(xd, 2048) # even more bits + lg1pM2k <- log1p(xM2k) - xM2k + reE00 <- asNumeric( relErrV(lg1pM2k, lg1pM.) ) + print(signif(rrre <- range(reE00)), 2) # [-1.5e-151, 4.8e-151] ==> 512 bits is sufficient + stopifnot(abs(rrre) < 1e-150) + plot(reE00 ~ xd, type="l") # see nothing but close to x ~= 0 + plot(abs(reE00) ~ xd, type="l", log="y") + ## zoom in close to x ~= 0: + plot(reE00 ~ xd, type="l", subset = abs(xd) <= 1/64) + plot(abs(reE00) ~ xd, type="l", log="y", subset = abs(xd) <= 1) + rE.log1pm <- asNumeric(relErrV(lg1pM2k, lg1pM)) + print(summary( rE.log1pm ) ) # -1.1e-26 7.9e-28 {"matching" tol_logcf = 1e-25 above} + } > showProc.time() Time (user system elapsed): 1.54 0.03 1.57 > re <- asNumeric(relErrV(lg1pM., log1pmx(xd))) > > ## MM: From around here, "move" to vignette <<<<<<<<<<<<<<<<<<<<<<< ../vignettes/log1pmx-etc.Rnw <<< > > plot(xd, re, type="b", cex=1/4, xlab = quote(x), + main = paste("relErrV(log1p(xM) - xM, log1pmx(x)), xM <- mpfr(x,",min(.getPrec(xM)),")")) > abline(h = (-2:2)*2^-52, lty=2, col=adjustcolor("gray20", 1/2)) > xl <- -0.84; xr <- -0.4 ## the "relevant x-range": > colR <- adjustcolor("tomato", 1/2) > abline(v= c(xl,xr), col=colR, lty=2, lwd=2) > axis(1,at=c(xl,xr), labels=NULL, col.axis="tomato", col=colR, lwd=2, cex.axis=3/4) > text(xl, adj=c(0,-1/2), par("usr")[3], "zoom in", col=colR) > showProc.time() Time (user system elapsed): 0.02 0 0.02 > > # only "relevant" x-range: values inside (xl, xr) > iN <- xl < xd & xd < xr > x.iN <- xd[iN] > plot(x.iN, re[iN], type="b", cex=1/4, main="rel.Error of log1pmx(x)", + xlab=quote(x), ylab=quote(rE(x)), sub = "(in relevant x-range)") > abline(h = (-2:2)*2^-52, lty=2, col=adjustcolor("gray30", 1/2)) > abline(v = c(xl,xr), col=colR, lty=2, lwd=2) > > ## *absolute* relative errors from here on: > plot(x.iN, abs(re[iN]), type="b", cex=1/2, log="y", + main = "| relErr( log1pmx(x) ) | {via 'Rmpfr'}", + ylab=quote(abs(rE(x))), xlab=quote(x), + ylim = c(4e-17, max(abs(asNumeric(re)[iN])))) > abline(h = c(1:2,4)*2^-52, lty=2, col=adjustcolor("gray20", 1/2)) > mL1 <- eval(formals(log1pmx)$minL1) > abline(v = mL1, lwd=3, col=adjustcolor(2, 1/2)) > axis(3, at=mL1, sprintf("minL1=%7.4f",mL1), col.axis=2, col=2) > > mL1 <- -0.7; re2 <- asNumeric(relErrV(lg1pM.[iN], log1pmx(x.iN, minL1 = mL1))) > cat("at x=", x.iN[im <- which.max(abs(re2))], "rel.Err.=", format(re2[im], digits = 3),"\n") at x= -0.6679688 rel.Err.= -5.37e-16 > ## at x= -0.6677246 rel.Err.= -6.48e-16 > lines(x.iN, abs(re2), col=adjustcolor(4, 1/2), lwd=2) > abline(v = mL1, lwd=3, col=adjustcolor(4, 2/3), lty=3) > axis(3, line=-1, at=mL1, sprintf("mL1 = %g",mL1), col.axis=4, col=4) > R <- x.iN >= mL1; stopifnot(re2[R] == re[iN][R]) > > mL1 <- -0.66; re3 <- asNumeric(relErrV(lg1pM.[iN], log1pmx(x.iN, minL1 = mL1))) > re3.2 <- asNumeric(relErrV(lg1pM.[iN], log1pmx(x.iN, minL1 = mL1, eps2 = .02))) > print(max(abs(re3))) == # 3.7029e-16 + max(abs(re3.2)) # TRUE ! of course, eps2 does not matter as long as it is far from our x-range !! [1] 3.165487e-16 [1] TRUE > cat("at x=", x.iN[im <- which.max(abs(re3))], "rel.Err.=", format(re3[im], digits = 3),"\n") at x= -0.65625 rel.Err.= -3.17e-16 > ## at x= -0.5634766 rel.Err.= -3.70e-16 > lines(x.iN, abs(re3), col=adjustcolor(6, 1/3), lwd=2) > abline(v = mL1, lwd=3, col=adjustcolor(6, 2/3), lty=3) > axis(3, line=-2, at=mL1, sprintf("mL1 = %g",mL1), col.axis=6, col=6) > lines(lowess(x.iN, abs(re[iN]), f=1/50), col=adjustcolor("gray", 1/2), lwd=6) > lines(lowess(x.iN, abs(re2), f=1/50), col=adjustcolor(4, 1/2), lwd=6) > lines(lowess(x.iN, abs(re3), f=1/50), col=adjustcolor(6, 1/2), lwd=6) > R <- x.iN > mL1; stopifnot(re3[R] == re2[R]) > ## MM: for vignette, stop here? > ## ---> since max(|re4|) {below} == max(|re3|) > > mL1 <- -0.64; re4 <- asNumeric(relErrV(lg1pM.[iN], log1pmx(x.iN, minL1 = mL1, tol_logcf = 3e-16))) > max(abs(re4)) == max(abs(re3)) # TRUE ! i.e., are *not* better (in minimax-sense) [1] FALSE > cat("at x=", x.iN[im <- which.max(abs(re4))], "rel.Err.=", format(re4[im], digits = 3),"\n") at x= -0.6445312 rel.Err.= 2.75e-16 > ## at x= -0.5634766 rel.Err.= -3.7e-16 -- of course unchanged from re3 (-0.66) > lines(x.iN, abs(re4), col=adjustcolor(7, 1/3), lwd=2) > abline(v = mL1, lwd=3, col=adjustcolor(7, 2/3), lty=3) > axis(3, line=-2, at=mL1, sprintf("mL1 = %g",mL1), col.axis=7, col=7) > lines(lowess(x.iN, abs(re[iN]), f=1/50), col=adjustcolor("gray", 1/2), lwd=6) > lines(lowess(x.iN, abs(re2), f=1/50), col=adjustcolor(4, 1/2), lwd=6) > lines(lowess(x.iN, abs(re3), f=1/50), col=adjustcolor(6, 1/2), lwd=6) > lines(lowess(x.iN, abs(re4), f=1/50), col=adjustcolor(7, 1/2), lwd=6) > R <- x.iN > mL1; table(re4[R] == re3[R]) # no longer FALSE TRUE 5 56 > > mL1 <- -0.6; re5 <- asNumeric(relErrV(lg1pM.[iN], log1pmx(x.iN, minL1 = mL1))) > max(abs(re5)) == max(abs(re3)) # TRUE (as still < -0.564.. !) [1] FALSE > cat("at x=", x.iN[im <- which.max(abs(re5))], "rel.Err.=", format(re5[im], digits = 3),"\n") at x= -0.5703125 rel.Err.= -3.06e-16 > lines(x.iN, abs(re5), col=adjustcolor(8, 1/3), lwd=2) > abline(v = mL1, lwd=3, col=adjustcolor(8, 2/3), lty=3) > axis(3, line=-1, at=mL1, sprintf("mL1 = %g",mL1), col.axis=8, col=8) > lines(lowess(x.iN, abs(re[iN]), f=1/50), col=adjustcolor("gray", 1/2), lwd=6) > lines(lowess(x.iN, abs(re2), f=1/50), col=adjustcolor(4, 1/2), lwd=6) > lines(lowess(x.iN, abs(re3), f=1/50), col=adjustcolor(6, 1/2), lwd=6) > lines(lowess(x.iN, abs(re4), f=1/50), col=adjustcolor(7, 1/2), lwd=6) > lines(lowess(x.iN, abs(re5), f=1/50), col=adjustcolor(8, 1/2), lwd=6) > R <- x.iN > mL1; table(re5[R] == re4[R]) # no longer equal FALSE TRUE 4 47 > ## -0.6 now is clearly too large, -0.7 was better {MM: why the > showProc.time() Time (user system elapsed): 0.05 0 0.05 > > > ##--- separate: is eps2 = 0.01 "optimal" / "good enough" ?? > > k. <- if(doExtras) 3 else 0 > str(xS <- unique(sort(c(seq(-1/32, 1/32, by=2^-(12+k.)), + seq(-1,1, by=2^-(7+k.))/2^(4+k.))))) ## 257 or 3841 values num [1:385] -0.0625 -0.062 -0.0615 -0.061 -0.0605 ... > plot(xS, log1pmx(xS), type="l", col=2, main = "log1pmx(x)") > abline(h=0, v = c(-.01, .01), lwd=3, col=adjustcolor(7, 2/3), lty=3) > ## nothing visible at +- .01 > > xSm <- mpfr(xS, 1024) > xSM <- mpfr(xS, 8192) > lg1pM. <- log1p(xSm) - xSm ## the direct formally using many bits ==> not much cancellation > summary(abs(asNumeric(relErrV(log1p(xSM) - xSM, lg1pM.))))# <= 4.1e-304 : 1024 is enough! Min. 1st Qu. Median Mean 3rd Qu. Max. 0.000e+00 6.609e-308 1.345e-307 4.574e-307 3.077e-307 1.695e-305 > reS <- asNumeric(relErrV(lg1pM., log1pmx(xS))) > summary(reS) ## Lnx 64b gcc 11.2.1, 20210728: Min. 1st Qu. Median Mean 3rd Qu. Max. -1.924e-16 -4.766e-17 4.682e-18 4.446e-18 5.684e-17 2.299e-16 > ## Min. 1st Qu. Median Mean 3rd Qu. Max. > ## -2.625e-16 -5.246e-17 1.973e-18 1.211e-18 5.286e-17 2.688e-16 > plot(xS, abs(reS), type="b", cex=1/2, log="y", + main = "| relErr( log1pmx() ) | {via 'Rmpfr'}" + , ylim = c(4e-17, max(abs(asNumeric(reS)))) + , panel.last= abline(h=0, v = c(-.01, .01), lwd=3, col=adjustcolor(7, 2/3), lty=3)) Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > ## no problem visible at +/- eps2 = +/- 0.01 > stopifnot(abs(reS) < 1e-15, mean(abs(reS)) < 6e-16) # seen 2.688e-16, 6.35e-17 > showProc.time() Time (user system elapsed): 0.2 0 0.2 > > ## in spite of the above; ... what's the approximation error of accurate log1pmx() ? > ## vvvvv (or 1e-18 above ?) > p.Err <- function(eps2 = .01, tol_logcf = 1e-17, xM=xSm) { + stopifnot(eps2 > 0, tol_logcf > 0, + inherits(xM, "mpfr"), (prec <- .getPrec(xM)) >= 64) + lg1pM <- log1p(xM) - xM # direct formula suffering from cancellation + lg1pM.<- log1pmx(xM, tol_logcf=tol_logcf, eps2=eps2) + sys.call() + reM <- asNumeric(relErrV(lg1pM., lg1pM)) + plot(asNumeric(xM), abs(reM), type="b", cex=1/2, log="y", xlab = quote(x), + ylim = quantile(abs(reM), c(1/20,1), names=FALSE), + sub = paste("min(getPrec(x)) =", min(prec)), + main = sprintf("relative error of R log1pmx(eps2=%g, tol_logcf=%g)", eps2, tol_logcf)) + abline(h=0, v = eps2*c(-1,1), lwd=3, col=adjustcolor(7, 2/3), lty=3) + invisible(reM) + } > > ### --- for tol_logcf = 1e-17 --------------------------------------------- > reSM <- p.Err() # aha! Here, eps2 = .01 was definitely too large! Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > reSM001 <- p.Err(.001) # still too large Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > xM2 <- mpfr(seq(-1,1, length=1+1024)/2^8, 1024) > reSM3_4 <- p.Err(3e-4, xM=xM2)# still too large Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > reSM294 <- p.Err(2.9e-4, xM=xM2)# still Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > reSM2854<- p.Err(2.85e-4, xM=xM2)# still (and practically identical to 2.9e-4) Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > ## cbind(reSM294, reSM2854)[reSM294 != reSM2854 , ] > ## reSM294 reSM2854 > ## [1,] 2.572818e-36 2.988004e-46 > ## [2,] -2.565617e-36 -2.977915e-46 > reSM284 <- p.Err(2.8e-4, xM=xM2)# now too small -- bad! Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > reSM274 <- p.Err(2.7e-4, xM=xM2)# (ditto) too small Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > showProc.time() Time (user system elapsed): 1.98 0 1.98 > > ### --- for tol_logcf = 1e-14 == R's default! > reSM <- p.Err( tol = 1e-14) ## eps2=.01 is *still* too large Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > reSM002 <- p.Err(.002, tol = 1e-14) # still too large Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > reSM0017 <- p.Err(.0017, tol = 1e-14) # still .. Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > reSM00165<- p.Err(.00165,tol = 1e-14) # still.. Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > reSM00163<- p.Err(.00163,tol = 1e-14) # still.. <<<<<<<<<<<<< SHOULD WE CHANGE R TO USE THIS? Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > reSM00162<- p.Err(.00162,tol = 1e-14) # already (very barely!) too small Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > reSM0016 <- p.Err(.0016, tol = 1e-14) # already (barely!) too small Warning message: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > showProc.time() Time (user system elapsed): 0.97 0 0.97 > if(do.pdf) dev.off() null device 1 > > > summary(warnings()) 1 identical warnings: In xy.coords(x, y, xlabel, ylabel, log) : 1 y value <= 0 omitted from logarithmic plot > > proc.time() user system elapsed 6.29 0.15 6.43