CRAN Package Check Results for Package DPQ

Last updated on 2019-11-26 00:51:50 CET.

Flavor Version Tinstall Tcheck Ttotal Status Flags
r-devel-linux-x86_64-debian-clang 0.3-5 9.40 147.95 157.35 OK
r-devel-linux-x86_64-debian-gcc 0.3-5 7.17 110.52 117.69 OK
r-devel-linux-x86_64-fedora-clang 0.3-5 183.97 OK
r-devel-linux-x86_64-fedora-gcc 0.3-5 180.49 OK
r-devel-windows-ix86+x86_64 0.3-5 32.00 361.00 393.00 OK
r-devel-windows-ix86+x86_64-gcc8 0.3-5 30.00 277.00 307.00 OK
r-patched-linux-x86_64 0.3-5 8.26 143.20 151.46 OK
r-patched-solaris-x86 0.3-5 245.90 OK
r-release-linux-x86_64 0.3-5 9.25 319.13 328.38 OK
r-release-windows-ix86+x86_64 0.3-5 35.00 296.00 331.00 OK
r-release-osx-x86_64 0.3-5 OK
r-oldrel-windows-ix86+x86_64 0.3-5 18.00 361.00 379.00 ERROR
r-oldrel-osx-x86_64 0.3-5 ERROR

Check Details

Version: 0.3-5
Check: whether package can be installed
Result: WARN
    Found the following significant warnings:
     Warning: unable to re-encode 'pnbeta.R' lines 12, 35
     Warning: unable to re-encode 'pnchisq.R' lines 263, 264, 265, 266, 648, 652, 654, 667, 675, 676, 716
     Warning: unable to re-encode 'qnchisq.R' lines 135, 136, 146, 147, 148, 149, 243, 248, 250
Flavor: r-oldrel-windows-ix86+x86_64

Version: 0.3-5
Check: running tests for arch ‘i386’
Result: ERROR
     Running 'chisq-nonc-ex.R' [34s]
     Running 'dnchisq-tst.R' [1s]
     Running 'pnbeta-tst.R' [1s]
     Running 'pnt-precision-problem.R' [21s]
     Running 'ppois-ex.R' [3s]
     Running 'qPoisBinom-ex.R' [1s]
     Running 'qbeta-dist.R' [18s]
     Running 'qbeta-tst.R' [1s]
     Running 'qgamma-ex.R' [23s]
     Running 't-nonc-tst.R' [9s]
     Running 'wienergerm-accuracy.R' [9s]
     Running 'wienergerm-pchisq-tst.R' [1s]
     Running 'wienergerm_nchisq.R' [3s]
    Running the tests in 'tests/qgamma-ex.R' failed.
    Complete output:
     > library(DPQ)
     >
     > ###---> Automatically find places where qgamma() is not so precise (PR#2214) :
     > ### For PR#2214, had '1e-8' below and found quite a bit
     > ## see /u/maechler/R/MM/NUMERICS/dpq-functions/beta-gamma-etc/qgamma-ex.R ..
     >
     > ## FIXME: Timing ! --- partly these matplot() partly get quite slow ~?
     > source(system.file(package="Matrix", "test-tools-1.R", mustWork=TRUE))
     Loading required package: tools
     > ##--> showProc.time(), assertError(), relErrV(), ...
     > showProc.time()
     Time elapsed: 1.61 0.06 1.68
     >
     > (doExtras <- DPQ:::doExtras())
     [1] FALSE
     > (sdir <- system.file("safe", package="DPQ")) ## save directory (to read from)
     [1] "D:/temp/RtmpIpG5h1/RLIBS_cc0c300314ea/DPQ/safe"
     >
     > ### Nowadays finds cases in a special region for really small p and cutoff 1e-11 :
     > set.seed(47)
     > n <- if(doExtras) 100 else 32
     > res <- cbind(p=1,df=1,rE=1)[-1,]
     > for(M in 1:(if(doExtras) 20 else 10))
     + for(p in runif(n)) for(df in rlnorm(n)) {
     + r <- 1- pchisq(qchisq(p, df),df)/p
     + if(abs(r) > 1e-11) res <- rbind(res, c(p,df,r))
     + }
     >
     > ### use df in U[0,1]: finds two cases with bound 1e-11
     > for(p in runif(n)/2) for(df in runif(n)) {
     + qq <- qchisq(p, df)
     + if(qq > 0 && p > 0) {
     + r <- 1- pchisq(qq, df) / p
     + if(abs(r) > 1e-11) res <- rbind(res, c(p,df,r))
     + }
     + }
     >
     > ### now df very close to 0 : ==> finds more cases
     > for(p in sort(c(runif(64)/2, exp(-(1+rlnorm(256))))))
     + for(df in 2^-rlnorm(256, mean=2, sdlog=1.5)) {
     + qq <- qchisq(p, df)
     + if(qq > 0 && p > 0) {
     + r <- 1- pchisq(qq, df) / p
     + if(abs(r) > 1e-11) res <- rbind(res, c(p,df,r))
     + }
     + }
     > showProc.time()
     Time elapsed: 0.78 0 0.78
     >
     > require(graphics)
     > if(!dev.interactive(orNone=TRUE)) pdf("qgamma-appr.pdf")
     > eaxis <- sfsmisc::eaxis
     >
     > showProc.time()
     Time elapsed: 0.09 0 0.1
     > ## if(nrow(res) > 0) {
     > cat("Found inaccurate examples where pchisq(qchisq(p, df),df) != p\n")
     Found inaccurate examples where pchisq(qchisq(p, df),df) != p
     > ## sort in p, then df:
     > res <- res[order(res[,"p"], res[,"df"]), ]
     > rE <- res[,"rE"]
     > if(nrow(res) > 20) { hist(rE, breaks = 30); rug(rE) }
     > plot(res[,1:2])##--> quite interesting : all along one curve
     > ## p <= 1/2 and df <= 1 (about) !!
     > res <- cbind(res, nDig = round(-log10(abs(rE)), 1))
     > print(res, digits=12)
     p df rE nDig
     [1,] 0.000194375438651 0.02334079639198 6.46592989151e-08 7.2
     [2,] 0.000605300028912 0.02041606754775 -1.99820160418e-11 10.7
     [3,] 0.001012316063255 0.01855615147677 -2.59145555105e-04 3.6
     [4,] 0.001248285290785 0.01838201076117 3.94217658517e-10 9.4
     [5,] 0.001682388899865 0.01720736646288 5.53088974600e-04 3.3
     [6,] 0.001746787400790 0.01731189518997 -5.86897856980e-08 7.2
     [7,] 0.002664451237518 0.01599317398629 1.48421342013e-04 3.8
     [8,] 0.002664451237518 0.01618024201222 7.75564700239e-08 7.1
     [9,] 0.003159421860255 0.01557612780310 -7.92117005632e-06 5.1
     [10,] 0.003159421860255 0.01568183691729 -4.52237560733e-08 7.3
     [11,] 0.004055462418244 0.01493858731306 -7.17404703909e-06 5.1
     [12,] 0.004400694140827 0.01459101672970 9.07907026435e-04 3.0
     [13,] 0.004458811277768 0.01457506850867 9.03139988525e-05 4.0
     [14,] 0.004481882165743 0.01468883074316 -3.23309491179e-07 6.5
     [15,] 0.004939609905705 0.01440168350452 -2.81810098879e-06 5.6
     [16,] 0.008824465120182 0.01276352706510 1.21107345756e-04 3.9
     [17,] 0.009040265960535 0.01273711629661 1.38964402733e-05 4.9
     [18,] 0.010839089634828 0.01242499920422 2.63408184153e-10 9.6
     [19,] 0.011642124851282 0.01201471267173 -3.00271327153e-04 3.5
     [20,] 0.014753716559535 0.01155624353203 1.52961088240e-10 9.8
     [21,] 0.015499213434879 0.01125420134457 1.00492447767e-04 4.0
     [22,] 0.015499213434879 0.01135920381800 -9.55739012376e-08 7.0
     [23,] 0.018603016576955 0.01071716109330 -2.07537936111e-03 2.7
     [24,] 0.018603016576955 0.01073655493589 2.14388784341e-04 3.7
     [25,] 0.022624242394389 0.01033379525113 -3.37865668776e-09 8.5
     [26,] 0.022624242394389 0.01034206121729 4.43079739565e-08 7.4
     [27,] 0.023730217356634 0.01016252135853 -1.07732682664e-06 6.0
     [28,] 0.032427027472295 0.00942923095016 5.11205522358e-11 10.3
     [29,] 0.044753525441333 0.00839626444749 1.22224173549e-05 4.9
     [30,] 0.081818424963746 0.00686007746204 -1.78633419168e-09 8.7
     [31,] 0.081818424963746 0.00689856335721 -2.30915286892e-11 10.6
     [32,] 0.082800309102258 0.00681234719059 4.17997558788e-09 8.4
     [33,] 0.083507718914457 0.00680676700443 9.77167236016e-11 10.0
     [34,] 0.090821658072474 0.00655269761981 -7.16033632386e-09 8.1
     [35,] 0.102294760453517 0.00623563107239 -6.63889809793e-09 8.2
     [36,] 0.110869751789691 0.00603268830251 8.95578944338e-10 9.0
     [37,] 0.123950804624116 0.00571305309327 2.84683721041e-10 9.5
     [38,] 0.127405857731893 0.00562369059572 6.60541454867e-09 8.2
     [39,] 0.135229634154169 0.00540073357520 -2.34762594200e-05 4.6
     [40,] 0.137732279982451 0.00533092076413 2.99285844990e-04 3.5
     [41,] 0.138112917548194 0.00535138710974 -2.05335777981e-06 5.7
     [42,] 0.141100635980184 0.00527305771429 4.31593832968e-05 4.4
     [43,] 0.141100635980184 0.00537073537183 8.96455243371e-10 9.0
     [44,] 0.142905299416015 0.00523680041306 3.48180824883e-04 3.5
     [45,] 0.145624557210331 0.00526923971034 -1.94501770245e-09 8.7
     [46,] 0.154606872884529 0.00506806894407 -4.59924667240e-07 6.3
     [47,] 0.154606872884529 0.00507366168703 2.72301046933e-07 6.6
     [48,] 0.163535630067488 0.00497650928578 3.39664962823e-11 10.5
     [49,] 0.169741036539408 0.00484181845356 5.31400978776e-09 8.3
     [50,] 0.177327576288650 0.00465956102839 5.53404362603e-05 4.3
     [51,] 0.178169157856761 0.00471949961255 4.79807193976e-10 9.3
     [52,] 0.190094017358772 0.00450373552308 1.78565421871e-06 5.7
     [53,] 0.190147641510530 0.00453468705710 5.66235636157e-09 8.2
     [54,] 0.200112534472267 0.00442273120514 7.20473680715e-11 10.1
     [55,] 0.201518808589718 0.00439936964342 1.58750790291e-11 10.8
     [56,] 0.201518808589718 0.00439976887947 -9.97180116258e-11 10.0
     [57,] 0.210803673024037 0.00427351441034 -1.70232716812e-10 9.8
     [58,] 0.213058614771766 0.00426179831847 -1.39808165045e-11 10.9
     [59,] 0.214780951412088 0.00419869272965 -1.91879370171e-08 7.7
     [60,] 0.232805106603566 0.00395399315002 -9.17581020055e-08 7.0
     [61,] 0.249102914025652 0.00380019404026 -1.15818687974e-10 9.9
     [62,] 0.249102914025652 0.00382493512126 -1.39672717836e-11 10.9
     [63,] 0.252076511947811 0.00374903834738 2.14569900181e-07 6.7
     [64,] 0.253082914021191 0.00375259362798 3.65436092498e-09 8.4
     [65,] 0.253922058700076 0.00371237348323 -5.09014937222e-06 5.3
     [66,] 0.254289278570932 0.00374343873151 -1.05664899053e-09 9.0
     [67,] 0.260017499519858 0.00366179605930 -2.90430199001e-07 6.5
     [68,] 0.270323906831467 0.00351999192121 -1.56164756277e-04 3.8
     [69,] 0.271699356057456 0.00355068132680 5.13092990317e-09 8.3
     [70,] 0.275516196070002 0.00346804047756 -4.35171547588e-04 3.4
     [71,] 0.280722231049885 0.00348224101220 -6.07109695849e-10 9.2
     [72,] 0.284601233201101 0.00344936339590 -2.63026489478e-10 9.6
     [73,] 0.290188543054775 0.00336613521112 -5.64443074502e-08 7.2
     [74,] 0.290579022038283 0.00334423496113 1.02667567781e-07 7.0
     [75,] 0.290579022038283 0.00336764858994 2.26061565023e-08 7.6
     [76,] 0.291850198713803 0.00333552811650 -1.27338760580e-06 5.9
     [77,] 0.296521136452775 0.00330308865102 -4.48927431229e-07 6.3
     [78,] 0.298034174946132 0.00330462333485 8.42470393447e-09 8.1
     [79,] 0.300556783277253 0.00323922530004 4.66003314391e-05 4.3
     [80,] 0.303182283998467 0.00328704590597 1.82253101499e-11 10.7
     [81,] 0.303182283998467 0.00328723648952 -2.23643326080e-11 10.7
     [82,] 0.322319846303892 0.00306134512927 -1.15130830540e-05 4.9
     [83,] 0.322319846303892 0.00310689001755 8.57751647487e-11 10.1
     [84,] 0.325071272052651 0.00302343293053 -2.47088704493e-04 3.6
     [85,] 0.325071272052651 0.00304146419577 -7.71376615361e-06 5.1
     [86,] 0.331888412404218 0.00300837121343 -4.96098895297e-09 8.3
     [87,] 0.362278153188527 0.00278204202032 4.53939330569e-10 9.3
     [88,] 0.385389476781711 0.00260981704384 -1.77430203863e-09 8.8
     [89,] 0.425333956955001 0.00232995789362 1.82823025607e-08 7.7
     [90,] 0.439503709203564 0.00222452690840 -4.53585193982e-06 5.3
     [91,] 0.439503709203564 0.00224964327069 -3.02331937263e-10 9.5
     [92,] 0.450804624124430 0.00216770324934 -4.59455036239e-08 7.3
     >
     > if(requireNamespace("scatterplot3d")) {
     + scatterplot3d::scatterplot3d(res[,1:3], type ='h') ## quite interesting:
     + ## the inaccurate (p,df) points are on nice monotone curve !!!
     + ## this is *less* revealing
     + scatterplot3d::scatterplot3d(res[,c("p","df","nDig")], type ='h')
     + }
     Loading required namespace: scatterplot3d
     > rL <- res[abs(res[,'rE']) > 1e-9,]
     > rL <- rL[order(rL[,1],rL[,2]),]
     > rL
     p df rE nDig
     [1,] 0.0001943754 0.023340796 6.465930e-08 7.2
     [2,] 0.0010123161 0.018556151 -2.591456e-04 3.6
     [3,] 0.0016823889 0.017207366 5.530890e-04 3.3
     [4,] 0.0017467874 0.017311895 -5.868979e-08 7.2
     [5,] 0.0026644512 0.015993174 1.484213e-04 3.8
     [6,] 0.0026644512 0.016180242 7.755647e-08 7.1
     [7,] 0.0031594219 0.015576128 -7.921170e-06 5.1
     [8,] 0.0031594219 0.015681837 -4.522376e-08 7.3
     [9,] 0.0040554624 0.014938587 -7.174047e-06 5.1
     [10,] 0.0044006941 0.014591017 9.079070e-04 3.0
     [11,] 0.0044588113 0.014575069 9.031400e-05 4.0
     [12,] 0.0044818822 0.014688831 -3.233095e-07 6.5
     [13,] 0.0049396099 0.014401684 -2.818101e-06 5.6
     [14,] 0.0088244651 0.012763527 1.211073e-04 3.9
     [15,] 0.0090402660 0.012737116 1.389644e-05 4.9
     [16,] 0.0116421249 0.012014713 -3.002713e-04 3.5
     [17,] 0.0154992134 0.011254201 1.004924e-04 4.0
     [18,] 0.0154992134 0.011359204 -9.557390e-08 7.0
     [19,] 0.0186030166 0.010717161 -2.075379e-03 2.7
     [20,] 0.0186030166 0.010736555 2.143888e-04 3.7
     [21,] 0.0226242424 0.010333795 -3.378657e-09 8.5
     [22,] 0.0226242424 0.010342061 4.430797e-08 7.4
     [23,] 0.0237302174 0.010162521 -1.077327e-06 6.0
     [24,] 0.0447535254 0.008396264 1.222242e-05 4.9
     [25,] 0.0818184250 0.006860077 -1.786334e-09 8.7
     [26,] 0.0828003091 0.006812347 4.179976e-09 8.4
     [27,] 0.0908216581 0.006552698 -7.160336e-09 8.1
     [28,] 0.1022947605 0.006235631 -6.638898e-09 8.2
     [29,] 0.1274058577 0.005623691 6.605415e-09 8.2
     [30,] 0.1352296342 0.005400734 -2.347626e-05 4.6
     [31,] 0.1377322800 0.005330921 2.992858e-04 3.5
     [32,] 0.1381129175 0.005351387 -2.053358e-06 5.7
     [33,] 0.1411006360 0.005273058 4.315938e-05 4.4
     [34,] 0.1429052994 0.005236800 3.481808e-04 3.5
     [35,] 0.1456245572 0.005269240 -1.945018e-09 8.7
     [36,] 0.1546068729 0.005068069 -4.599247e-07 6.3
     [37,] 0.1546068729 0.005073662 2.723010e-07 6.6
     [38,] 0.1697410365 0.004841818 5.314010e-09 8.3
     [39,] 0.1773275763 0.004659561 5.534044e-05 4.3
     [40,] 0.1900940174 0.004503736 1.785654e-06 5.7
     [41,] 0.1901476415 0.004534687 5.662356e-09 8.2
     [42,] 0.2147809514 0.004198693 -1.918794e-08 7.7
     [43,] 0.2328051066 0.003953993 -9.175810e-08 7.0
     [44,] 0.2520765119 0.003749038 2.145699e-07 6.7
     [45,] 0.2530829140 0.003752594 3.654361e-09 8.4
     [46,] 0.2539220587 0.003712373 -5.090149e-06 5.3
     [47,] 0.2542892786 0.003743439 -1.056649e-09 9.0
     [48,] 0.2600174995 0.003661796 -2.904302e-07 6.5
     [49,] 0.2703239068 0.003519992 -1.561648e-04 3.8
     [50,] 0.2716993561 0.003550681 5.130930e-09 8.3
     [51,] 0.2755161961 0.003468040 -4.351715e-04 3.4
     [52,] 0.2901885431 0.003366135 -5.644431e-08 7.2
     [53,] 0.2905790220 0.003344235 1.026676e-07 7.0
     [54,] 0.2905790220 0.003367649 2.260616e-08 7.6
     [55,] 0.2918501987 0.003335528 -1.273388e-06 5.9
     [56,] 0.2965211365 0.003303089 -4.489274e-07 6.3
     [57,] 0.2980341749 0.003304623 8.424704e-09 8.1
     [58,] 0.3005567833 0.003239225 4.660033e-05 4.3
     [59,] 0.3223198463 0.003061345 -1.151308e-05 4.9
     [60,] 0.3250712721 0.003023433 -2.470887e-04 3.6
     [61,] 0.3250712721 0.003041464 -7.713766e-06 5.1
     [62,] 0.3318884124 0.003008371 -4.960989e-09 8.3
     [63,] 0.3853894768 0.002609817 -1.774302e-09 8.8
     [64,] 0.4253339570 0.002329958 1.828230e-08 7.7
     [65,] 0.4395037092 0.002224527 -4.535852e-06 5.3
     [66,] 0.4508046241 0.002167703 -4.594550e-08 7.3
     > plot(rL[,1:2], type = "b", main = "inaccurate pchisq/qchisq pairs")
     >
     > plot(rL[,1:2], type = "b", log = "x", ylim = range(0, rL[,"df"]),
     + xaxt = "n",
     + main = "inaccurate pchisq/qchisq pairs"); abline(h = 0, lty=2)
     > ## aha -- a perfect line !!
     > lines(res[,1:2], col = adjustcolor(1, 0.5))
     > eaxis(1); axis(1, at = 1/2)
     >
     > d <- as.data.frame(res)
     > plot (df ~ log(p), data = d, type = "b", cex=1/4, col="gray")
     > points(df ~ log(p), data = as.data.frame(rL), col=2, cex = 1/2)
     >
     > summary(fm <- lm (df ~ log(p), data = d, weights = -log(abs(rE))))
    
     Call:
     lm(formula = df ~ log(p), data = d, weights = -log(abs(rE)))
    
     Weighted Residuals:
     Min 1Q Median 3Q Max
     -6.912e-04 -1.441e-04 -1.969e-05 7.710e-05 1.082e-03
    
     Coefficients:
     Estimate Std. Error t value Pr(>|t|)
     (Intercept) 6.184e-06 1.133e-05 0.546 0.587
     log(p) -2.725e-03 3.660e-06 -744.530 <2e-16 ***
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     Residual standard error: 0.0002557 on 90 degrees of freedom
     Multiple R-squared: 0.9998, Adjusted R-squared: 0.9998
     F-statistic: 5.543e+05 on 1 and 90 DF, p-value: < 2.2e-16
    
     > ## R^2 = 0.9998
     >
     > p0 <- 2^seq(-50,-1, by=1/8)
     > dN <- data.frame(p = p0,
     + df = predict(fm, newdata = data.frame(p = p0)))
     > rE <- with(dN, 1- pchisq(qchisq(p, df),df)/p)
     > dN <- cbind(dN, rE = rE, nDig = round(-log10(abs(rE)), 1))
     > print(dN, digits=10)
     p df rE nDig
     1 8.881784197e-16 0.094448581017 -7.463565388e-08 7.1
     2 9.685654347e-16 0.094212475025 1.225666370e-06 5.9
     3 1.056228096e-15 0.093976369034 4.405535633e-07 6.4
     4 1.151824906e-15 0.093740263042 -3.124345540e-07 6.5
     5 1.256073967e-15 0.093504157050 1.062787617e-06 6.0
     6 1.369758374e-15 0.093268051059 3.686723420e-07 6.4
     7 1.493732098e-15 0.093031945067 -2.933092969e-07 6.5
     8 1.628926404e-15 0.092795839075 1.156863211e-06 5.9
     9 1.776356839e-15 0.092559733084 -1.520860305e-06 5.8
     10 1.937130869e-15 0.092323627092 -1.717405684e-08 7.8
     11 2.112456192e-15 0.092087521100 1.507978934e-06 5.8
     12 2.303649813e-15 0.091851415109 -1.062629703e-06 6.0
     13 2.512147934e-15 0.091615309117 5.160568806e-07 6.3
     14 2.739516748e-15 0.091379203125 6.832590693e-08 7.2
     15 2.987464197e-15 0.091143097134 -3.472479155e-07 6.5
     16 3.257852808e-15 0.090906991142 1.306469104e-06 5.9
     17 3.552713679e-15 0.090670885150 -1.081912813e-06 6.0
     18 3.874261739e-15 0.090434779159 6.253705814e-07 6.2
     19 4.224912384e-15 0.090198673167 3.330732011e-07 6.5
     20 4.607299625e-15 0.089962567175 7.294750648e-08 7.1
     21 5.024295868e-15 0.089726461183 -1.550036943e-07 6.8
     22 5.479033495e-15 0.089490355192 -3.507775794e-07 6.5
     23 5.974928394e-15 0.089254249200 1.485185219e-06 5.8
     24 6.515705616e-15 0.089018143208 -6.457821511e-07 6.2
     25 7.105427358e-15 0.088782037217 1.243792663e-06 5.9
     26 7.748523477e-15 0.088545931225 -8.120437036e-07 6.1
     27 8.449824769e-15 0.088309825233 1.131156029e-06 5.9
     28 9.214599250e-15 0.088073719242 -8.495397843e-07 6.1
     29 1.004859174e-14 0.087837613250 1.147297759e-06 5.9
     30 1.095806699e-14 0.087601507258 -7.582479513e-07 6.1
     31 1.194985679e-14 0.087365401267 1.292240281e-06 5.9
     32 1.303141123e-14 0.087129295275 -5.381458088e-07 6.3
     33 1.421085472e-14 0.086893189283 -3.797838852e-07 6.4
     34 1.549704695e-14 0.086657083292 -1.892109929e-07 6.7
     35 1.689964954e-14 0.086420977300 3.357565381e-08 7.5
     36 1.842919850e-14 0.086184871308 2.885788447e-07 6.5
     37 2.009718347e-14 0.085948765317 -1.348411111e-06 5.9
     38 2.191613398e-14 0.085712659325 8.952460000e-07 6.0
     39 2.389971358e-14 0.085476553333 -6.665524503e-07 6.2
     40 2.606282246e-14 0.085240447341 -2.772836227e-07 6.6
     41 2.842170943e-14 0.085004341350 1.442152372e-07 6.8
     42 3.099409391e-14 0.084768235358 -1.299326189e-06 5.9
     43 3.379929908e-14 0.084532129366 1.083914191e-06 6.0
     44 3.685839700e-14 0.084296023375 -2.844134916e-07 6.5
     45 4.019436694e-14 0.084059917383 2.714025346e-07 6.6
     46 4.383226796e-14 0.083823811391 8.594620383e-07 6.1
     47 4.779942715e-14 0.083587705400 -3.905788137e-07 6.4
     48 5.212564492e-14 0.083351599408 2.673437002e-07 6.6
     49 5.684341886e-14 0.083115493416 9.575175537e-07 6.0
     50 6.198818782e-14 0.082879387425 -1.742195652e-07 6.8
     51 6.759859815e-14 0.082643281433 -1.262888024e-06 5.9
     52 7.371679400e-14 0.082407175441 1.378141781e-06 5.9
     53 8.038873388e-14 0.082171069450 3.647252201e-07 6.4
     54 8.766453592e-14 0.081934963458 -6.056185433e-07 6.2
     55 9.559885430e-14 0.081698857466 2.942142103e-07 6.5
     56 1.042512898e-13 0.081462751475 1.226316300e-06 5.9
     57 1.136868377e-13 0.081226645483 3.743139556e-07 6.4
     58 1.239763756e-13 0.080990539491 -4.346107259e-07 6.4
     59 1.351971963e-13 0.080754433499 -1.200456962e-06 5.9
     60 1.474335880e-13 0.080518327508 -1.231656597e-07 6.9
     61 1.607774678e-13 0.080282221516 9.864073824e-07 6.0
     62 1.753290718e-13 0.080046115524 3.389273247e-07 6.5
     63 1.911977086e-13 0.079810009533 -2.654694915e-07 6.6
     64 2.085025797e-13 0.079573903541 -8.267823173e-07 6.1
     65 2.273736754e-13 0.079337797549 4.280209399e-07 6.4
     66 2.479527513e-13 0.079101691558 -5.255530455e-08 7.3
     67 2.703943926e-13 0.078865585566 1.272196181e-06 5.9
     68 2.948671760e-13 0.078629479574 8.723627990e-07 6.1
     69 3.215549355e-13 0.078393373583 5.156188723e-07 6.3
     70 3.506581437e-13 0.078157267591 2.019651378e-07 6.7
     71 3.823954172e-13 0.077921161599 -6.859768598e-08 7.2
     72 4.170051594e-13 0.077685055608 -2.960688799e-07 6.5
     73 4.547473509e-13 0.077448949616 -4.804477280e-07 6.3
     74 4.959055026e-13 0.077212843624 -6.217335187e-07 6.2
     75 5.407887852e-13 0.076976737633 -7.199255436e-07 6.1
     76 5.897343520e-13 0.076740631641 -7.750230990e-07 6.1
     77 6.431098711e-13 0.076504525649 -7.870254888e-07 6.1
     78 7.013162874e-13 0.076268419657 -7.559320068e-07 6.1
     79 7.647908344e-13 0.076032313666 -6.817419649e-07 6.2
     80 8.340103188e-13 0.075796207674 -5.644546737e-07 6.2
     81 9.094947018e-13 0.075560101682 -4.040694535e-07 6.4
     82 9.918110051e-13 0.075323995691 -2.005856121e-07 6.7
     83 1.081577570e-12 0.075087889699 4.599752257e-08 7.3
     84 1.179468704e-12 0.074851783707 3.356806249e-07 6.5
     85 1.286219742e-12 0.074615677716 6.684643596e-07 6.2
     86 1.402632575e-12 0.074379571724 1.044349404e-06 6.0
     87 1.529581669e-12 0.074143465732 -1.904324292e-07 6.7
     88 1.668020638e-12 0.073907359741 2.770713794e-07 6.6
     89 1.818989404e-12 0.073671253749 -8.551932531e-07 6.1
     90 1.983622010e-12 0.073435147757 -2.960718908e-07 6.5
     91 2.163155141e-12 0.073199041766 3.061522486e-07 6.5
     92 2.358937408e-12 0.072962935774 9.514798076e-07 6.0
     93 2.572439484e-12 0.072726829782 1.875341071e-08 7.7
     94 2.805265149e-12 0.072490723791 -8.599822618e-07 6.1
     95 3.059163338e-12 0.072254617799 -7.452604778e-08 7.1
     96 3.336041275e-12 0.072018511807 7.540343130e-07 6.1
     97 3.637978807e-12 0.071782405815 2.630667650e-08 7.6
     98 3.967244020e-12 0.071546299824 -6.474386178e-07 6.2
     99 4.326310282e-12 0.071310193832 3.212443331e-07 6.5
     100 4.717874816e-12 0.071074087840 -2.500146477e-07 6.6
     101 5.144878969e-12 0.070837981849 -7.672983209e-07 6.1
     102 5.610530299e-12 0.070601875857 3.415009425e-07 6.5
     103 6.118326675e-12 0.070365769865 -7.330663121e-08 7.1
     104 6.672082550e-12 0.070129663874 -4.341459592e-07 6.4
     105 7.275957614e-12 0.069893557882 8.147632093e-07 6.1
     106 7.934488041e-12 0.069657451890 -9.939312102e-07 6.0
     107 8.652620563e-12 0.069421345899 3.519772401e-07 6.5
     108 9.435749632e-12 0.069185239907 2.015231014e-07 6.7
     109 1.028975794e-11 0.068949133915 1.050243966e-07 7.0
     110 1.122106060e-11 0.068713027924 6.247828421e-08 7.2
     111 1.223665335e-11 0.068476921932 7.388190371e-08 7.1
     112 1.334416510e-11 0.068240815940 1.392323935e-07 6.9
     113 1.455191523e-11 0.068004709949 2.585268818e-07 6.6
     114 1.586897608e-11 0.067768603957 -1.074910069e-06 6.0
     115 1.730524113e-11 0.067532497965 6.589363856e-07 6.2
     116 1.887149926e-11 0.067296391974 -5.557286085e-07 6.3
     117 2.057951587e-11 0.067060285982 -2.152397722e-07 6.7
     118 2.244212120e-11 0.066824179990 1.791772732e-07 6.7
     119 2.447330670e-11 0.066588073998 -8.518538748e-07 6.1
     120 2.668833020e-11 0.066351968007 -3.441463585e-07 6.5
     121 2.910383046e-11 0.066115862015 2.174792778e-07 6.7
     122 3.173795216e-11 0.065879756023 -6.299664652e-07 6.2
     123 3.461048225e-11 0.065643650032 4.492504890e-08 7.3
     124 3.774299853e-11 0.065407544040 7.737245498e-07 6.1
     125 4.115903175e-11 0.065171438048 1.098156234e-07 7.0
     126 4.488424239e-11 0.064935332057 -4.892612773e-07 6.3
     127 4.894661340e-11 0.064699226065 -1.023513942e-06 6.0
     128 5.337666040e-11 0.064463120073 -6.281881859e-08 7.2
     129 5.820766091e-11 0.064227014082 9.517654395e-07 6.0
     130 6.347590433e-11 0.063990908090 6.009614978e-07 6.2
     131 6.922096451e-11 0.063754802098 3.149542980e-07 6.5
     132 7.548599706e-11 0.063518696107 9.373602483e-08 7.0
     133 8.231806350e-11 0.063282590115 -6.270114161e-08 7.2
     134 8.976848478e-11 0.063046484123 -1.543650285e-07 6.8
     135 9.789322680e-11 0.062810378132 -1.812634687e-07 6.7
     136 1.067533208e-10 0.062574272140 -1.434043004e-07 6.8
     137 1.164153218e-10 0.062338166148 -4.079537130e-08 7.4
     138 1.269518087e-10 0.062102060156 1.265554705e-07 6.9
     139 1.384419290e-10 0.061865954165 3.586403706e-07 6.4
     140 1.509719941e-10 0.061629848173 6.554514619e-07 6.2
     141 1.646361270e-10 0.061393742181 1.016980880e-06 6.0
     142 1.795369696e-10 0.061157636190 8.972756516e-08 7.0
     143 1.957864536e-10 0.060921530198 -7.618360451e-07 6.1
     144 2.135066416e-10 0.060685424206 -1.952735760e-07 6.7
     145 2.328306437e-10 0.060449318215 4.359719799e-07 6.4
     146 2.539036173e-10 0.060213212223 -1.996387393e-07 6.7
     147 2.768838580e-10 0.059977106231 -7.596117781e-07 6.1
     148 3.019439882e-10 0.059741000240 7.654519618e-08 7.1
     149 3.292722540e-10 0.059504894248 9.773499157e-07 6.0
     150 3.590739391e-10 0.059268788256 -6.763557745e-07 6.2
     151 3.915729072e-10 0.059032682265 3.646181882e-07 6.4
     152 4.270132832e-10 0.058796576273 1.716164724e-07 6.8
     153 4.656612873e-10 0.058560470281 5.417191207e-08 7.3
     154 5.078072346e-10 0.058324364290 1.227028412e-08 7.9
     155 5.537677160e-10 0.058088258298 4.589736990e-08 7.3
     156 6.038879765e-10 0.057852152306 1.550389357e-07 6.8
     157 6.585445080e-10 0.057616046314 3.396807524e-07 6.5
     158 7.181478783e-10 0.057379940323 -6.657672476e-07 6.2
     159 7.831458144e-10 0.057143834331 9.354081640e-07 6.0
     160 8.540265665e-10 0.056907728339 9.185895744e-08 7.0
     161 9.313225746e-10 0.056671622348 -6.652309521e-07 6.2
     162 1.015614469e-09 0.056435516356 -9.234395604e-08 7.0
     163 1.107535432e-09 0.056199410364 5.559531132e-07 6.3
     164 1.207775953e-09 0.055963304373 4.705635348e-08 7.3
     165 1.317089016e-09 0.055727198381 -3.754633768e-07 6.4
     166 1.436295757e-09 0.055491092389 -7.116280665e-07 6.1
     167 1.566291629e-09 0.055254986398 2.545663202e-07 6.6
     168 1.708053133e-09 0.055018880406 8.553017516e-08 7.1
     169 1.862645149e-09 0.054782774414 2.785924047e-09 8.6
     170 2.031228938e-09 0.054546668423 6.311568379e-09 8.2
     171 2.215070864e-09 0.054310562431 9.608511653e-08 7.0
     172 2.415551906e-09 0.054074456439 2.720845634e-07 6.6
     173 2.634178032e-09 0.053838350448 -6.486866118e-07 6.2
     174 2.872591513e-09 0.053602244456 -2.948026030e-07 6.5
     175 3.132583258e-09 0.053366138464 1.452384917e-07 6.8
     176 3.416106266e-09 0.053130032472 6.714146586e-07 6.2
     177 3.725290298e-09 0.052893926481 1.227572537e-07 6.9
     178 4.062457877e-09 0.052657820489 -3.287850758e-07 6.5
     179 4.430141728e-09 0.052421714497 4.666330337e-07 6.3
     180 4.831103812e-09 0.052185608506 2.036972987e-07 6.7
     181 5.268356064e-09 0.051949502514 3.778649948e-08 7.4
     182 5.745183026e-09 0.051713396522 -3.113050129e-08 7.5
     183 6.265166516e-09 0.051477290531 -3.084833278e-09 8.5
     184 6.832212532e-09 0.051241184539 1.218923642e-07 6.9
     185 7.450580597e-09 0.051005078547 3.437699647e-07 6.5
     186 8.124915754e-09 0.050768972556 -4.487175920e-07 6.3
     187 8.860283457e-09 0.050532866564 -2.762538487e-08 7.6
     188 9.662207623e-09 0.050296760572 4.902705625e-07 6.3
     189 1.053671213e-08 0.050060654581 1.026316276e-08 8.0
     190 1.149036605e-08 0.049824548589 -3.619620201e-07 6.4
     191 1.253033303e-08 0.049588442597 -6.264466410e-07 6.2
     192 1.366442506e-08 0.049352336606 -7.832323377e-07 6.1
     193 1.490116119e-08 0.049116230614 2.401634375e-07 6.6
     194 1.624983151e-08 0.048880124622 -7.738734984e-07 6.1
     195 1.772056691e-08 0.048644018630 -6.078122046e-07 6.2
     196 1.932441525e-08 0.048407912639 -3.342184687e-07 6.5
     197 2.107342426e-08 0.048171806647 4.686610289e-08 7.3
     198 2.298073210e-08 0.047935700655 5.353999194e-07 6.3
     199 2.506066606e-08 0.047699594664 9.197667672e-08 7.0
     200 2.732885013e-08 0.047463488672 -2.330247564e-07 6.6
     201 2.980232239e-08 0.047227382680 5.886113930e-07 6.2
     202 3.249966302e-08 0.046991276689 -5.279761615e-07 6.3
     203 3.544113383e-08 0.046755170697 5.191576196e-07 6.3
     204 3.864883049e-08 0.046519064705 -3.498820569e-07 6.5
     205 4.214684851e-08 0.046282958714 -8.357652814e-08 7.1
     206 4.596146421e-08 0.046046852722 3.008301066e-07 6.5
     207 5.012133212e-08 0.045810746730 -1.917574097e-07 6.7
     208 5.465770025e-08 0.045574640739 -5.552768079e-07 6.3
     209 5.960464478e-08 0.045338534747 1.941360427e-07 6.7
     210 6.499932603e-08 0.045102428755 8.300436694e-08 7.1
     211 7.088226765e-08 0.044866322764 1.007455417e-07 7.0
     212 7.729766099e-08 0.044630216772 -7.200064360e-07 6.1
     213 8.429369702e-08 0.044394110780 -4.391848569e-07 6.4
     214 9.192292842e-08 0.044158004788 -2.969477242e-08 7.5
     215 1.002426642e-07 0.043921898797 5.083973372e-07 6.3
     216 1.093154005e-07 0.043685792805 2.298753619e-07 6.6
     217 1.192092896e-07 0.043449686813 9.089335073e-08 7.0
     218 1.299986521e-07 0.043213580822 9.137020085e-08 7.0
     219 1.417645353e-07 0.042977474830 2.312248708e-07 6.6
     220 1.545953220e-07 0.042741368838 -4.125440396e-07 6.4
     221 1.685873940e-07 0.042505262847 1.135778693e-08 7.9
     222 1.838458568e-07 0.042269156855 5.743894688e-07 6.2
     223 2.004853285e-07 0.042033050863 3.701760187e-07 6.4
     224 2.186308010e-07 0.041796944872 3.160064802e-07 6.5
     225 2.384185791e-07 0.041560838880 4.117840361e-07 6.4
     226 2.599973041e-07 0.041324732888 6.574119313e-07 6.2
     227 2.835290706e-07 0.041088626897 1.687303717e-07 6.8
     228 3.091906439e-07 0.040852520905 -1.591887169e-07 6.8
     229 3.371747881e-07 0.040616414913 5.464674205e-07 6.3
     230 3.676917137e-07 0.040380308922 -3.331954996e-07 6.5
     231 4.009706570e-07 0.040144202930 -1.795109488e-07 6.7
     232 4.372616020e-07 0.039908096938 1.344805564e-07 6.9
     233 4.768371582e-07 0.039671990946 -2.420257981e-07 6.6
     234 5.199946082e-07 0.039435884955 -4.473476847e-07 6.3
     235 5.670581412e-07 0.039199778963 -4.816173824e-07 6.3
     236 6.183812879e-07 0.038963672971 -3.449670454e-07 6.5
     237 6.743495762e-07 0.038727566980 -3.752869437e-08 7.4
     238 7.353834273e-07 0.038491460988 4.405657880e-07 6.4
     239 8.019413140e-07 0.038255354996 -5.454495733e-07 6.3
     240 8.745232040e-07 0.038019249005 2.846233194e-07 6.5
     241 9.536743164e-07 0.037783143013 -3.274450626e-07 6.5
     242 1.039989216e-06 0.037547037021 5.340015219e-08 7.3
     243 1.134116282e-06 0.037310931030 -1.797782818e-07 6.7
     244 1.236762576e-06 0.037074825038 5.685058231e-07 6.2
     245 1.348699152e-06 0.036838719046 -7.036316796e-08 7.2
     246 1.470766855e-06 0.036602613055 -5.064464588e-07 6.3
     247 1.603882628e-06 0.036366507063 3.281364691e-08 7.5
     248 1.749046408e-06 0.036130401071 -3.854574882e-09 8.4
     249 1.907348633e-06 0.035894295080 1.616862507e-07 6.8
     250 2.079978433e-06 0.035658189088 -2.268208801e-07 6.6
     251 2.268232565e-06 0.035422083096 -4.023808344e-07 6.4
     252 2.473525152e-06 0.035185977105 -3.652116916e-07 6.4
     253 2.697398305e-06 0.034949871113 -1.155312188e-07 6.9
     254 2.941533709e-06 0.034713765121 3.464431350e-07 6.5
     255 3.207765256e-06 0.034477659129 -4.359180803e-07 6.4
     256 3.498092816e-06 0.034241553138 4.610668849e-07 6.3
     257 3.814697266e-06 0.034005447146 1.355166542e-07 6.9
     258 4.159956866e-06 0.033769341154 4.357489369e-08 7.4
     259 4.536465130e-06 0.033533235163 1.849733936e-07 6.7
     260 4.947050303e-06 0.033297129171 -1.408841688e-07 6.9
     261 5.394796609e-06 0.033061023179 -2.228060141e-07 6.7
     262 5.883067419e-06 0.032824917188 -6.108725370e-08 7.2
     263 6.415530512e-06 0.032588811196 3.439775413e-07 6.5
     264 6.996185632e-06 0.032352705204 3.140542086e-07 6.5
     265 7.629394531e-06 0.032116599213 -1.344452274e-07 6.9
     266 8.319913732e-06 0.031880493221 -3.182555455e-07 6.5
     267 9.072930260e-06 0.031644387229 4.235536309e-07 6.4
     268 9.894100606e-06 0.031408281238 1.067806606e-07 7.0
     269 1.078959322e-05 0.031172175246 6.478579662e-08 7.2
     270 1.176613484e-05 0.030936069254 2.971861416e-07 6.5
     271 1.283106102e-05 0.030699963263 -4.743406139e-07 6.3
     272 1.399237126e-05 0.030463857271 3.167996901e-07 6.5
     273 1.525878906e-05 0.030227751279 1.255002516e-07 6.9
     274 1.663982746e-05 0.029991645287 2.292534416e-07 6.6
     275 1.814586052e-05 0.029755539296 1.095587332e-08 8.0
     276 1.978820121e-05 0.029519433304 9.794742695e-08 7.0
     277 2.157918644e-05 0.029283327312 -1.157428633e-07 6.9
     278 2.353226967e-05 0.029047221321 -1.397550675e-08 7.9
     279 2.566212205e-05 0.028811115329 4.027337354e-07 6.4
     280 2.798474253e-05 0.028575009337 -4.365435546e-08 7.4
     281 3.051757812e-05 0.028338903346 -1.538802443e-07 6.8
     282 3.327965493e-05 0.028102797354 7.146577696e-08 7.1
     283 3.629172104e-05 0.027866691362 5.978360873e-08 7.2
     284 3.957640242e-05 0.027630585371 3.936098777e-07 6.4
     285 4.315837288e-05 0.027394479379 -4.939450760e-08 7.3
     286 4.706453935e-05 0.027158373387 -1.258946654e-07 6.9
     287 5.132424410e-05 0.026922267396 -3.862692297e-07 6.4
     288 5.596948506e-05 0.026686161404 -2.704251889e-07 6.6
     289 6.103515625e-05 0.026450055412 2.208878911e-07 6.7
     290 6.655930986e-05 0.026213949421 2.101241836e-08 7.7
     291 7.258344208e-05 0.025977843429 2.173189355e-07 6.7
     292 7.915280485e-05 0.025741737437 -2.345979131e-07 6.6
     293 8.631674575e-05 0.025505631445 -2.697597914e-07 6.6
     294 9.412907870e-05 0.025269525454 1.109045821e-07 7.0
     295 1.026484882e-04 0.025033419462 -1.036342423e-07 7.0
     296 1.119389701e-04 0.024797313470 1.180052477e-07 6.9
     297 1.220703125e-04 0.024561207479 -2.129809857e-07 6.7
     298 1.331186197e-04 0.024325101487 -8.761589876e-08 7.1
     299 1.451668842e-04 0.024088995495 1.024239626e-08 8.0
     300 1.583056097e-04 0.023852889504 9.613828733e-08 7.0
     301 1.726334915e-04 0.023616783512 1.855632931e-07 6.7
     302 1.882581574e-04 0.023380677520 2.939531960e-07 6.5
     303 2.052969764e-04 0.023144571529 -2.374739316e-08 7.6
     304 2.238779402e-04 0.022908465537 1.742085209e-07 6.8
     305 2.441406250e-04 0.022672359545 -1.221569534e-08 7.9
     306 2.662372394e-04 0.022436253554 -1.074149700e-07 7.0
     307 2.903337683e-04 0.022200147562 -9.628875608e-08 7.0
     308 3.166112194e-04 0.021964041570 3.618904520e-08 7.4
     309 3.452669830e-04 0.021727935579 3.049675551e-07 6.5
     310 3.765163148e-04 0.021491829587 3.034669259e-07 6.5
     311 4.105939528e-04 0.021255723595 6.315431433e-08 7.2
     312 4.477558805e-04 0.021019617603 2.569936886e-08 7.6
     313 4.882812500e-04 0.020783511612 -1.990711820e-07 6.7
     314 5.324744788e-04 0.020547405620 -1.808471892e-07 6.7
     315 5.806675366e-04 0.020311299628 9.470999129e-08 7.0
     316 6.332224388e-04 0.020075193637 2.537703486e-07 6.6
     317 6.905339660e-04 0.019839087645 -5.540736026e-08 7.3
     318 7.530326296e-04 0.019602981653 -3.165362417e-08 7.5
     319 8.211879055e-04 0.019366875662 -3.256026071e-08 7.5
     320 8.955117609e-04 0.019130769670 -2.782856834e-08 7.6
     321 9.765625000e-04 0.018894663678 1.266636318e-08 7.9
     322 1.064948958e-03 0.018658557687 -2.358667284e-07 6.6
     323 1.161335073e-03 0.018422451695 -2.867862725e-08 7.5
     324 1.266444878e-03 0.018186345703 -4.015371813e-08 7.4
     325 1.381067932e-03 0.017950239712 -2.243580433e-07 6.6
     326 1.506065259e-03 0.017714133720 1.295179664e-07 6.9
     327 1.642375811e-03 0.017478027728 5.258487723e-08 7.3
     328 1.791023522e-03 0.017241921737 -7.191778106e-08 7.1
     329 1.953125000e-03 0.017005815745 1.168092258e-07 6.9
     330 2.129897915e-03 0.016769709753 2.579736935e-08 7.6
     331 2.322670146e-03 0.016533603761 2.070165972e-08 7.7
     332 2.532889755e-03 0.016297497770 1.453019420e-07 6.8
     333 2.762135864e-03 0.016061391778 -1.448045621e-07 6.8
     334 3.012130518e-03 0.015825285786 9.153762492e-08 7.0
     335 3.284751622e-03 0.015589179795 3.166067630e-08 7.5
     336 3.582047044e-03 0.015353073803 2.793761855e-08 7.6
     337 3.906250000e-03 0.015116967811 -1.335779407e-07 6.9
     338 4.259795831e-03 0.014880861820 -1.127023992e-07 6.9
     339 4.645340293e-03 0.014644755828 1.472359789e-07 6.8
     340 5.065779510e-03 0.014408649836 1.913901931e-07 6.7
     341 5.524271728e-03 0.014172543845 1.078438131e-07 7.0
     342 6.024261037e-03 0.013936437853 -1.618990075e-08 7.8
     343 6.569503244e-03 0.013700331861 -9.444800253e-08 7.0
     344 7.164094088e-03 0.013464225870 -4.169422585e-08 7.4
     345 7.812500000e-03 0.013228119878 -1.890127832e-09 8.7
     346 8.519591661e-03 0.012992013886 -9.842786630e-08 7.0
     347 9.290680586e-03 0.012755907895 1.241692305e-10 9.9
     348 1.013155902e-02 0.012519801903 -3.338975185e-08 7.5
     349 1.104854346e-02 0.012283695911 1.347443089e-07 6.9
     350 1.204852207e-02 0.012047589919 1.089463753e-08 8.0
     351 1.313900649e-02 0.011811483928 -5.345811105e-08 7.3
     352 1.432818818e-02 0.011575377936 7.911680988e-08 7.1
     353 1.562500000e-02 0.011339271944 -1.082518675e-08 8.0
     354 1.703918332e-02 0.011103165953 3.672092086e-08 7.4
     355 1.858136117e-02 0.010867059961 3.491128742e-08 7.5
     356 2.026311804e-02 0.010630953969 4.962150579e-09 8.3
     357 2.209708691e-02 0.010394847978 -1.451724585e-08 7.8
     358 2.409704415e-02 0.010158741986 3.193182674e-08 7.5
     359 2.627801298e-02 0.009922635994 6.318548917e-08 7.2
     360 2.865637635e-02 0.009686530003 3.477948718e-08 7.5
     361 3.125000000e-02 0.009450424011 8.024744635e-08 7.1
     362 3.407836665e-02 0.009214318019 7.269387325e-08 7.1
     363 3.716272234e-02 0.008978212028 -6.037044309e-08 7.2
     364 4.052623608e-02 0.008742106036 4.128745124e-08 7.4
     365 4.419417382e-02 0.008506000044 -1.954348461e-09 8.7
     366 4.819408829e-02 0.008269894053 6.703489763e-09 8.2
     367 5.255602595e-02 0.008033788061 6.777183081e-08 7.2
     368 5.731275270e-02 0.007797682069 3.394610371e-08 7.5
     369 6.250000000e-02 0.007561576077 6.149548670e-08 7.2
     370 6.815673329e-02 0.007325470086 -2.079005479e-08 7.7
     371 7.432544469e-02 0.007089364094 2.085197826e-08 7.7
     372 8.105247217e-02 0.006853258102 -3.090769152e-08 7.5
     373 8.838834765e-02 0.006617152111 -5.460620156e-08 7.3
     374 9.638817659e-02 0.006381046119 4.531310338e-08 7.3
     375 1.051120519e-01 0.006144940127 1.717334919e-08 7.8
     376 1.146255054e-01 0.005908834136 1.463540500e-08 7.8
     377 1.250000000e-01 0.005672728144 3.113683844e-08 7.5
     378 1.363134666e-01 0.005436622152 1.390943816e-08 7.9
     379 1.486508894e-01 0.005200516161 -3.333184884e-08 7.5
     380 1.621049443e-01 0.004964410169 -1.011861062e-08 8.0
     381 1.767766953e-01 0.004728304177 -2.307868430e-08 7.6
     382 1.927763532e-01 0.004492198186 -3.451609221e-09 8.5
     383 2.102241038e-01 0.004256092194 -1.417570927e-08 7.8
     384 2.292510108e-01 0.004019986202 -5.527056146e-09 8.3
     385 2.500000000e-01 0.003783880211 -9.496292419e-09 8.0
     386 2.726269332e-01 0.003547774219 -4.487596073e-09 8.3
     387 2.973017788e-01 0.003311668227 -4.862096725e-09 8.3
     388 3.242098887e-01 0.003075562235 -7.645436728e-09 8.1
     389 3.535533906e-01 0.002839456244 3.773244051e-09 8.4
     390 3.855527064e-01 0.002603350252 5.646113910e-09 8.2
     391 4.204482076e-01 0.002367244260 -3.754701661e-09 8.4
     392 4.585020216e-01 0.002131138269 -2.839019464e-09 8.5
     393 5.000000000e-01 0.001895032277 -5.991975804e-10 9.2
     >
     > ## } ## only when we find inaccurate regions
     > showProc.time()
     Time elapsed: 0.1 0.03 0.17
     >
     >
     > ## Oops: another qgamma() / qchisq() problem: mostly NaN's == all solved now
     > curve(qgamma(x, 20), 1e-16, 1e-10, log='x')
     > curve(qgamma(x, 20), 1e-300, .99 , log='xy') # and add the critical region from above:
     > abline(v=c(1e-16, 1e-10), col="light blue")
     > curve(qgamma(x, 20), 1e-26, 1e-07, log='x')
     > ##-> now using log=TRUE in same region:
     > curve(qgamma(x, 20, log=TRUE), -38, -16)## no problem!!
     > curve(qgamma(exp(x), 20), add=TRUE, col="green3", n=2001)
     > ## had problem here, but no longer !
     >
     > ##--> Further fix for qgamma: when 'x' is very small: use "log=TRUE of log(x)"!
     >
     > ## had bug (gave NaN), but no longer:
     > (q_12 <- qgamma(1e-12, 20))
     [1] 2.330042
     > all.equal(1e-12, pgamma(q_12, 20), tol=0)# show rel.err (Lnx 64-bit: 4.04e-16)
     [1] "Mean relative difference: 8.279884e-15"
     > stopifnot(
     + all.equal(1e-12, pgamma(q_12, 20), tolerance = 1e-14)
     + )
     >
     >
     > ## --- Nice graphic : --- but amazingly *S..L..O..W*
     >
     > p.qgammaSml <- function(from= 1e-110, to = 1e-5, ylim = c(0.4, 1000),
     + n = 201, k.lab = 3,
     + a1 = c(10, seq(10.1,20, by=.2), 21:105),
     + a2 = seq(110,330, by=10),
     + a3 = seq(350,1600, by=50))
     + {
     + ## Purpose: nice qgamma() lines ``for small x'' aka p
     + ## ----------------------------------------------------------------------
     + ## Author: Martin Maechler, Date: 22 Mar 2004, 14:23
     + x <- exp(seq(log(from), log(to), length = n))
     +
     + op <- par(las=1, lab = c(10,10, 7), xaxs = "i", mex = 0.8)
     + on.exit(par(op))
     + plot(x, qgamma(x, a1[1]), log="xy", ylim=ylim, type='l', xaxt = "n",
     + main = paste("qgamma(x, a) for very small x, a in [",
     + formatC(a1[1]),", ",formatC(max(a1,a2,a3)),"] - log-log", sep=''),
     + sub = R.version.string)
     + lab.x <- pretty(log10(c(from,to)), 20)
     + axis(1, at=10^lab.x, lab = paste("10^",formatC(lab.x),sep=''))
     + if(is.nan(qgamma(1e-12, 20)))
     + text(1e-60, 20, "all NaN", cex = 2)
     + if(!is.finite(qgamma(1e-140, 155)))
     + text(1e-240, 5, "all +Inf", cex = 2)
     +
     + lines.txt <- function(a.s, col = par("col")) {
     + col <- rep(col, length=length(a.s))
     + for(i in seq(along=a.s)) {
     + qx <- qgamma(x, (a <- a.s[i]))
     + if(i %% k.lab == 0 &&
     + any(ifi <- is.finite(qx) & qx >= ylim[1])) {
     + ik <- (i%%(2*k.lab))/k.lab # = 0 or 1
     + j <- quantile(which(ifi), c(.02,(1:3)/4+ ik/10, .98))
     + ## "segments" around the labels :
     + i0 <- 1
     + for(jj in j) {
     + ii <- i0:(jj-1)
     + i2 <- jj + -1:1
     + lines(x[ii], qx[ii], col=col[i])
     + lines(x[i2], qx[i2], col=col[i], type = 'c')
     + i0 <- jj+1
     + }
     + text(x[j], qx[j], formatC(a), col= "gray40", cex = 0.8)
     + }
     + else
     + lines(x, qx, col=col[i])
     +
     + }
     + }
     + oo <- options(warn = -1)
     + lines.txt(a1[-1])
     + lines.txt(a2, col= 2)
     + lines.txt(a3, col= rainbow(length(a3), .8, .8,
     + start = (max(a3)-min(a3))/(1+max(a3))))
     + invisible(options(oo))
     + }
     >
     > showProc.time()
     Time elapsed: 0.03 0.02 0.06
     >
     > p.qgammaSml()
     > p.qgammaSml(1e-300)
     > p.qgammaSml(1e-300,1e-50, a2= seq(100,360, by=4), a3=seq(350,1500, by=10))
     >
     > showProc.time()
     Time elapsed: 1.5 0.01 1.51
     >
     > ## The "upper" problematic corner:
     > p.qgammaSml(1e-19, 1e-3, a2=NULL,a3=NULL, ylim=c(.1,20))
     > p.qgammaSml(1e-19, 1e-3, a2=seq(1,12, by=.04), ylim=c(.1,20),a3=NULL,k.lab=10)
     > ## now shows the problem (quite well):
     > ## could it be in pgamma()'s inaccuracy, leading to qgamma() bias ?
     > aa <- c(seq(9, 22, by=0.25),seq(22.3,40,by=0.4))
     > caa <- formatC(range(aa))
     > sfsmisc::mult.fig(2)
     > curve(pgamma(x, sh=aa[1]), 0.5, 20, log = 'xy', ylim = c(1e-60, .2),
     + main = sprintf("pgamma(x, a) for a in [%s,%s]", caa[1],caa[2]))
     > for(sh in aa) curve(pgamma(x, sh), add = TRUE, col=2)
     > abline(h=c(1e-15), col="light blue", lty=2)
     >
     > curve(pgamma(x, sh=aa[1]), 0.5, 20, log = 'xy', ylim = c(1e-15, .8),
     + main = sprintf("pgamma(x, a) for a in [%s,%s]", caa[1],caa[2]))
     > for(sh in aa) curve(pgamma(x, sh), add = TRUE, col=2)
     > ## the "border curve" between "Pearson" and "Continued fraction (upper tail)"
     > ## in pgamma.c :
     > curve(pgamma(max(1,x), x), add = TRUE, col=4)
     > ## ==> pgamma() is perfect here {series expansion up to eps_C accuracy}!
     >
     > aa <- c(seq(9, 22, by=0.25),seq(22.3,40.4,by=0.4))
     > p.qgammaSml(1e-24, 1e-5, a1=aa, a2=NULL,a3=NULL, ylim=c(.8,8))
     > ## -------- save the above?
     > aa1 <- c(aa,seq(40.5,90, by=0.5))
     > p.qgammaSml(1e-60, 1e-5, a1=aa1, a2=NULL,a3=NULL, ylim=c(.9, 16))
     > aa2 <- c(aa1, seq(91,150, by= 1))
     > p.qgammaSml(1e-90, 1e-5, a1=aa2, a2=NULL,a3=NULL, ylim=c(.9, 35))
     > aa3 <- c(aa2, seq(150,250, by= 2), seq(253, 400, by=5))
     > p.qgammaSml(1e-200, 1e-5, a1=aa3, a2=NULL,a3=NULL, ylim=c(.9, 100))
     > p.qgammaSml(1e-200, 1e-5, a1=aa3, a2=NULL,a3=NULL, ylim=c(.9, 200),k.lab=9e9)
     > p.qgammaSml(1e-60, 1e-5, a1=aa3, a2=NULL,a3=NULL, ylim=c(.9, 200),k.lab=9e9)
     >
     > showProc.time()
     Time elapsed: 4.21 0.04 4.29
     >
     > ## lower a \> 10
     >
     > curve(qgamma(x, 19), 1e-14, 1e-9, log='x')
     > curve(qgamma(x, 18), 1e-14, 1e-9, log='x')
     > curve(qgamma(x, 15), 1e-11, 5e-9, log='x')
     > curve(qgamma(x, 13), 5e-10, 1e-8, log='x')
     > curve(qgamma(x, 11), 1e-8, 5e-8, log='x')
     > curve(qgamma(x, 10.5), 4.2e-8, 6e-8, log='x')
     > curve(qgamma(x, 10.3), 6e-8, 7e-8, log='x')
     > curve(qgamma(x, 10.2), 7.1e-8, 7.6e-8, log='x')
     > curve(qgamma(x, 10.15),7.7e-8, 7.9e-8, log='x')
     > curve(qgamma(x, 10.14),7.88e-8,7.92e-8, log='x',n=10001)
     >
     > ## no more problems for smaller a!! here:
     > curve(qgamma(x, 10.13), 1e-10, 5e-4, log='x',n=20001)
     > curve(qgamma(x, 10.12), 1e-10, 5e-4, log='x',n=20001)
     > curve(qgamma(x, 10.1), 1e-10, 5e-4, log='x',n=20001)
     >
     > showProc.time()
     Time elapsed: 0.79 0.01 0.85
     >
     > ##--- the "+Inf" / premature "0" case:
     > curve(qgamma(x, 155, log=TRUE), -1500, 0, log='y', n=2001,col=2)
     > curve(qgamma(x, 1e3, log=TRUE), -1500, 0, log='y', n=2001,col=2)
     > ## now works, but slowly and with kink
     > curve(qgamma (x, 1e5, log=TRUE), -3e5, 0, log='y', n=2001,col=2,lwd=3)
     > curve(qgammaAppr(x, 1e5, log=TRUE), add = TRUE, n=2001, col="blue",lwd=.4)
     > ## --- curves are almost "identical"
     > ## ===> the kink *does* come from the initial approx... hmm
     >
     > ## still "identical"
     > curve(qgamma (x, 1e4, log=TRUE), -3e4, 0, log='y', n=2001,col=2)
     > curve(qgammaAppr(x, 1e4, log=TRUE), add = TRUE, n=2001, col="tomato3")
     >
     > ## now see some difference (approx. has kink at ~ -165)
     > curve(qgamma (x, 100, log=TRUE), -200, 0, log='y', n=2001,col=2)
     > curve(qgammaAppr(x, 100, log=TRUE), add = TRUE, n=2001, col="tomato3")
     > ##
     > (kk <- 100 * 2/1.24)# 161.29
     [1] 161.2903
     > curve(qgamma (x, 100, log=TRUE), -1.1*kk, -.95*kk, log='y', n=2001,col=2)
     > curve(qgammaAppr(x, 100, log=TRUE), add = TRUE, n=2001, col="tomato3")
     > abline(v = -kk, col='blue', lty=2)# exactly: kink is at a * 2 / 1.24 = a / .62
     > curve(qgammaAppr(x - 100/.62, 100,log=TRUE), -1e-3, +1e-3)
     >
     > showProc.time()
     Time elapsed: 0.19 0 0.22
     >
     > p.qgammaLog <- function(alpha, xl.f = 1.5, xr.f = 0.4, n = 2001)
     + {
     + ## Purpose:
     + ## ----------------------------------------------------------------------
     + ## Arguments:
     + ## ----------------------------------------------------------------------
     + ## Author: Martin Maechler, Date: 30 Mar 2004, 18:44
     + kk <- -alpha / .62 # = (alpha * 2) / (-1.24)
     + curve(qgamma(x, alpha, log=TRUE), xl.f*kk, xr.f*kk, log='y',
     + n=n, col=2, lwd=3.6, lty = 4,
     + main= paste("qgamma(x, alpha=",formatC(alpha,digits=10),", log = TRUE)"))
     + lines(kk, qgamma(kk, alpha, log=TRUE), type = 'h', lty = 3)
     + curve(qgamma (exp(x), alpha), add = TRUE, col="orange", n=n, lwd= 2)
     + curve(qgammaAppr(x, alpha, log=TRUE), add = TRUE, col=3, n=n,lwd = .4)
     + }
     >
     > showProc.time()
     Time elapsed: 0 0 0
     >
     > p.qgammaLog(25)
     > p.qgammaLog(16)# ~ [-25, -20]
     > p.qgammaLog(12, 1.2, 0.8)# small problem remaining
     > p.qgammaLog(11, 1.2, 0.8)# even smaller
     > p.qgammaLog(10.5, 1.1, 0.9)# even smaller
     > p.qgammaLog(10.25, 1.1, 0.9)# even smaller
     > ## 2019-08: __nothing__ visible from here on:
     > p.qgammaLog(10.18, 1.02, 0.98)# even smaller
     > p.qgammaLog(10.15, 1.02, 0.98)# even smaller
     > p.qgammaLog(10.14, 1.001, 0.999)# even smaller
     > p.qgammaLog(10.139, 1.0002, 0.9998)#
     > p.qgammaLog(10.138, 1.0002, 0.9998)#
     > p.qgammaLog(10.137, 1.00001, 0.99999)#
     > p.qgammaLog(10.13699, 1.0000001, 0.9999999)#
     > p.qgammaLog(10.1369899, 1.0000001, 0.9999999)#
     > p.qgammaLog(10.1369894, 1.0000001, 0.9999999)#
     > p.qgammaLog(10.1369893, 1.0000001, 0.9999999)# even smaller at -16.34998
     >
     > showProc.time()
     Time elapsed: 0.67 0.03 0.71
     >
     > ##-- here is the boundary --- for 64-bit AMD Opteron ---
     > ## and for 32-bit AMD Athlon
     >
     > p.qgammaLog(10.1369892, 1.0000001, 0.9999999)# no more
     > p.qgammaLog(10.136989, 1.0000001, 0.9999999)#
     > p.qgammaLog(10.136988, 1.0000001, 0.9999999)#
     > p.qgammaLog(10.136985, 1.0000001, 0.9999999)#
     > p.qgammaLog(10.13698, 1.0000001, 0.9999999)#
     > p.qgammaLog(10.13697, 1.0000001, 0.9999999)#
     > p.qgammaLog(10.13695, 1.0000001, 0.9999999)#
     > p.qgammaLog(10.1368, 1.000001, 0.999999)#
     > p.qgammaLog(10.1365, 1.000001, 0.999999)#
     > p.qgammaLog(10.136, 1.000001, 0.999999)#
     > p.qgammaLog(10.125, 1.1, 0.9)# --- see it now
     > p.qgammaLog(10, 1.2, 0.8)
     > p.qgammaLog(9)
     >
     > showProc.time()
     Time elapsed: 0.52 0.03 0.61
     >
     > ## For large alpha: show difference to see problem better
     > ## ---> for alpha >= 10, the x problem starts *roughly* at x = -0.8*alpha
     > ##
     >
     > sfsmisc::mult.fig(2)
     > curve(qgammaAppr(x, 5, log=TRUE), - 8.1, -8, n=2001)
     > curve(qgammaAppr(x- 5/.62, 5, log=TRUE), -1e-15, 0)
     >
     > ## is the kink from pgamma() ? : no: this looks fine,
     > curve(pgamma(x, 1e5, log=TRUE), 1, 2e5, log='x', n=2001,col=2)
     > ## and this does too:
     > curve( dgamma(x, 1e5), .5e5, 2e5); par(new=TRUE)
     > curve( dgamma(x, 1e5, log=TRUE), .5e5, 2e5, col=2, yaxt="n")
     > axis(4,col.axis=2); par(new=TRUE)
     > curve( pgamma(x, 1e5), .5e5, 2e5, n=2001, col=3); par(new=TRUE)
     > curve( pgamma(x, 1e5, log=TRUE), .5e5, 2e5, n=2001, col=4); par(new=TRUE)
     > curve(-pgamma(x, 1e5, log=TRUE,lower=FALSE), .5e5, 2e5, n=2001, col=4)
     > ## all looking nice
     >
     >
     > x <- 10^seq(2,6, length=4001)
     > qx <- qgamma(pgamma(x, 1e5, log=TRUE), 1e5, log=TRUE)
     > plot(x, qx, type ='l', col=2, asp = 1); abline(0,1, lty=3)
     >
     > showProc.time()
     Time elapsed: 0.09 0.02 0.13
     > <0c>
     > ###------------- Approximations of qgamma() ------
     > ##
     >
     > ## source("/u/maechler/R/MM/NUMERICS/dpq-functions/qchisqAppr.R")
     > ##--> qchisqAppr()
     > ##--> qchisqWH [ = Wilson Hilferty ]
     > ##--> qchisqKG [ = Kennedy & Gentle's improvements "a la AS 91" ]
     > ## dyn.load("/u/maechler/R/MM/NUMERICS/dpq-functions/qchisq_appr.so")
     >
     > ## Consider the two different implementations of
     > ## lgamma1p(a) := lgamma(1+a) == log(gamma(1+a) == log(a*gamma(a)) "stable":
     >
     > if(!exists("lseq", mode="function"))
     + lseq <- if(requireNamespace("sfsmisc")) sfsmisc::lseq else
     + function(from, to, length) exp(seq(log(from), log(to), length.out = length))
     >
     > if(require("Rmpfr")) { ##---------------- MPFR numbers -------------------------
     +
     + .mpfr.all.eq <- Rmpfr::all.equal
     + AllEq <- function(target, current, ...)
     + .mpfr.all.eq(target, current, ...,
     + formatFUN = function(x, ...) Rmpfr::format(x, digits = 9))
     +
     + print(gammaE <- Const("gamma",200)); pi. <- Const("pi",200)
     + print(a0 <- (gammaE^2 + pi.^2/6)/2)
     + print(psi2.1 <- -2*zeta(mpfr(3,200)))# == psigamma(1,2) =~ -2.4041138
     + print(a1 <- (psi2.1 - gammaE*(pi.^2/2 + gammaE^2))/6)
     +
     + x <- lseq(1e-30, 0.8, length = if(doExtras) 1000 else 125)
     + x. <- mpfr(x, 200)
     + xct. <- log(x. * gamma(x.)) ## using MPFR arithmetic .. no overflow ...
     + xc2. <- log(x.) + lgamma(x.)## (ditto)
     + print(AllEq(xct., xc2., tol = 0)) # 3.15779......e-57
     + xct <- as.numeric(xct.)
     + stopifnot(exprs = {
     + AllEq(xct., xc2., tol = 1e-45)
     + AllEq(xct , xc2., tol = 1e-15)
     + ##
     + all.equal(lgamma1p(x), lgamma1p(x, tol= 1e-16), tol=0)
     + ## -> no difference; i.e., default tol = 1e-14 seems fine enough!
     + })
     + showProc.time()
     +
     + m.appr <- cbind(log(x*gamma(x)), lgamma(1+x), log(x) + lgamma(x),
     + lgamma1p.(x, k=1, cut=3e-6),
     + lgamma1p.(x, k=2, cut=1e-4),
     + lgamma1p.(x, k=3, cut=8e-4),
     + lgamma1p(x))#, tol= 1e-14), # = default
     +
     + eMat <- m.appr - xct # absolute error
     + ## Relative errors:
     + str(reMat. <- m.appr/xct. - 1)
     + str(reMat <- as(reMat., "array")) # as(., "matrix") fails in older versions
     +
     + matplot(x, eMat , log="x", type="l", lty=1) #-> problematic log(x) + lgamma(x) for "large"
     + matplot(x, abs( eMat), log="xy", type="l", lty=1) #-> but good for small; lgamma1p is much better
     + matplot(x, abs(reMat), log="xy", type="l", lty=1)
     + abline(v= 3.47548562941137e-08, col = "gray80", lwd=3)#<- the cutoff value of lgamma1p()
     + ##---> should use earlier cutoff!
     + ## zoom in:
     + matplot(x, abs(reMat), log="xy", type="l", lty=1, xlim=c(8e-9, 1e-3))
     + abline(v= 3.47548562941137e-08, col = "gray80", lwd=3)#<- the cutoff value of lgamma1p()
     +
     + ## rm(x., xct., xc2., reMat., eMat, AllEq)
     + detach("package:Rmpfr")
     + showProc.time()
     +
     + } ## if( MPFR ) ----------------------------------------------------------------
     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 32 bits per limb
    
    
     Attaching package: 'Rmpfr'
    
     The following object is masked from 'package:gmp':
    
     outer
    
     The following objects are masked from 'package:stats':
    
     dbinom, dnorm, dpois, pnorm
    
     The following objects are masked from 'package:base':
    
     cbind, pmax, pmin, rbind
    
     1 'mpfr' number of precision 200 bits
     [1] 0.5772156649015328606065120900824024310421593359399235988057672
     1 'mpfr' number of precision 200 bits
     [1] 0.9890559953279725553953956515006347079391835207282140904431957
     1 'mpfr' number of precision 200 bits
     [1] -2.404113806319188570799476323022899981529972584680997763584544
     1 'mpfr' number of precision 200 bits
     [1] -0.9074790760808862890165601673562751149286114490725637609413306
     [1] "Mean relative difference: 3.181053350957963600882729201783933343641695627742022642032500e-57"
     Time elapsed: 1.87 0.03 1.96
     Class 'mpfrArray' [package "Rmpfr"] of dimension c(125, 7) and precision 200
     -1 -1 ...
     num [1:125, 1:7] -1.00 -1.00 6.34e+13 3.64e+13 -1.00 ...
     - attr(*, "dimnames")=List of 2
     ..$ : NULL
     ..$ : NULL
     Time elapsed: 0.08 0 0.1
     Warning messages:
     1: In xy.coords(x, y, xlabel, ylabel, log = log) :
     348 y values <= 0 omitted from logarithmic plot
     2: In xy.coords(x, y, xlabel, ylabel, log) :
     1 y value <= 0 omitted from logarithmic plot
     >
     >
     > ## ../R/qchisqAppr.R -- talks about the "small shape" qgamma() approxmation
     > ## ----------------- --> .qgammaApprBnd() :
     > curve(.qgammaApprBnd, 1e-18, 1e-15, col=2)
     > abline(h=0, col="gray70", lty=2)
     > eps.c <- .Machine$double.eps
     > axis(3,at=(1:3)* eps.c,
     + label=expression(epsilon[c], 2*epsilon[c], 3*epsilon[c]))
     > (rt.b <- uniroot(.qgammaApprBnd, c(1,3)*eps.c, tol=1e-12))
     $root
     [1] 2.220446e-16
    
     $f.root
     [1] 1.281676e-16
    
     $iter
     [1] 0
    
     $init.it
     [1] NA
    
     $estim.prec
     [1] 4.440892e-16
    
     > rt.b$root ## 3.954775e-16
     [1] 2.220446e-16
     > rt.b$root / eps.c ## 1.781072
     [1] 1
     > ##==> for a < 1.781*eps, bnd > 0 ==> we have log(p) < bnd for all p
     > ## otherwise, we should effectively 'test'
     > curve(.qgammaApprBnd, 1e-16, 1e-10, log="x", col=2)
     > showProc.time()
     Time elapsed: 0.02 0 0.01
     >
     >
     > ## source ("/u/maechler/R/MM/NUMERICS/dpq-functions/beta-gamma-etc/qgamma-fn.R")
     > ## ##--> qchisqAppr.R() -- which has 'kind = ' argument!
     > ## ##--> qgamma.R()
     >
     > p.qchi.appr <-
     + function(x, qm= { m <- cbind(qchisq(x, df, log=TRUE),
     + sapply(knds, function(kind)
     + qchisqAppr.R(x,df,log=TRUE,kind=kind)))
     + colnames(m) <- c("True", "default", knds[-1])
     + m },
     + df,
     + knds = list(NULL,"chi.small", "WH", "p1WH", "df.small"),
     + call = match.call(), main = deparse(call), log = "y", ...)
     + {
     + ## Purpose:
     + ## ----------------------------------------------------------------------
     + ## Arguments:
     + ## ----------------------------------------------------------------------
     + ## Author: Martin Maechler, Date: 25 Mar 2004, 22:08
     +
     + col <- c(2,1,3:6)
     + lty <- c(1,3,1,1,1,1)
     + lwd <- c(2,2,1,1,1,1)
     + matplot(x, qm, col=col, lty=lty, lwd=lwd, log = log,
     + type = 'l', main = main, ...)
     + y0 <- c( .02, .98) %*% par('usr')[3:4]
     + if(par("ylog")) y0 <- 10^y0
     + legend(min(x), y0, colnames(qm), col=col, lty=lty, lwd=lwd)
     + invisible(list(x=x, qm=qm, call=match.call()))
     + }
     >
     >
     > pD.chi.appr <- function(pqr, err.kind=c("relative", "absolute"),
     + type = "l", log = "y",
     + lwds = c(2, rep(1, k-1)),
     + cols = seq(along=lwds),
     + ltys = rep(1,k),
     + ...)
     + {
     + ## Purpose: Plot Difference from "True" qchisq()
     + ## ----------------------------------------------------------------------
     + ## Arguments: pqr: a list as resulting from p.chi.appr()
     + ## ----------------------------------------------------------------------
     + ## Author: Martin Maechler, Date: 31 Mar 2004, 09:38
     + err.kind <- match.arg(err.kind)
     + if(!is.list(pqr) || !is.numeric(k <- ncol(pqr$qm)-1) || k <= 0)
     + stop("invalid first argument 'pqr'")
     + with(pqr, {
     + err <- abs(if(err.kind == "relative")
     + (1- qm[,-1] / qm[,1]) else (qm[,-1] - qm[,1]))
     + matplot(x, err, ylab = "",
     + main = paste(err.kind,"Error from", deparse(call)),
     + type= type, log= log, lty=ltys, col=cols, lwd=lwds, ...)
     + legend(par("xaxp")[1], par("yaxp")[2], colnames(qm)[-1],
     + lty=ltys, col=cols, lwd=lwds)
     + })
     + }
     >
     > ## if(FALSE) # you can manually set
     > ## do.pqchi <- TRUE # before source()ing this file
     > ## if(!exists("do.pqchi") || !is.logical(do.pqchi))
     > ## do.pqchi <- interactive()
     >
     > ## if(do.pqchi) { #------------- FIXME look at speed up or cache indeed !?? <<<<<
     >
     > ## pqFile <- "/u/maechler/R/MM/NUMERICS/dpq-functions/pqchi.rda"
     > ## ## ls -l ................... 1325446 Nov 2 2009 pqchi.rda
     > ## if(file.exists(pqFile)) {
     > ## attach(pqFile) ## it loads more than we create here __FIXME__
     > ## print(ls(2, all.names=TRUE))
     > ## ## [1] "pq.1" "pq.25" "pq.25." "pq.31" "pq.33" "pq.33." "pq.33.2"
     > ## ## [8] "pq.33.3" "pq.33.4" "pq.5" "pq.5." "pq1" "pq1." "pq1.95"
     > ## ## [15] "pq1.95." "pq1.95.2" "pq10" "pq10." "pq10.2" "pq100" "pq2"
     > ## ## [22] "pq2." "pq2.05" "pq2.05." "pq2.05.2" "pq2.5" "pq2.5." "pq2.5.2"
     > ## ## [29] "pq200" "pq2L" "pq4" "pq4." "pq4.2" "pqFile"
     > ## }
     >
     > showProc.time()
     Time elapsed: 0 0.02 0.02
     > x <- seq(-300, -.01, length=501)
     >
     > all.equal(qchisqAppr.R(x, 200, log=TRUE),
     + qchisqAppr (x, 200, log=TRUE),tol=0)
     [1] "Mean relative difference: 3.820202e-16"
     > ## 4.48 e-16 / TRUE (Opteron)
     >
     > all.equal(qchisqAppr.R(x, 2, log=TRUE),
     + qchisqAppr (x, 2, log=TRUE),tol=0)
     [1] "Mean relative difference: 3.903776e-16"
     > ## 3.90 e-16 / TRUE (Opteron)
     >
     > all.equal(qchisqAppr.R(x, 0.1, log=TRUE),
     + qchisqAppr (x, 0.1, log=TRUE),tol=0)
     [1] "Mean relative difference: 1.72789e-15"
     > ## 7.15 e-15 / 1.179e-8 !!!!! (Opteron)
     >
     > pq200 <- p.qchi.appr(x = seq(-300, -.01, length=501), df = 200)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     939 y values <= 0 omitted from logarithmic plot
     > pq100 <- p.qchi.appr(x = seq(-160, -.01, length=501), df = 100)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     941 y values <= 0 omitted from logarithmic plot
     > ## after (slow) computing, quickly repeat:
     > with(pq200, p.qchi.appr(x=x, qm=qm, call=call))
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     939 y values <= 0 omitted from logarithmic plot
     > with(pq100, p.qchi.appr(x=x, qm=qm, call=call))
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     941 y values <= 0 omitted from logarithmic plot
     >
     > ## this "hangs forever" -- before I introduced 'maxit' (for 'nu.small'):
     > pq10 <- p.qchi.appr(x = seq(-12, -.005, length=501), df = 10)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     891 y values <= 0 omitted from logarithmic plot
     > ## want to see the jump:
     > pq10. <- p.qchi.appr(x = seq(-10, -4, length=501), df = 10)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     1002 y values <= 0 omitted from logarithmic plot
     > pq10.2<- p.qchi.appr(x = seq(-8.5,-7.5, length=501), df = 10)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     1002 y values <= 0 omitted from logarithmic plot
     > with(pq10.2, p.qchi.appr(x=x, qm=qm, call=call))
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     1002 y values <= 0 omitted from logarithmic plot
     >
     >
     > pq2.5 <- p.qchi.appr(x = seq(-3.4, -.01, length=501), df = 2.5)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     244 y values <= 0 omitted from logarithmic plot
     > pq2.5. <- p.qchi.appr(x = seq(-2.1, -1.8, length=901), df = 2.5)#the jump
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     901 y values <= 0 omitted from logarithmic plot
     > ## what about p1WH (which is fantastic for df=2)?
     > pq2.5.2<- p.qchi.appr(x = seq(-0.5, -1e-3, length=901), df = 2.5)
     > with(pq2.5, p.qchi.appr(x=x, qm=qm, call=call))
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     244 y values <= 0 omitted from logarithmic plot
     > with(pq2.5.2, p.qchi.appr(x=x, qm=qm, call=call))
     > pD.chi.appr(pq2.5.2)# nothing special
     >
     > pq2.05 <- p.qchi.appr(x = seq(-3.4, -.01, length=501), df = 2.05)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     81 y values <= 0 omitted from logarithmic plot
     > pq2.05. <- p.qchi.appr(x = seq(-2.5, -1.5, length=901), df = 2.05)#the jump
     > ## ^^ the jump from chi.small to WH is much too late here
     > ## what about p1WH (which is fantastic for df=2)?
     > pq2.05.2<- p.qchi.appr(x = seq(-0.4, -1e-5, length=901), df = 2.05)
     > pD.chi.appr(pq2.05.2) # p1WH is starting to become better
     > # and the jump (WH -> p1WH) is too late
     >
     > with(pq2.05, p.qchi.appr(x=x, qm=qm, call=call))
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     81 y values <= 0 omitted from logarithmic plot
     > with(pq2.05.2, p.qchi.appr(x=x, qm=qm, call=call))
     >
     >
     > ## Here, all are 'ok' (but "nu.small"):
     > pq2L <- p.qchi.appr(seq(-300, -.01, length=201), df = 2)
     There were 50 or more warnings (use warnings() to see the first 50)
     >
     > pq2 <- p.qchi.appr(x = seq(-5, -.001, length=501), df = 2)
     > pq2. <- p.qchi.appr(x = seq(-2.5, -1, length=901), df = 2)
     > with(pq2., p.qchi.appr(x=x, qm=qm, call=call))
     >
     > pq4 <- p.qchi.appr(x = seq(-8, -0.01, length = 501), df = 4)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     845 y values <= 0 omitted from logarithmic plot
     > summary(warnings())
     1 identical warnings:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     845 y values <= 0 omitted from logarithmic plot
     > pq4.2 <- p.qchi.appr(x = seq(-0.1, -1e-04, length = 901), df = 4)
     >
     > pq1.95 <- p.qchi.appr(x = seq(-3., -.01, length=501), df = 1.95)
     > pq1.95. <- p.qchi.appr(x = seq(-2.2, -1.5, length=901), df = 1.95)#the jump -1.57
     > ## ^^ the jump from chi.small to WH is *much* too late here
     > ## what about p1WH (which is fantastic for df=2)?
     > pq1.95.2<- p.qchi.appr(x = seq(-0.4, -1e-7, length=901), df = 1.95)
     > pD.chi.appr(pq1.95.2) # p1WH is starting to become better
     > # and the jump (WH -> p1WH) is too late
     >
     > with(pq1.95, p.qchi.appr(x=x, qm=qm, call=call))
     > with(pq1.95.2, p.qchi.appr(x=x, qm=qm, call=call))
     >
     >
     > pq1 <- p.qchi.appr(x = seq(-4, -.001, length=501), df = 1)
     There were 50 or more warnings (use warnings() to see the first 50)
     > pq1. <- p.qchi.appr(x = seq(-1.8, -.5, length=901), df = 1)
     > with(pq1., p.qchi.appr(x=x, qm=qm, call=call))
     >
     > pq.5 <- p.qchi.appr(x = seq(-1.5, -.001, length=501), df = 0.5)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     243 y values <= 0 omitted from logarithmic plot
     > pq.5. <- p.qchi.appr(x = seq(-0.8, -0.2, length=901), df = 0.5, ylim=c(.04,.6))
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     40 y values <= 0 omitted from logarithmic plot
     >
     > pq.33 <- p.qchi.appr(x= seq(-0.9, -.001,length=501), df= 0.33)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     246 y values <= 0 omitted from logarithmic plot
     > pq.33. <- p.qchi.appr(x= seq(-0.4, -0.02,length=901), df= 0.33)
     > pq.33.2<- p.qchi.appr(x= seq(-0.4, -0.2, length=901), df= 0.33, ylim=c(.15,.60))
     > with(pq.33.2, p.qchi.appr(x=x, qm=qm, call=call,ylim=c(.15,.45)))
     >
     > pq.33.3<- p.qchi.appr(x= seq(-0.4, -0.005, length=901), df= 0.33, ylim=c(.15, 4.00))
     > with(pq.33.3, p.qchi.appr(x=x, qm=qm, call=call))#,ylim=c(.15, 8)))
     >
     > pq.33.4<- p.qchi.appr(x= seq(-0.003, -1e-6, length=901), df= 0.33,ylim=c(5,25))
     > with(pq.33.4, p.qchi.appr(x=x, qm=qm, call=call,ylim=c(5,25)))
     >
     > ## nu <= 0.32 is the "magic" border of Best & Roberts
     >
     > pq.31 <- p.qchi.appr(x = seq(-0.45,-.010, length=501), df = 0.31)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     27 y values <= 0 omitted from logarithmic plot
     > with(pq.31, p.qchi.appr(x=x, qm=qm, call=call))
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     27 y values <= 0 omitted from logarithmic plot
     >
     > pq.25 <- p.qchi.appr(x = seq(-0.3, -0.02, length=901), df = 0.25)
     >
     > pq.1 <- p.qchi.appr(x = seq(-0.16, -0.01, length=901), df = 0.1)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     226 y values <= 0 omitted from logarithmic plot
     > with(pq.1, p.qchi.appr(x=x, qm=qm, call=call))
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     226 y values <= 0 omitted from logarithmic plot
     > showProc.time()
     Time elapsed: 7.92 0.09 8.28
     >
     > ## if(!file.exists(pqFile)) # don't overwrite for now (as it contains pq2L ,
     > ## save(list=ls(pat="^pq"), file = pqFile)
     >
     > ##}## end if(do.pqchi){ only if interactive } ======================================
     >
     > pD.chi.appr(pq2L, "abs")
     Warning messages:
     1: In xy.coords(x, y, xlabel, ylabel, log = log) :
     363 y values <= 0 omitted from logarithmic plot
     2: In xy.coords(x, y, xlabel, ylabel, log) :
     180 y values <= 0 omitted from logarithmic plot
     > pD.chi.appr(pq2L, "rel")
     Warning messages:
     1: In xy.coords(x, y, xlabel, ylabel, log = log) :
     363 y values <= 0 omitted from logarithmic plot
     2: In xy.coords(x, y, xlabel, ylabel, log) :
     180 y values <= 0 omitted from logarithmic plot
     > ## --> want only much smaller x-range:
     > pD.chi.appr(pq2,"abs")#--> fantastic p1WH
     Warning messages:
     1: In xy.coords(x, y, xlabel, ylabel, log = log) :
     363 y values <= 0 omitted from logarithmic plot
     2: In xy.coords(x, y, xlabel, ylabel, log) :
     1 y value <= 0 omitted from logarithmic plot
     > pD.chi.appr(pq2) # (ditto)
     Warning messages:
     1: In xy.coords(x, y, xlabel, ylabel, log = log) :
     363 y values <= 0 omitted from logarithmic plot
     2: In xy.coords(x, y, xlabel, ylabel, log) :
     1 y value <= 0 omitted from logarithmic plot
     >
     > pD.chi.appr(pq4.2)# p1WH: only at very right
     > pD.chi.appr(pq4.2, xlim=c(-.016,0))# p1WH: only at very right
     >
     >
     > ## no Newton step here, eg:
     > (qgA100 <- qgammaAppr(1e-100, 100))
     [1] 3.799269
     > (qg.100 <- qgamma (1e-100, 100))
     [1] 3.950799
     > all.equal(qgA100, qg.100)
     [1] "Mean relative difference: 0.03988396"
     > ## too much different
     > dgamma(1e-100, 100, log = TRUE)# -23154.7 i.e., "non-log" is 0
     [1] -23154.73
     >
     > qgamma.R(1e-100, 100, verbose = TRUE)#-> final Newton fails!
     qgamma(p= 1e-100, alpha= 100, scale= 1, l.t.= 1, log.p= 0): (p1,df)=(-230.259,200) ==> small chi-sq., ch0 = 7.59854
    
     ==> ch = 7.59854: Ph.II iter; ch=7.59854, p2=9.76738e-101
     it=2, ch=2.45916e+09, p2=-1
     --> Del became NA
    
     it=1: p=-230.259, x = 3.79927, p.=-234.019; p1:=D{p}=-3.76094
     new t= 3.94774005, p.= -230.332933 it=2, d{p}=-0.074424
     new t= 3.95079758, p.= -230.258539 it=3, d{p}=-2.99766e-05
     new t= 3.95079881, p.= -230.258509 it=4, d{p}=-4.88853e-12
     new t= 3.95079881, p.= -230.258509 it=5, d{p}=0
     [1] 3.950799
     >
     > ## but here, the final Newton iterations do work :
     > x <- qgamma.R(log(1e-100), 100, log = TRUE, verbose = TRUE)
     qgamma(p=-230.259, alpha= 100, scale= 1, l.t.= 1, log.p= 1): (p1,df)=(-230.259,200) ==> small chi-sq., ch0 = 7.59854
    
     it=1: p=-230.259, x = 3.79927, p.=-234.019; p1:=D{p}=-3.76094
     new t= 3.94774005, p.= -230.332933 it=2, d{p}=-0.074424
     new t= 3.95079758, p.= -230.258539 it=3, d{p}=-2.99766e-05
     new t= 3.95079881, p.= -230.258509 it=4, d{p}=-4.88853e-12
     new t= 3.95079881, p.= -230.258509 it=5, d{p}=0
     > pgamma(x, 100) # = 1e-100 ! perfect
     [1] 1e-100
     > showProc.time()
     Time elapsed: 0.11 0 0.11
     >
     > ###--> Use this to devise an improved final Newton algorithm !!!
     >
     > <0c>
     > ## From: Prof Brian Ripley <ripley@stats.ox.ac.uk>
     > ## To: skylab.gupta@gmail.com
     > ## Cc: R-bugs@biostat.ku.dk, r-devel@stat.math.ethz.ch
     > ## Subject: Re: [Rd] qgamma inaccuracy (PR#12324)
     > ## Date: Tue, 12 Aug 2008 20:50:50 +0100 (BST)
     >
     > ## This is a really extreme usage. AFAICS the code works well enough down to
     > ## shape=1e-10 or so, e.g.
     >
     > qgamma(1e-10, 5e-11, lower.tail=FALSE)
     [1] 0.08237203
     > ## [1] 0.08237203
     > ## in R 2.9.. .. 2.10.0 -- with an accuracy warning {which is *wrong*!}
     >
     > ## I would be interested to know what substantive problem you were trying to
     > ## solve here that required such values.
     >
     > ## I am pretty sure that a completely different algorithm will be required.
     >
     > ## MM: It looks like this is basis for a new algo:
     > a <- 1e-14
     > gammaE <- 0.57721566490153286060651209008240243079 # Euler's gamma
     > curve(pgamma(x, a, lower.tail=FALSE)/a + log(x) + gammaE, 1e-300, 1e-1, log="",n=1000)
     > curve(pgamma(x, a, lower.tail=FALSE)/a + log(x) + gammaE - x, 1e-300, 1e-1, log="",n=1000)
     > ## ==> Q = 1 - P = pgamma(x,a, lower=FALSE) ~= a*(-log(x) - gammaE + x - x^2/4)
     > ## i.e., Q ~= -a(log(x) + gammaE { -x + x^2/4 }
     > ## -Q/a - gammaE ~= log(x) { -x + x^2/4 }
     > ## ==> x ~= exp(- (Q/a + gammaE))
     > ## e.g., example below:
     > Q <- 1e-100; a <- 5e-101
     > ## MM: Find inverse :
     > str(r.a <- uniroot(function(x) pgamma(x,a,lower.tail=FALSE) - Q,
     + int = c(0.01, 0.1), tol=1e-20))
     List of 5
     $ root : num 0.0824
     $ f.root : num -1.27e-116
     $ iter : int 8
     $ init.it : int NA
     $ estim.prec: num 4.16e-17
     > dput(x0 <- r.a$root) ## 0.0823720296207203
     0.0823720296207203
     > (x1 <- exp(- (Q/a + gammaE)))## 0.07598528 .. not so good
     [1] 0.07598528
     > qgammaApprSmallP(Q, a, lower.tail=FALSE)## ~= 0.07598528 -- the same!!
     [1] 0.07598528
     >
     > pgamma(x0, a, lower.tail=FALSE) ## 1.00000e-100.
     [1] 1e-100
     > pgamma(x1, a, lower.tail=FALSE) ## 1.03728e-100 ... hmm "close"
     [1] 1.037283e-100
     > ##
     >
     > ## MM: -- now look at the bigger picture
     > p.qg.2a <- function(l2x.min= -15, l2x.max = -100, n = 501,
     + do.offset = FALSE,
     + type = "o", log = "x", cex = 0.6, ...) {
     + x.log <- any("x" == strsplit(log,"")[[1]])
     + x <- if(x.log) 2^seq(l2x.min, l2x.max, length=n)
     + else seq(2^l2x.min,2^l2x.max, length=n)
     + if(do.offset)
     + plot(x, qgamma(2*x, x, lower.tail=FALSE) - 0.0823720296206873,
     + type=type, cex=cex, log=log, ...)
     + else plot(x, qgamma(2*x, x, lower.tail=FALSE),
     + type=type, cex=cex, log=log, ...)
     + }
     > p.qg.2a() # was "very bad" in R <= 2.10.0, now --> 0.082372... "perfect" smooth
     > ## still a little ---but acceptable--- remaining inaccuracy ...zooming in:
     > p.qg.2a(-43,-55, do.offset=TRUE)
     > p.qg.2a(,-1024)
     > p.qg.2a(,-1024, log="", pch=".")## linear in x !!
     > ## zoom in at the limit
     > p.qg.2a(-30,-1024, do.offset=TRUE, ylim = 1e-11*c(-1,1))
     > p.qg.2a(-33,-1024, do.offset=TRUE, ylim = 1e-12*c(-1,1))
     > p.qg.2a(-33,-1024, do.offset=TRUE, ylim = 1e-13*c(-1,1))
     >
     > a <- 2^-(7:900)
     > qg <- qgamma(2*a, a, lower.tail=FALSE)
     > re <- 1-pgamma(qg, a, lower.tail=FALSE)/(2*a)
     > plot(a, re, log="x", type="b", col=2)
     > stopifnot(abs(re) < 2e-12) # but really, *should be a bit better
     >
     > showProc.time()
     Time elapsed: 0.11 0.02 0.16
     > ## For completeness we may write that in due course, but for now (R 2.7.2) I
     > ## suggest just issuing a warning for miniscule 'shape'.
     >
     > ## On Thu, 7 Aug 2008, skylab.gupta@gmail.com wrote:
     >
     > ## > Full_Name:
     > ## > Version: 2.7.1 (2008-06-23)
     > ## > OS: windows vista
     > ## > Submission from: (NULL) (216.82.144.137)
     > ## >
     > ## >
     > ## > Hello,
     > ## >
     > ## > I have been working with various probability distributions in R, and it seems
     > ## > the gamma distribution is inaccurate for some inputs.
     >
     > ## > For example, qgamma(1e-100, 5e-101, lower.tail=FALSE) gives: 1.0. However, it
     > ## > seems this is incorrect; I think the correct answer should be
     > ## > 0.082372029620717283. When I check these numbers using pgamma, I get:
     >
     > (qg <- qgamma(1e-100, 5e-101, lower.tail=FALSE))# 1 (wrong, originally)
     [1] 0.08237203
     > ## 0.08237203 now (2009-11-04), i.e. ok
     > pgamma(qg, 5e-101, lower.tail=FALSE)# now -> 1e-100 : ok
     [1] 1e-100
     >
     > pgamma(0.082372029620717283,5e-101, lower.tail=FALSE)
     [1] 1e-100
     > ## 1.0000000000000166e-100.
     >
     > RE.pqgamma <- function(p, shape, lower.tail = TRUE, log.p = FALSE) {
     + ## Relative Error of pgamma(qgamma(*), ..):
     + 1 - pgamma(qgamma(p, shape, lower.tail=lower.tail, log.p=log.p),
     + shape=shape, lower.tail=lower.tail, log.p=log.p) / p
     + }
     > RE.qpgamma <- function(q, shape, lower.tail = TRUE, log.p = FALSE) {
     + ## Relative Error of qgamma(pgamma(*), ..):
     + 1 - qgamma(pgamma(q, shape, lower.tail=lower.tail, log.p=log.p),
     + shape=shape, lower.tail=lower.tail, log.p=log.p) / q
     + }
     >
     > ## Ok, how extreme can we get -- let a := alpha := shape --> 0 :
     > x <- 1e-100
     > a <- 2^-(7:300)# is still "ok":
     > plot(a, (re <- RE.pqgamma(x, a, lower.tail=FALSE)), type="b", col=2, log="x")# oops!
     > a <- 2^-(7:400)
     > plot(a, (re <- RE.pqgamma(x, a, lower.tail=FALSE)), type="b", col=2, log="x")# oops!
     > ## Oops!
     > ## but, it is clear
     > qgamma(x, 2^-400, lower.tail=FALSE)## is exactly 0
     [1] 0
     >
     > ## -> it goes to 0 quickly . .. zooming in:
     > curve(qgamma(1e-100, x, lower.tail=FALSE), 1e-110, 1e-70, log="xy", col=2, axes=FALSE)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     18 y values <= 0 omitted from logarithmic plot
     > eaxis(1);eaxis(2)
     >
     > ## from when on is it exactly 0:
     > uniroot(function(u) qgamma(1e-100, 2^u,lower.tail=FALSE)-1e-315, c(-400, -300))$root
     [1] -341.6941
     > ## -341.6941
     > ## => use
     > a <- 2^-(7:341)
     > plot(a, (re <- RE.pqgamma(x, a, lower.tail=FALSE)),
     + type="l", col=2, log="x")# small glips
     > ## zoom in:
     > curve(abs(RE.pqgamma(1e-100, x, lower.tail=FALSE)), 2^-341, 1e-90, log="xy", n=2000)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     125 y values <= 0 omitted from logarithmic plot
     > curve(RE.pqgamma(1e-100, x, lower.tail=FALSE), 1e-100, 10e-100, n=2000)
     > curve(RE.pqgamma(1e-100, x, lower.tail=FALSE), 4e-100, 6e-100, n=2000)
     > ## Ok: at least here is a problem
     > RE.pqgamma(1e-100, 5e-100, lower.tail=FALSE)# -0.1538
     [1] 1.687539e-14
     >
     > ## more general
     > curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 1e-100, 1e-20, n=10000, log="x")
     > ## problem *everywhere* , starting quite early: (lesser problem at ~ 1e-16 !)
     > curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 1e-25, 1e-15, n=1000, log="x")
     > curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 1e-21, 10e-21,n=1000)
     > curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 2e-21, 6e-21, n=1000)
     > curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 4e-21, 4.5e-21, n=1000)
     > ## and indeed, it's qgamma() that jumps here
     > curve(qgamma(x, 5*x, lower.tail=FALSE), 4e-21, 4.5e-21, ylim=c(.97, 1.1))
     >
     > ## well, looking at pgamma(), finally reveals the buglet is there first:
     > ## There's a jump at x = 1 !!
     > curve(pgamma(x, 1e-30, lower=FALSE), .9999, 1.0001)
     > curve(pgamma(x, 1e-20, lower=FALSE), .9999, 1.0001)
     > curve(pgamma(x, 1e-17, lower=FALSE), .9999, 1.0001)
     > curve(pgamma(x, 1e-15, lower=FALSE), .9999, 1.0001)
     > curve(pgamma(x, 1e-13, lower=FALSE), .9999, 1.0001)
     > curve(pgamma(x, 1e-12, lower=FALSE), .9999, 1.0001)# still clearly visible
     > curve(pgamma(x, 1e-11, lower=FALSE), .9999, 1.0001)# barely visible
     > ## for larger alpha == shape, must zoom in more and more:
     > curve(pgamma(x, 1e-11, lower=FALSE), .99999, 1.00001)#
     > curve(pgamma(x, 1e-10, lower=FALSE), .999999, 1.000001)
     > curve(pgamma(x, 1e-9, lower=FALSE), .9999999, 1.0000001)
     > curve(pgamma(x, 1e-8, lower=FALSE), .99999999, 1.00000001)
     > curve(pgamma(x, 1e-7, lower=FALSE), .999999999, 1.000000001)
     > curve(pgamma(x, 1e-6, lower=FALSE), .9999999999, 1.0000000001)
     > curve(pgamma(x, 1e-5, lower=FALSE), .99999999999, 1.00000000001)
     > curve(pgamma(x, 1e-4, lower=FALSE), .999999999999, 1.000000000001)
     > ## now we get close to noise level:
     > curve(pgamma(x, 1e-3, lower=FALSE), .9999999999999, 1.0000000000001)
     > curve(pgamma(x, 1e-2, lower=FALSE), .99999999999999, 1.00000000000001)
     Warning message:
     In plot.window(...) :
     relative range of values = 90 * EPS, is small (axis 1)
     >
     > showProc.time()
     Time elapsed: 0.34 0.01 0.45
     >
     > del.pgamma <- function(a, eps = 1e-13)
     + {
     + ## Purpose: *relative* jump size at x = 1 of pgamma(x, a, lower=FALSE)
     + ## ----------------------------------------------------------------------
     + ## Author: Martin Maechler, Date: 5 Nov 2009, 16:08
     + stopifnot(a > 0, length(a) == 1, eps > 0, length(eps) == 1)
     + pp <- pgamma(1+c(-eps,eps), a, lower.tail = FALSE)
     + (pp[2] - pp[1])*2/(pp[2] + pp[1])
     + }
     >
     > a <- lseq(1e-300, 1e-3, length=400)
     > dpa <- sapply(a, del.pgamma)
     > plot(a, -dpa, log="xy", type="l", col=2, yaxt="n");eaxis(2)
     > ## ok, it remains constant all the way to 1e-300
     > ## --> focus
     > a <- lseq(1e-40, 1e-5, length=400)
     > dpa <- sapply(a, del.pgamma)
     > plot(a, -dpa, log="xy", type="l", col=2, axes=FALSE)
     > eaxis(1, at = 10^-seq(5,40, by=5))
     > eaxis(2)
     >
     >
     > xm <- .Machine$double.xmin
     > pgamma(xm, shape= 1e-20)# is "practically 1" --> *most* qgamma() will be exactly 0
     [1] 1
     > ## how close to 1 ? ---> use upper_tail [possibly on log scale:]
     > pgamma(xm, shape= 1e-20, lower.tail=FALSE) # 7.078192e-18
     [1] 7.078192e-18
     > pgamma(xm, shape= 1e-20, lower.tail=FALSE, log=TRUE) # -39.48951
     [1] -39.48951
     >
     > ## Where is the 'boundary' (from which on qgamma() must return 0, since it can't give
     > ## xm = 2.2e-308 ) ?
     > a <- 2^-(7:1030)
     > plot(a, (p <- pgamma(xm, a, lower.tail=FALSE, log=TRUE)),
     + cex=.5,type ="l", col=2, log="x")
     > summary(lm. <- lm(p ~ log(a) + a + I(a^2))) ## coeff. of log(a) *is* 1
    
     Call:
     lm(formula = p ~ log(a) + a + I(a^2))
    
     Residuals:
     Min 1Q Median 3Q Max
     -0.0081566 -0.0000098 0.0000171 0.0000440 0.0107065
    
     Coefficients:
     Estimate Std. Error t value Pr(>|t|)
     (Intercept) 6.562e+00 3.327e-05 197240 <2e-16 ***
     log(a) 1.000e+00 8.024e-08 12463221 <2e-16 ***
     a -3.403e+02 2.039e-01 -1669 <2e-16 ***
     I(a^2) 1.551e+04 2.911e+01 533 <2e-16 ***
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     Residual standard error: 0.0005225 on 1020 degrees of freedom
     Multiple R-squared: 1, Adjusted R-squared: 1
     F-statistic: 5.249e+13 on 3 and 1020 DF, p-value: < 2.2e-16
    
     > summary(lm. <- lm(p ~ offset(log(a)) + a + I(a^2)))
    
     Call:
     lm(formula = p ~ offset(log(a)) + a + I(a^2))
    
     Residuals:
     Min 1Q Median 3Q Max
     -0.0081788 0.0000179 0.0000179 0.0000179 0.0107351
    
     Coefficients:
     Estimate Std. Error t value Pr(>|t|)
     (Intercept) 6.562e+00 1.639e-05 400476.3 <2e-16 ***
     a -3.404e+02 2.033e-01 -1674.4 <2e-16 ***
     I(a^2) 1.552e+04 2.907e+01 533.7 <2e-16 ***
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     Residual standard error: 0.0005232 on 1021 degrees of freedom
     Multiple R-squared: 1, Adjusted R-squared: 1
     F-statistic: 7.853e+13 on 2 and 1021 DF, p-value: < 2.2e-16
    
     >
     > p.l <- pgamma(xm, a, lower.tail=FALSE, log=TRUE) - log(a)
     > dput(mean(tail(p.l))) ## 6.56218869790132
     6.56218869790132
     > ##=> pgamma(xm, a) ~= log(a) + 6.5621886979 +
     > ## ok, to get this better, now need different a:
     > al <- seq(1e-300, 1e-3, length=200)
     > plot(al, (pl <- pgamma(xm, al, lower.tail=FALSE, log=TRUE)) - (log(al)+6.56218869790132),
     + cex=.5,type ="l", col=2)
     > summary(lm.. <- lm(pl ~ offset(log(al) + 6.56218869790132) + 0 + al + I(al^2)))
    
     Call:
     lm(formula = pl ~ offset(log(al) + 6.56218869790132) + 0 + al +
     I(al^2))
    
     Residuals:
     Min 1Q Median 3Q Max
     -1.201e-05 -4.268e-06 -1.336e-06 3.075e-06 5.247e-06
    
     Coefficients:
     Estimate Std. Error t value Pr(>|t|)
     al -3.539e+02 2.032e-03 -174123 <2e-16 ***
     I(al^2) 2.075e+04 2.617e+00 7929 <2e-16 ***
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     Residual standard error: 4.153e-06 on 198 degrees of freedom
     Multiple R-squared: 1, Adjusted R-squared: 1
     F-statistic: 1.359e+16 on 2 and 198 DF, p-value: < 2.2e-16
    
     > confint(lm..)
     2.5 % 97.5 %
     al -353.8628 -353.8547
     I(al^2) 20745.6597 20755.9814
     > coef(lm..)
     al I(al^2)
     -353.8587 20750.8205
     > ## al I(al^2)
     > ## -353.8587 20750.8205
     > ##=> pgamma(xm, a) ~= log(a) + 6.5621886979 - 353.858745 * a
     >
     > ## > Similarly, for example: -- now (2009-11-04) ok
     > qgamma(1e-100,0.005,lower.tail=FALSE) # = 109.36757177 instead of 219.5937..
     [1] 219.5937
     > pgamma(109.36757177007101, 0.005, lower.tail=FALSE)# = 1.4787306506694758e-52.
     [1] 1.478731e-52
     >
     > ## > This looks completely wrong. The correct value, I think, should be
     > ## > 219.59373661415756. In fact,
     > pgamma(219.59373661415756, 0.005, lower.tail=FALSE)# = 9.9999999999999558e-101.
     [1] 1e-100
     > ## >
     > ## > In fact, when I do the following in R, the results are completely wrong,
     > ## >
     > a <- 5*10^-(1:40)
     > z1 <- qgamma(1e-100,a,lower.tail=FALSE)
     > (y <- pgamma(z1,a,lower.tail=FALSE))
     [1] 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100
     [11] 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100
     [21] 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100
     [31] 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100
     > ## The value of y that I get should be close to 1e-100, but they are not:
     > ## [1] 1.000000e-100 1.871683e-51 1.478731e-52 1.444034e-53 1.440606e-54
     > ## [6] 1.440264e-55 1.440230e-56 1.440226e-57 1.440226e-58 1.440226e-59
     > summary(abs(1 - y/1e-100))
     Min. 1st Qu. Median Mean 3rd Qu. Max.
     0.000e+00 6.328e-15 1.293e-14 1.575e-14 2.465e-14 3.741e-14
     > plot(a, abs(1 - y/1e-100), log="xy", type="b")
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     1 y value <= 0 omitted from logarithmic plot
     > stopifnot(abs(1 - y/1e-100) < 2e-13)# max (32b Linux, P.M) = 4.186e-14
     >
     > ## > The correct values of z1 should be:
     > z1true <- c(226.97154111939946, 222.15218724493326, 219.59373661415756,
     + 217.27485383840451, 214.98015408183574, 212.68797118872064,
     + 210.39614286838227, 208.10445550564617, 205.81289009100664,
     + 203.52144711679352)
     > all.equal(z1[1:10], z1true, tol=0)# 1.307e-16 on 32-bit (Pentium M); 1.686e-16 (64b Lnx, 2019)
     [1] "Mean relative difference: 1.302676e-16"
     > stopifnot(all.equal(z1[1:10], z1true, tol = 1e-15))
     > showProc.time()
     Time elapsed: 0.16 0 0.17
     > ##>
     > ##> With these values of z1true, we get the expected values:
     >
     > (y <- pgamma(z1true, a, lower.tail=FALSE))
     [1] 1.000000e-100 1.000000e-100 1.000000e-100 1.000000e-100 1.000000e-100
     [6] 1.000000e-100 1.000000e-100 1.000000e-100 1.000000e-100 1.000000e-100
     [11] 5.869617e-112 7.428625e-111 9.705940e-111 9.970232e-111 9.997024e-111
     [16] 9.999703e-111 9.999970e-111 9.999997e-111 1.000000e-110 1.000000e-110
     [21] 5.869617e-122 7.428625e-121 9.705940e-121 9.970232e-121 9.997024e-121
     [26] 9.999703e-121 9.999970e-121 9.999997e-121 1.000000e-120 1.000000e-120
     [31] 5.869617e-132 7.428625e-131 9.705940e-131 9.970232e-131 9.997024e-131
     [36] 9.999703e-131 9.999970e-131 9.999997e-131 1.000000e-130 1.000000e-130
     > ## [1] 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100
     >
     > ## > I am using the precompiled binary version of R, under Windows Vista.
     > ## > -----------
     > ## >> version
     > ## > _
     > ## > platform i386-pc-mingw32
     > ## > arch i386
     > ## > os mingw32
     > ## > system i386, mingw32
     > ## > status
     > ## > major 2
     > ## > minor 7.1
     > ## > year 2008
     > ## > month 06
     > ## > day 23
     > ## > svn rev 45970
     > ## > language R
     > ## > version.string R version 2.7.1 (2008-06-23)
     > ## > ------------
     > ## >
     > ## > So, it seems qgamma is inaccurate for small probability values in the upper
     > ## > tail, when the shape parameter is also small.
     >
     >
     > ###_-- MM: Still wrong:
     >
     > (xm <- 2^-1074.9999) # is less than .Machine $ double.xmin == the really x > 0
     [1] 4.940656e-324
     >
     > pgamma(xm, .00001)# 0.992589
     [1] 0.992589
     > qgamma(.99, .00001)##--> NaN -- should give 0 or "xmin" or so
     [1] 0
     > ## FIXME -- ok, now
     >
     > ## but
     > curve(qgamma(x, .001, lower=FALSE), .4, .8, n=1001, log="y")
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     688 y values <= 0 omitted from logarithmic plot
     > curve(qgamma(x, 1e-5, lower=FALSE), .002, .2, n=1001, log="xy")
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     716 y values <= 0 omitted from logarithmic plot
     > curve(qgamma(x, 1e-7, lower=FALSE), 1e-5, .04, n=1001, log="xy")
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     759 y values <= 0 omitted from logarithmic plot
     > curve(qgamma(x, 1e-12, lower=FALSE), 1e-12, 1e-2, n=1001, log="xy")
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     713 y values <= 0 omitted from logarithmic plot
     >
     > ## or
     > curve(qgamma(x, 1e-121, lower=FALSE), 7e-119, 8e-119,
     + n=2001, log="y", yaxt="n")
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     1117 y values <= 0 omitted from logarithmic plot
     > try( # reveals eaxis() bug ? -- for the *subnormal* numbers
     + eaxis(2, at = 10^-seq(304,324, by=2))
     + )
     Error in eaxis(2, at = 10^-seq(304, 324, by = 2)) :
     invalid 'log=TRUE' for at <= 0: not a true log scale plot?
     >
     > curve(qgamma(x, .001, lower=FALSE), .4, .6, n=1001, log="y")
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     376 y values <= 0 omitted from logarithmic plot
     > curve(qgamma(x, .001, lower=FALSE), .5, .55, n=1001, log="y")
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     503 y values <= 0 omitted from logarithmic plot
     > ## gives a *warning* from axis() because of subnormal y-range {was error, fixed in R-devel (2018-08)}
     > curve(qgamma(x, .001, lower=FALSE), .52, .53, n=1001, log="y")
     Error in axis(side = side, at = at, labels = labels, ...) :
     CreateAtVector [log-axis()]: axp[0] = 0 < 0!
     Calls: curve ... plot.default -> localAxis -> Axis -> Axis.default -> axis
     In addition: Warning messages:
     1: In xy.coords(x, y, xlabel, ylabel, log) :
     514 y values <= 0 omitted from logarithmic plot
     2: In plot.window(...) : Internal(pretty()): very small range.. corrected
     3: In axis(side = side, at = at, labels = labels, ...) :
     CreateAtVector "log"(from axis()): axp[0] = 0 !
     Execution halted
Flavor: r-oldrel-windows-ix86+x86_64

Version: 0.3-5
Check: running tests for arch ‘x64’
Result: ERROR
     Running 'chisq-nonc-ex.R' [33s]
     Running 'dnchisq-tst.R' [1s]
     Running 'pnbeta-tst.R' [1s]
     Running 'pnt-precision-problem.R' [20s]
     Running 'ppois-ex.R' [2s]
     Running 'qPoisBinom-ex.R' [1s]
     Running 'qbeta-dist.R' [16s]
     Running 'qbeta-tst.R' [1s]
     Running 'qgamma-ex.R' [24s]
     Running 't-nonc-tst.R' [10s]
     Running 'wienergerm-accuracy.R' [9s]
     Running 'wienergerm-pchisq-tst.R' [1s]
     Running 'wienergerm_nchisq.R' [3s]
    Running the tests in 'tests/qgamma-ex.R' failed.
    Complete output:
     > library(DPQ)
     >
     > ###---> Automatically find places where qgamma() is not so precise (PR#2214) :
     > ### For PR#2214, had '1e-8' below and found quite a bit
     > ## see /u/maechler/R/MM/NUMERICS/dpq-functions/beta-gamma-etc/qgamma-ex.R ..
     >
     > ## FIXME: Timing ! --- partly these matplot() partly get quite slow ~?
     > source(system.file(package="Matrix", "test-tools-1.R", mustWork=TRUE))
     Loading required package: tools
     > ##--> showProc.time(), assertError(), relErrV(), ...
     > showProc.time()
     Time elapsed: 1.93 0.08 2.11
     >
     > (doExtras <- DPQ:::doExtras())
     [1] FALSE
     > (sdir <- system.file("safe", package="DPQ")) ## save directory (to read from)
     [1] "D:/temp/RtmpIpG5h1/RLIBS_cc0c300314ea/DPQ/safe"
     >
     > ### Nowadays finds cases in a special region for really small p and cutoff 1e-11 :
     > set.seed(47)
     > n <- if(doExtras) 100 else 32
     > res <- cbind(p=1,df=1,rE=1)[-1,]
     > for(M in 1:(if(doExtras) 20 else 10))
     + for(p in runif(n)) for(df in rlnorm(n)) {
     + r <- 1- pchisq(qchisq(p, df),df)/p
     + if(abs(r) > 1e-11) res <- rbind(res, c(p,df,r))
     + }
     >
     > ### use df in U[0,1]: finds two cases with bound 1e-11
     > for(p in runif(n)/2) for(df in runif(n)) {
     + qq <- qchisq(p, df)
     + if(qq > 0 && p > 0) {
     + r <- 1- pchisq(qq, df) / p
     + if(abs(r) > 1e-11) res <- rbind(res, c(p,df,r))
     + }
     + }
     >
     > ### now df very close to 0 : ==> finds more cases
     > for(p in sort(c(runif(64)/2, exp(-(1+rlnorm(256))))))
     + for(df in 2^-rlnorm(256, mean=2, sdlog=1.5)) {
     + qq <- qchisq(p, df)
     + if(qq > 0 && p > 0) {
     + r <- 1- pchisq(qq, df) / p
     + if(abs(r) > 1e-11) res <- rbind(res, c(p,df,r))
     + }
     + }
     > showProc.time()
     Time elapsed: 0.88 0.03 0.92
     >
     > require(graphics)
     > if(!dev.interactive(orNone=TRUE)) pdf("qgamma-appr.pdf")
     > eaxis <- sfsmisc::eaxis
     >
     > showProc.time()
     Time elapsed: 0.09 0.02 0.11
     > ## if(nrow(res) > 0) {
     > cat("Found inaccurate examples where pchisq(qchisq(p, df),df) != p\n")
     Found inaccurate examples where pchisq(qchisq(p, df),df) != p
     > ## sort in p, then df:
     > res <- res[order(res[,"p"], res[,"df"]), ]
     > rE <- res[,"rE"]
     > if(nrow(res) > 20) { hist(rE, breaks = 30); rug(rE) }
     > plot(res[,1:2])##--> quite interesting : all along one curve
     > ## p <= 1/2 and df <= 1 (about) !!
     > res <- cbind(res, nDig = round(-log10(abs(rE)), 1))
     > print(res, digits=12)
     p df rE nDig
     [1,] 0.000194375438651 0.02334079639198 6.46593024678e-08 7.2
     [2,] 0.000605300028912 0.02041606754775 -1.99857908001e-11 10.7
     [3,] 0.001012316063255 0.01855615147677 -2.59145555106e-04 3.6
     [4,] 0.001248285290785 0.01838201076117 3.94217658517e-10 9.4
     [5,] 0.001682388899865 0.01720736646288 5.53088974600e-04 3.3
     [6,] 0.001746787400790 0.01731189518997 -5.86897839217e-08 7.2
     [7,] 0.002664451237518 0.01599317398629 1.48421342013e-04 3.8
     [8,] 0.002664451237518 0.01618024201222 7.75564700239e-08 7.1
     [9,] 0.003159421860255 0.01557612780310 -7.92117005632e-06 5.1
     [10,] 0.003159421860255 0.01568183691729 -4.52237520765e-08 7.3
     [11,] 0.004055462418244 0.01493858731306 -7.17404703909e-06 5.1
     [12,] 0.004400694140827 0.01459101672970 9.07907026434e-04 3.0
     [13,] 0.004458811277768 0.01457506850867 9.03139988533e-05 4.0
     [14,] 0.004481882165743 0.01468883074316 -3.23309491179e-07 6.5
     [15,] 0.004939609905705 0.01440168350452 -2.81810098879e-06 5.6
     [16,] 0.008824465120182 0.01276352706510 1.21107345756e-04 3.9
     [17,] 0.009040265960535 0.01273711629661 1.38964402733e-05 4.9
     [18,] 0.010839089634828 0.01242499920422 2.63413624246e-10 9.6
     [19,] 0.011642124851282 0.01201471267173 -3.00271327148e-04 3.5
     [20,] 0.014753716559535 0.01155624353203 1.52962087441e-10 9.8
     [21,] 0.015499213434879 0.01125420134457 1.00492447767e-04 4.0
     [22,] 0.015499213434879 0.01135920381800 -9.55739012376e-08 7.0
     [23,] 0.018603016576955 0.01071716109330 -2.07537936112e-03 2.7
     [24,] 0.018603016576955 0.01073655493589 2.14388784340e-04 3.7
     [25,] 0.022624242394389 0.01033379525113 -3.37865757594e-09 8.5
     [26,] 0.022624242394389 0.01034206121729 4.43079705148e-08 7.4
     [27,] 0.023730217356634 0.01016252135853 -1.07732682708e-06 6.0
     [28,] 0.032427027472295 0.00942923095016 5.11205522358e-11 10.3
     [29,] 0.044753525441333 0.00839626444749 1.22224173549e-05 4.9
     [30,] 0.081818424963746 0.00686007746204 -1.78633419168e-09 8.7
     [31,] 0.081818424963746 0.00689856335721 -2.30915286892e-11 10.6
     [32,] 0.082800309102258 0.00681234719059 4.17997558788e-09 8.4
     [33,] 0.083507718914457 0.00680676700443 9.77167236016e-11 10.0
     [34,] 0.090821658072474 0.00655269761981 -7.16033632386e-09 8.1
     [35,] 0.102294760453517 0.00623563107239 -6.63889809793e-09 8.2
     [36,] 0.110869751789691 0.00603268830251 8.95578944338e-10 9.0
     [37,] 0.123950804624116 0.00571305309327 2.84683721041e-10 9.5
     [38,] 0.127405857731893 0.00562369059572 6.60541454867e-09 8.2
     [39,] 0.135229634154169 0.00540073357520 -2.34762594200e-05 4.6
     [40,] 0.137732279982451 0.00533092076413 2.99285844990e-04 3.5
     [41,] 0.138112917548194 0.00535138710974 -2.05335777981e-06 5.7
     [42,] 0.141100635980184 0.00527305771429 4.31593832968e-05 4.4
     [43,] 0.141100635980184 0.00537073537183 8.96455243371e-10 9.0
     [44,] 0.142905299416015 0.00523680041306 3.48180824883e-04 3.5
     [45,] 0.145624557210331 0.00526923971034 -1.94501770245e-09 8.7
     [46,] 0.154606872884529 0.00506806894407 -4.59924667240e-07 6.3
     [47,] 0.154606872884529 0.00507366168703 2.72301046933e-07 6.6
     [48,] 0.163535630067488 0.00497650928578 3.39664962823e-11 10.5
     [49,] 0.169741036539408 0.00484181845356 5.31400978776e-09 8.3
     [50,] 0.177327576288650 0.00465956102839 5.53404362603e-05 4.3
     [51,] 0.178169157856761 0.00471949961255 4.79807527043e-10 9.3
     [52,] 0.190094017358772 0.00450373552308 1.78565421871e-06 5.7
     [53,] 0.190147641510530 0.00453468705710 5.66235636157e-09 8.2
     [54,] 0.200112534472267 0.00442273120514 7.20473680715e-11 10.1
     [55,] 0.201518808589718 0.00439936964342 1.58748569845e-11 10.8
     [56,] 0.201518808589718 0.00439976887947 -9.97182336704e-11 10.0
     [57,] 0.210803673024037 0.00427351441034 -1.70232938856e-10 9.8
     [58,] 0.213058614771766 0.00426179831847 -1.39812605937e-11 10.9
     [59,] 0.214780951412088 0.00419869272965 -1.91879370171e-08 7.7
     [60,] 0.232805106603566 0.00395399315002 -9.17581020055e-08 7.0
     [61,] 0.249102914025652 0.00380019404026 -1.15818465929e-10 9.9
     [62,] 0.249102914025652 0.00382493512126 -1.39670497390e-11 10.9
     [63,] 0.252076511947811 0.00374903834738 2.14569900181e-07 6.7
     [64,] 0.253082914021191 0.00375259362798 3.65436092498e-09 8.4
     [65,] 0.253922058700076 0.00371237348323 -5.09014937222e-06 5.3
     [66,] 0.254289278570932 0.00374343873151 -1.05664899053e-09 9.0
     [67,] 0.260017499519858 0.00366179605930 -2.90430199223e-07 6.5
     [68,] 0.270323906831467 0.00351999192121 -1.56164756277e-04 3.8
     [69,] 0.271699356057456 0.00355068132680 5.13092990317e-09 8.3
     [70,] 0.275516196070002 0.00346804047756 -4.35171547588e-04 3.4
     [71,] 0.280722231049885 0.00348224101220 -6.07109695849e-10 9.2
     [72,] 0.284601233201101 0.00344936339590 -2.63026489478e-10 9.6
     [73,] 0.290188543054775 0.00336613521112 -5.64443074502e-08 7.2
     [74,] 0.290579022038283 0.00334423496113 1.02667567781e-07 7.0
     [75,] 0.290579022038283 0.00336764858994 2.26061565023e-08 7.6
     [76,] 0.291850198713803 0.00333552811650 -1.27338760580e-06 5.9
     [77,] 0.296521136452775 0.00330308865102 -4.48927431229e-07 6.3
     [78,] 0.298034174946132 0.00330462333485 8.42470393447e-09 8.1
     [79,] 0.300556783277253 0.00323922530004 4.66003314391e-05 4.3
     [80,] 0.303182283998467 0.00328704590597 1.82253101499e-11 10.7
     [81,] 0.303182283998467 0.00328723648952 -2.23672191879e-11 10.7
     [82,] 0.322319846303892 0.00306134512927 -1.15130830540e-05 4.9
     [83,] 0.322319846303892 0.00310689001755 8.57751647487e-11 10.1
     [84,] 0.325071272052651 0.00302343293053 -2.47088704493e-04 3.6
     [85,] 0.325071272052651 0.00304146419577 -7.71376615361e-06 5.1
     [86,] 0.331888412404218 0.00300837121343 -4.96098895297e-09 8.3
     [87,] 0.362278153188527 0.00278204202032 4.53939330569e-10 9.3
     [88,] 0.385389476781711 0.00260981704384 -1.77430203863e-09 8.8
     [89,] 0.425333956955001 0.00232995789362 1.82823025607e-08 7.7
     [90,] 0.439503709203564 0.00222452690840 -4.53585193982e-06 5.3
     [91,] 0.439503709203564 0.00224964327069 -3.02331937263e-10 9.5
     [92,] 0.450804624124430 0.00216770324934 -4.59455036239e-08 7.3
     >
     > if(requireNamespace("scatterplot3d")) {
     + scatterplot3d::scatterplot3d(res[,1:3], type ='h') ## quite interesting:
     + ## the inaccurate (p,df) points are on nice monotone curve !!!
     + ## this is *less* revealing
     + scatterplot3d::scatterplot3d(res[,c("p","df","nDig")], type ='h')
     + }
     Loading required namespace: scatterplot3d
     > rL <- res[abs(res[,'rE']) > 1e-9,]
     > rL <- rL[order(rL[,1],rL[,2]),]
     > rL
     p df rE nDig
     [1,] 0.0001943754 0.023340796 6.465930e-08 7.2
     [2,] 0.0010123161 0.018556151 -2.591456e-04 3.6
     [3,] 0.0016823889 0.017207366 5.530890e-04 3.3
     [4,] 0.0017467874 0.017311895 -5.868978e-08 7.2
     [5,] 0.0026644512 0.015993174 1.484213e-04 3.8
     [6,] 0.0026644512 0.016180242 7.755647e-08 7.1
     [7,] 0.0031594219 0.015576128 -7.921170e-06 5.1
     [8,] 0.0031594219 0.015681837 -4.522375e-08 7.3
     [9,] 0.0040554624 0.014938587 -7.174047e-06 5.1
     [10,] 0.0044006941 0.014591017 9.079070e-04 3.0
     [11,] 0.0044588113 0.014575069 9.031400e-05 4.0
     [12,] 0.0044818822 0.014688831 -3.233095e-07 6.5
     [13,] 0.0049396099 0.014401684 -2.818101e-06 5.6
     [14,] 0.0088244651 0.012763527 1.211073e-04 3.9
     [15,] 0.0090402660 0.012737116 1.389644e-05 4.9
     [16,] 0.0116421249 0.012014713 -3.002713e-04 3.5
     [17,] 0.0154992134 0.011254201 1.004924e-04 4.0
     [18,] 0.0154992134 0.011359204 -9.557390e-08 7.0
     [19,] 0.0186030166 0.010717161 -2.075379e-03 2.7
     [20,] 0.0186030166 0.010736555 2.143888e-04 3.7
     [21,] 0.0226242424 0.010333795 -3.378658e-09 8.5
     [22,] 0.0226242424 0.010342061 4.430797e-08 7.4
     [23,] 0.0237302174 0.010162521 -1.077327e-06 6.0
     [24,] 0.0447535254 0.008396264 1.222242e-05 4.9
     [25,] 0.0818184250 0.006860077 -1.786334e-09 8.7
     [26,] 0.0828003091 0.006812347 4.179976e-09 8.4
     [27,] 0.0908216581 0.006552698 -7.160336e-09 8.1
     [28,] 0.1022947605 0.006235631 -6.638898e-09 8.2
     [29,] 0.1274058577 0.005623691 6.605415e-09 8.2
     [30,] 0.1352296342 0.005400734 -2.347626e-05 4.6
     [31,] 0.1377322800 0.005330921 2.992858e-04 3.5
     [32,] 0.1381129175 0.005351387 -2.053358e-06 5.7
     [33,] 0.1411006360 0.005273058 4.315938e-05 4.4
     [34,] 0.1429052994 0.005236800 3.481808e-04 3.5
     [35,] 0.1456245572 0.005269240 -1.945018e-09 8.7
     [36,] 0.1546068729 0.005068069 -4.599247e-07 6.3
     [37,] 0.1546068729 0.005073662 2.723010e-07 6.6
     [38,] 0.1697410365 0.004841818 5.314010e-09 8.3
     [39,] 0.1773275763 0.004659561 5.534044e-05 4.3
     [40,] 0.1900940174 0.004503736 1.785654e-06 5.7
     [41,] 0.1901476415 0.004534687 5.662356e-09 8.2
     [42,] 0.2147809514 0.004198693 -1.918794e-08 7.7
     [43,] 0.2328051066 0.003953993 -9.175810e-08 7.0
     [44,] 0.2520765119 0.003749038 2.145699e-07 6.7
     [45,] 0.2530829140 0.003752594 3.654361e-09 8.4
     [46,] 0.2539220587 0.003712373 -5.090149e-06 5.3
     [47,] 0.2542892786 0.003743439 -1.056649e-09 9.0
     [48,] 0.2600174995 0.003661796 -2.904302e-07 6.5
     [49,] 0.2703239068 0.003519992 -1.561648e-04 3.8
     [50,] 0.2716993561 0.003550681 5.130930e-09 8.3
     [51,] 0.2755161961 0.003468040 -4.351715e-04 3.4
     [52,] 0.2901885431 0.003366135 -5.644431e-08 7.2
     [53,] 0.2905790220 0.003344235 1.026676e-07 7.0
     [54,] 0.2905790220 0.003367649 2.260616e-08 7.6
     [55,] 0.2918501987 0.003335528 -1.273388e-06 5.9
     [56,] 0.2965211365 0.003303089 -4.489274e-07 6.3
     [57,] 0.2980341749 0.003304623 8.424704e-09 8.1
     [58,] 0.3005567833 0.003239225 4.660033e-05 4.3
     [59,] 0.3223198463 0.003061345 -1.151308e-05 4.9
     [60,] 0.3250712721 0.003023433 -2.470887e-04 3.6
     [61,] 0.3250712721 0.003041464 -7.713766e-06 5.1
     [62,] 0.3318884124 0.003008371 -4.960989e-09 8.3
     [63,] 0.3853894768 0.002609817 -1.774302e-09 8.8
     [64,] 0.4253339570 0.002329958 1.828230e-08 7.7
     [65,] 0.4395037092 0.002224527 -4.535852e-06 5.3
     [66,] 0.4508046241 0.002167703 -4.594550e-08 7.3
     > plot(rL[,1:2], type = "b", main = "inaccurate pchisq/qchisq pairs")
     >
     > plot(rL[,1:2], type = "b", log = "x", ylim = range(0, rL[,"df"]),
     + xaxt = "n",
     + main = "inaccurate pchisq/qchisq pairs"); abline(h = 0, lty=2)
     > ## aha -- a perfect line !!
     > lines(res[,1:2], col = adjustcolor(1, 0.5))
     > eaxis(1); axis(1, at = 1/2)
     >
     > d <- as.data.frame(res)
     > plot (df ~ log(p), data = d, type = "b", cex=1/4, col="gray")
     > points(df ~ log(p), data = as.data.frame(rL), col=2, cex = 1/2)
     >
     > summary(fm <- lm (df ~ log(p), data = d, weights = -log(abs(rE))))
    
     Call:
     lm(formula = df ~ log(p), data = d, weights = -log(abs(rE)))
    
     Weighted Residuals:
     Min 1Q Median 3Q Max
     -6.912e-04 -1.441e-04 -1.969e-05 7.710e-05 1.082e-03
    
     Coefficients:
     Estimate Std. Error t value Pr(>|t|)
     (Intercept) 6.184e-06 1.133e-05 0.546 0.587
     log(p) -2.725e-03 3.660e-06 -744.530 <2e-16 ***
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     Residual standard error: 0.0002557 on 90 degrees of freedom
     Multiple R-squared: 0.9998, Adjusted R-squared: 0.9998
     F-statistic: 5.543e+05 on 1 and 90 DF, p-value: < 2.2e-16
    
     > ## R^2 = 0.9998
     >
     > p0 <- 2^seq(-50,-1, by=1/8)
     > dN <- data.frame(p = p0,
     + df = predict(fm, newdata = data.frame(p = p0)))
     > rE <- with(dN, 1- pchisq(qchisq(p, df),df)/p)
     > dN <- cbind(dN, rE = rE, nDig = round(-log10(abs(rE)), 1))
     > print(dN, digits=10)
     p df rE nDig
     1 8.881784197e-16 0.094448579669 -5.693322411e-07 6.2
     2 9.685654347e-16 0.094212473681 7.322734105e-07 6.1
     3 1.056228096e-15 0.093976367693 -5.153676885e-08 7.3
     4 1.151824906e-15 0.093740261705 -8.032222418e-07 6.1
     5 1.256073967e-15 0.093504155717 5.733035913e-07 6.2
     6 1.369758374e-15 0.093268049728 -1.195090105e-07 6.9
     7 1.493732098e-15 0.093031943740 -7.801879547e-07 6.1
     8 1.628926404e-15 0.092795837752 6.712882394e-07 6.2
     9 1.776356839e-15 0.092559731764 6.950104159e-08 7.2
     10 1.937130869e-15 0.092323625776 -5.001435726e-07 6.3
     11 2.112456192e-15 0.092087519788 1.026313132e-06 6.0
     12 2.303649813e-15 0.091851413800 -1.542993703e-06 5.8
     13 2.512147934e-15 0.091615307811 3.699662154e-08 7.4
     14 2.739516748e-15 0.091379201823 -4.094315547e-07 6.4
     15 2.987464197e-15 0.091143095835 -8.237025693e-07 6.1
     16 3.257852808e-15 0.090906989847 8.313182130e-07 6.1
     17 3.552713679e-15 0.090670883859 4.759887695e-07 6.3
     18 3.874261739e-15 0.090434777871 1.528253788e-07 6.8
     19 4.224912384e-15 0.090198671883 -1.381691319e-07 6.9
     20 4.607299625e-15 0.089962565895 -3.969919471e-07 6.4
     21 5.024295868e-15 0.089726459906 1.386675560e-06 5.9
     22 5.479033495e-15 0.089490353918 -8.181112241e-07 6.1
     23 5.974928394e-15 0.089254247930 1.019155402e-06 6.0
     24 6.515705616e-15 0.089018141942 -1.110509928e-06 6.0
     25 7.105427358e-15 0.088782035954 7.803687333e-07 6.1
     26 7.748523477e-15 0.088545929966 -1.274165560e-06 5.9
     27 8.449824769e-15 0.088309823978 6.703380386e-07 6.2
     28 9.214599250e-15 0.088073717989 -1.309055653e-06 5.9
     29 1.004859174e-14 0.087837612001 6.890857738e-07 6.2
     30 1.095806699e-14 0.087601506013 -1.215157786e-06 5.9
     31 1.194985679e-14 0.087365400025 8.366343469e-07 6.1
     32 1.303141123e-14 0.087129294037 -9.924495552e-07 6.0
     33 1.421085472e-14 0.086893188049 1.113006151e-06 6.0
     34 1.549704695e-14 0.086657082061 -6.409085871e-07 6.2
     35 1.689964954e-14 0.086420976073 -4.168188492e-07 6.4
     36 1.842919850e-14 0.086184870084 -1.605125546e-07 6.8
     37 2.009718347e-14 0.085948764096 1.280130847e-07 6.9
     38 2.191613398e-14 0.085712658108 4.487608530e-07 6.3
     39 2.389971358e-14 0.085476552120 -1.111735280e-06 6.0
     40 2.606282246e-14 0.085240446132 1.186933915e-06 5.9
     41 2.842170943e-14 0.085004340144 -2.983612493e-07 6.5
     42 3.099409391e-14 0.084768234156 1.566736150e-07 6.8
     43 3.379929908e-14 0.084532128167 6.439440862e-07 6.2
     44 3.685839700e-14 0.084296022179 -7.230811827e-07 6.1
     45 4.019436694e-14 0.084059916191 -1.659619306e-07 6.8
     46 4.383226796e-14 0.083823810203 4.234008103e-07 6.4
     47 4.779942715e-14 0.083587704215 -8.253375789e-07 6.1
     48 5.212564492e-14 0.083351598227 -1.661117941e-07 6.8
     49 5.684341886e-14 0.083115492239 5.253653413e-07 6.3
     50 6.198818782e-14 0.082879386251 -6.050692618e-07 6.2
     51 6.759859815e-14 0.082643280262 1.562851194e-07 6.8
     52 7.371679400e-14 0.082407174274 9.498987034e-07 6.0
     53 8.038873388e-14 0.082171068286 -6.221528537e-08 7.2
     54 8.766453592e-14 0.081934962298 -1.031256456e-06 6.0
     55 9.559885430e-14 0.081698856310 -1.301203456e-07 6.9
     56 1.042512898e-13 0.081462750322 8.032851107e-07 6.1
     57 1.136868377e-13 0.081226644334 -4.741459603e-08 7.3
     58 1.239763756e-13 0.080990538345 -8.550366157e-07 6.1
     59 1.351971963e-13 0.080754432357 1.859174034e-07 6.7
     60 1.474335880e-13 0.080518326369 1.259150868e-06 5.9
     61 1.607774678e-13 0.080282220381 5.698910264e-07 6.2
     62 1.753290718e-13 0.080046114393 -7.628630283e-08 7.1
     63 1.911977086e-13 0.079810008405 -6.793803771e-07 6.2
     64 2.085025797e-13 0.079573902417 -1.239390442e-06 5.9
     65 2.273736754e-13 0.079337796429 1.671630268e-08 7.8
     66 2.479527513e-13 0.079101690440 1.305117147e-06 5.9
     67 2.703943926e-13 0.078865584452 8.634978415e-07 6.1
     68 2.948671760e-13 0.078629478464 4.649672880e-07 6.3
     69 3.215549355e-13 0.078393372476 1.095262072e-07 7.0
     70 3.506581437e-13 0.078157266488 -2.028246688e-07 6.7
     71 3.823954172e-13 0.077921160500 -4.720846172e-07 6.3
     72 4.170051594e-13 0.077685054512 -6.982529182e-07 6.2
     73 4.547473509e-13 0.077448948523 -8.813288566e-07 6.1
     74 4.959055026e-13 0.077212842535 -1.021311716e-06 6.0
     75 5.407887852e-13 0.076976736547 -1.118200798e-06 6.0
     76 5.897343520e-13 0.076740630559 -1.171995389e-06 5.9
     77 6.431098711e-13 0.076504524571 -1.182694798e-06 5.9
     78 7.013162874e-13 0.076268418583 -1.150298324e-06 5.9
     79 7.647908344e-13 0.076032312595 -1.074805275e-06 6.0
     80 8.340103188e-13 0.075796206607 -9.562149603e-07 6.0
     81 9.094947018e-13 0.075560100618 -7.945266962e-07 6.1
     82 9.918110051e-13 0.075323994630 -5.897397957e-07 6.2
     83 1.081577570e-12 0.075087888642 -3.418535914e-07 6.5
     84 1.179468704e-12 0.074851782654 -5.086740451e-08 7.3
     85 1.286219742e-12 0.074615676666 2.832194348e-07 6.5
     86 1.402632575e-12 0.074379570678 6.604075928e-07 6.2
     87 1.529581669e-12 0.074143464690 1.080697732e-06 6.0
     88 1.668020638e-12 0.073907358701 -1.042647599e-07 7.0
     89 1.818989404e-12 0.073671252713 4.076445904e-07 6.4
     90 1.983622010e-12 0.073435146725 -6.748022854e-07 6.2
     91 2.163155141e-12 0.073199040737 -7.127495261e-08 7.1
     92 2.358937408e-12 0.072962934749 5.753558137e-07 6.2
     93 2.572439484e-12 0.072726828761 -3.560679425e-07 6.4
     94 2.805265149e-12 0.072490722773 3.821799130e-07 6.4
     95 3.059163338e-12 0.072254616785 -4.467414887e-07 6.3
     96 3.336041275e-12 0.072018510796 3.831221452e-07 6.4
     97 3.637978807e-12 0.071782404808 -3.433027738e-07 6.5
     98 3.967244020e-12 0.071546298820 -1.015745335e-06 6.0
     99 4.326310282e-12 0.071310192832 -4.575907031e-08 7.3
     100 4.717874816e-12 0.071074086844 9.673320129e-07 6.0
     101 5.144878969e-12 0.070837980856 -1.131696160e-06 5.9
     102 5.610530299e-12 0.070601874868 -2.159353318e-08 7.7
     103 6.118326675e-12 0.070365768879 1.131613904e-06 5.9
     104 6.672082550e-12 0.070129662891 -7.946347651e-07 6.1
     105 7.275957614e-12 0.069893556903 4.555778093e-07 6.3
     106 7.934488041e-12 0.069657450915 1.985071415e-07 6.7
     107 8.652620563e-12 0.069421344927 -4.602383807e-09 8.3
     108 9.435749632e-12 0.069185238939 -1.537536061e-07 6.8
     109 1.028975794e-11 0.068949132951 -2.489493711e-07 6.6
     110 1.122106060e-11 0.068713026963 -2.901925258e-07 6.5
     111 1.223665335e-11 0.068476920974 -2.774859360e-07 6.6
     112 1.334416510e-11 0.068240814986 -2.108324577e-07 6.7
     113 1.455191523e-11 0.068004708998 -9.023495839e-08 7.0
     114 1.586897608e-11 0.067768603010 8.430368981e-08 7.1
     115 1.730524113e-11 0.067532497022 3.127806092e-07 6.5
     116 1.887149926e-11 0.067296391034 -9.005818178e-07 6.0
     117 2.057951587e-11 0.067060285046 9.315377304e-07 6.0
     118 2.244212120e-11 0.066824179057 -1.630697595e-07 6.8
     119 2.447330670e-11 0.066588073069 2.865757222e-07 6.5
     120 2.668833020e-11 0.066351967081 7.901436261e-07 6.1
     121 2.910383046e-11 0.066115861093 -1.208588474e-07 6.9
     122 3.173795216e-11 0.065879755105 -9.670019019e-07 6.0
     123 3.461048225e-11 0.065643649117 -2.908072072e-07 6.5
     124 3.774299853e-11 0.065407543129 4.392954942e-07 6.4
     125 4.115903175e-11 0.065171437141 -2.233106784e-07 6.7
     126 4.488424239e-11 0.064935331152 -8.210848088e-07 6.1
     127 4.894661340e-11 0.064699225164 8.158930909e-08 7.1
     128 5.337666040e-11 0.064463119176 1.038156913e-06 6.0
     129 5.820766091e-11 0.064227013188 6.238512490e-07 6.2
     130 6.347590433e-11 0.063990907200 2.743501614e-07 6.6
     131 6.922096451e-11 0.063754801212 -1.035416930e-08 8.0
     132 7.548599706e-11 0.063518695224 -2.302695521e-07 6.6
     133 8.231806350e-11 0.063282589235 -3.854038086e-07 6.4
     134 8.976848478e-11 0.063046483247 -4.757647607e-07 6.3
     135 9.789322680e-11 0.062810377259 -5.013602509e-07 6.3
     136 1.067533208e-10 0.062574271271 -4.621981136e-07 6.3
     137 1.164153218e-10 0.062338165283 -3.582861912e-07 6.4
     138 1.269518087e-10 0.062102059295 -1.896323394e-07 6.7
     139 1.384419290e-10 0.061865953307 4.375558904e-08 7.4
     140 1.509719941e-10 0.061629847319 3.418697269e-07 6.5
     141 1.646361270e-10 0.061393741330 7.047022075e-07 6.2
     142 1.795369696e-10 0.061157635342 -2.212484238e-07 6.7
     143 1.957864536e-10 0.060921529354 2.764617620e-07 6.6
     144 2.135066416e-10 0.060685423366 -5.036437378e-07 6.3
     145 2.328306437e-10 0.060449317378 1.289049596e-07 6.9
     146 2.539036173e-10 0.060213211390 8.261288568e-07 6.1
     147 2.768838580e-10 0.059977105402 2.619449679e-07 6.6
     148 3.019439882e-10 0.059740999413 -2.266130656e-07 6.6
     149 3.292722540e-10 0.059504893425 6.754948677e-07 6.2
     150 3.590739391e-10 0.059268787437 -9.769083436e-07 6.0
     151 3.915729072e-10 0.059032681449 6.536886865e-08 7.2
     152 4.270132832e-10 0.058796575461 -1.263299489e-07 6.9
     153 4.656612873e-10 0.058560469473 -2.424715895e-07 6.6
     154 5.078072346e-10 0.058324363485 -2.830702768e-07 6.5
     155 5.537677160e-10 0.058088257497 -2.481402319e-07 6.6
     156 6.038879765e-10 0.057852151508 -1.376956813e-07 6.9
     157 6.585445080e-10 0.057616045520 4.824913535e-08 7.3
     158 7.181478783e-10 0.057379939532 3.096799803e-07 6.5
     159 7.831458144e-10 0.057143833544 6.465826115e-07 6.2
     160 8.540265665e-10 0.056907727556 -1.956638767e-07 6.7
     161 9.313225746e-10 0.056671621568 2.976208309e-07 6.5
     162 1.015614469e-09 0.056435515580 8.663322291e-07 6.1
     163 1.107535432e-09 0.056199409591 2.723392490e-07 6.6
     164 1.207775953e-09 0.055963303603 -2.352547019e-07 6.6
     165 1.317089016e-09 0.055727197615 -6.564715971e-07 6.2
     166 1.436295757e-09 0.055491091627 2.302113088e-07 6.6
     167 1.566291629e-09 0.055254985639 -2.383583730e-08 7.6
     168 1.708053133e-09 0.055018879651 -1.915690806e-07 6.7
     169 1.862645149e-09 0.054782773663 -2.730104096e-07 6.6
     170 2.031228938e-09 0.054546667675 -2.681818159e-07 6.6
     171 2.215070864e-09 0.054310561686 -1.771053018e-07 6.8
     172 2.415551906e-09 0.054074455698 1.971357522e-10 9.7
     173 2.634178032e-09 0.053838349710 2.637034958e-07 6.6
     174 2.872591513e-09 0.053602243722 -5.640842873e-07 6.2
     175 3.132583258e-09 0.053366137734 -1.227401389e-07 6.9
     176 3.416106266e-09 0.053130031746 4.047391041e-07 6.4
     177 3.725290298e-09 0.052893925758 -1.426154970e-07 6.8
     178 4.062457877e-09 0.052657819769 -5.928550000e-07 6.2
     179 4.430141728e-09 0.052421713781 2.038662528e-07 6.7
     180 4.831103812e-09 0.052185607793 -5.776660905e-08 7.2
     181 5.268356064e-09 0.051949501805 -2.223745084e-07 6.7
     182 5.745183026e-09 0.051713395817 8.433560312e-07 6.1
     183 6.265166516e-09 0.051477289829 -2.606399732e-07 6.6
     184 6.832212532e-09 0.051241183841 -1.343598028e-07 6.9
     185 7.450580597e-09 0.051005077853 8.882078772e-08 7.1
     186 8.124915754e-09 0.050768971864 -7.023640207e-07 6.2
     187 8.860283457e-09 0.050532865876 8.257587175e-07 6.1
     188 9.662207623e-09 0.050296759888 2.392302341e-07 6.6
     189 1.053671213e-08 0.050060653900 -2.394743450e-07 6.6
     190 1.149036605e-08 0.049824547912 -6.103966785e-07 6.2
     191 1.253033303e-08 0.049588441924 2.100143996e-07 6.7
     192 1.366442506e-08 0.049352335936 4.899553085e-08 7.3
     193 1.490116119e-08 0.049116229947 -4.362272321e-09 8.4
     194 1.624983151e-08 0.048880123959 4.989935698e-08 7.3
     195 1.772056691e-08 0.048644017971 2.117388018e-07 6.7
     196 1.932441525e-08 0.048407911983 -5.748355087e-07 6.2
     197 2.107342426e-08 0.048171805995 -1.924479196e-07 6.7
     198 2.298073210e-08 0.047935700007 2.973889422e-07 6.5
     199 2.506066606e-08 0.047699594019 -1.447314699e-07 6.8
     200 2.732885013e-08 0.047463488031 -4.684300434e-07 6.3
     201 2.980232239e-08 0.047227382042 3.545092193e-07 6.5
     202 3.249966302e-08 0.046991276054 -7.607756496e-07 6.1
     203 3.544113383e-08 0.046755170066 2.876612886e-07 6.5
     204 3.864883049e-08 0.046519064078 -5.800756444e-07 6.2
     205 4.214684851e-08 0.046282958090 6.936620679e-07 6.2
     206 4.596146421e-08 0.046046852102 7.324251272e-08 7.1
     207 5.012133212e-08 0.045810746114 -4.180421775e-07 6.4
     208 5.465770025e-08 0.045574640125 2.092254814e-07 6.7
     209 5.960464478e-08 0.045338534137 -2.954279021e-08 7.5
     210 6.499932603e-08 0.045102428149 -1.393715623e-07 6.9
     211 7.088226765e-08 0.044866322161 -1.203274604e-07 6.9
     212 7.729766099e-08 0.044630216173 2.752292738e-08 7.6
     213 8.429369702e-08 0.044394110185 -6.576521183e-07 6.2
     214 9.192292842e-08 0.044158004197 -2.468590250e-07 6.6
     215 1.002426642e-07 0.043921898209 2.925361181e-07 6.5
     216 1.093154005e-07 0.043685792220 1.531700988e-08 7.8
     217 1.192092896e-07 0.043449686232 -1.223621087e-07 6.9
     218 1.299986521e-07 0.043213580244 -1.205823330e-07 6.9
     219 1.417645353e-07 0.042977474256 2.057528392e-08 7.7
     220 1.545953220e-07 0.042741368268 -6.218908339e-07 6.2
     221 1.685873940e-07 0.042505262280 -1.966860006e-07 6.7
     222 1.838458568e-07 0.042269156292 3.676487106e-07 6.4
     223 2.004853285e-07 0.042033050303 1.647381379e-07 6.8
     224 2.186308010e-07 0.041796944315 1.118715082e-07 7.0
     225 2.384185791e-07 0.041560838327 2.089520004e-07 6.7
     226 2.599973041e-07 0.041324732339 4.558828597e-07 6.3
     227 2.835290706e-07 0.041088626351 -3.149587724e-08 7.5
     228 3.091906439e-07 0.040852520363 -3.581121097e-07 6.4
     229 3.371747881e-07 0.040616414375 3.488470744e-07 6.5
     230 3.676917137e-07 0.040380308387 -5.295130916e-07 6.3
     231 4.009706570e-07 0.040144202398 4.872893066e-07 6.3
     232 4.372616020e-07 0.039908096410 -5.923112045e-08 7.2
     233 4.768371582e-07 0.039671990422 -4.344346312e-07 6.4
     234 5.199946082e-07 0.039435884434 2.066676108e-07 6.7
     235 5.670581412e-07 0.039199778446 1.681372633e-07 6.8
     236 6.183812879e-07 0.038963672458 -5.334671598e-07 6.3
     237 6.743495762e-07 0.038727566470 -2.247258435e-07 6.6
     238 7.353834273e-07 0.038491460481 2.546716327e-07 6.6
     239 8.019413140e-07 0.038255354493 8.725913336e-08 7.1
     240 8.745232040e-07 0.038019248505 1.013349561e-07 7.0
     241 9.536743164e-07 0.037783142517 -5.094306235e-07 6.3
     242 1.039989216e-06 0.037547036529 -1.272824359e-07 6.9
     243 1.134116282e-06 0.037310930541 4.358925881e-07 6.4
     244 1.236762576e-06 0.037074824553 3.904291315e-07 6.4
     245 1.348699152e-06 0.036838718565 5.367821070e-07 6.3
     246 1.470766855e-06 0.036602612576 9.641457999e-08 7.0
     247 1.603882628e-06 0.036366506588 -1.413544166e-07 6.8
     248 1.749046408e-06 0.036130400600 -1.767197404e-07 6.8
     249 1.907348633e-06 0.035894294612 -9.875985141e-09 8.0
     250 2.079978433e-06 0.035658188624 -3.970802771e-07 6.4
     251 2.268232565e-06 0.035422082636 1.791465434e-07 6.7
     252 2.473525152e-06 0.035185976648 -5.328653061e-07 6.3
     253 2.697398305e-06 0.034949870659 -2.818818945e-07 6.5
     254 2.941533709e-06 0.034713764671 1.813954295e-07 6.7
     255 3.207765256e-06 0.034477658683 1.285281512e-07 6.9
     256 3.498092816e-06 0.034241552695 2.986249925e-07 6.5
     257 3.814697266e-06 0.034005446707 -2.562239221e-08 7.6
     258 4.159956866e-06 0.033769340719 -1.162612726e-07 6.9
     259 4.536465130e-06 0.033533234731 2.644014629e-08 7.6
     260 4.947050303e-06 0.033297128743 4.022141099e-07 6.4
     261 5.394796609e-06 0.033061022754 -3.787335339e-07 6.4
     262 5.883067419e-06 0.032824916766 4.734672679e-07 6.3
     263 6.415530512e-06 0.032588810778 1.906558875e-07 6.7
     264 6.996185632e-06 0.032352704790 1.620354394e-07 6.8
     265 7.629394531e-06 0.032116598802 3.872828165e-07 6.4
     266 8.319913732e-06 0.031880492814 -4.676686218e-07 6.3
     267 9.072930260e-06 0.031644386826 2.754435448e-07 6.6
     268 9.894100606e-06 0.031408280837 -4.002658072e-08 7.4
     269 1.078959322e-05 0.031172174849 -8.071856361e-08 7.1
     270 1.176613484e-05 0.030936068861 1.529846981e-07 6.8
     271 1.283106102e-05 0.030699962873 2.171783053e-08 7.7
     272 1.399237126e-05 0.030463856885 1.752040166e-07 6.8
     273 1.525878906e-05 0.030227750897 -1.479256340e-08 7.8
     274 1.663982746e-05 0.029991644909 9.026352044e-08 7.0
     275 1.814586052e-05 0.029755538921 -1.267311944e-07 6.9
     276 1.978820121e-05 0.029519432932 -3.843675112e-08 7.4
     277 2.157918644e-05 0.029283326944 -2.508241883e-07 6.6
     278 2.353226967e-05 0.029047220956 -1.477539424e-07 6.8
     279 2.566212205e-05 0.028811114968 2.702582267e-07 6.6
     280 2.798474253e-05 0.028575008980 -1.748270420e-07 6.8
     281 3.051757812e-05 0.028338902992 -2.837500699e-07 6.5
     282 3.327965493e-05 0.028102797004 -5.710114848e-08 7.2
     283 3.629172104e-05 0.027866691015 -6.748044701e-08 7.2
     284 3.957640242e-05 0.027630585027 2.676487311e-07 6.6
     285 4.315837288e-05 0.027394479039 3.867922076e-07 6.4
     286 4.706453935e-05 0.027158373051 -2.492501328e-07 6.6
     287 5.132424410e-05 0.026922267063 4.134786513e-08 7.4
     288 5.596948506e-05 0.026686161075 -3.911749389e-07 6.4
     289 6.103515625e-05 0.026450055087 1.014410593e-07 7.0
     290 6.655930986e-05 0.026213949098 -9.713157123e-08 7.0
     291 7.258344208e-05 0.025977843110 1.004778277e-07 7.0
     292 7.915280485e-05 0.025741737122 -3.501362080e-07 6.5
     293 8.631674575e-05 0.025505631134 -3.839952303e-07 6.4
     294 9.412907870e-05 0.025269525146 -2.027959578e-09 8.7
     295 1.026484882e-04 0.025033419158 -2.152639500e-07 6.7
     296 1.119389701e-04 0.024797313170 7.678417235e-09 8.1
     297 1.220703125e-04 0.024561207182 -3.220049949e-07 6.5
     298 1.331186197e-04 0.024325101193 -1.953370445e-07 6.7
     299 1.451668842e-04 0.024088995205 -9.617589081e-08 7.0
     300 1.583056097e-04 0.023852889217 -8.977143029e-09 8.0
     301 1.726334915e-04 0.023616783229 8.175071808e-08 7.1
     302 1.882581574e-04 0.023380677241 1.914434747e-07 6.7
     303 2.052969764e-04 0.023144571253 -1.249542994e-07 6.9
     304 2.238779402e-04 0.022908465265 7.430447435e-08 7.1
     305 2.441406250e-04 0.022672359276 -1.108169194e-07 7.0
     306 2.662372394e-04 0.022436253288 2.389977569e-07 6.6
     307 2.903337683e-04 0.022200147300 2.458562984e-07 6.6
     308 3.166112194e-04 0.021964041312 -5.850366769e-08 7.2
     309 3.452669830e-04 0.021727935324 2.115776949e-07 6.7
     310 3.765163148e-04 0.021491829336 2.113798955e-07 6.7
     311 4.105939528e-04 0.021255723348 -2.762990770e-08 7.6
     312 4.477558805e-04 0.021019617360 -6.378203188e-08 7.2
     313 4.882812500e-04 0.020783511371 -2.872497766e-07 6.5
     314 5.324744788e-04 0.020547405383 -2.677229627e-07 6.6
     315 5.806675366e-04 0.020311299395 9.137057333e-09 8.0
     316 6.332224388e-04 0.020075193407 1.695002406e-07 6.8
     317 6.905339660e-04 0.019839087419 -1.383746775e-07 6.9
     318 7.530326296e-04 0.019602981431 2.636260105e-07 6.6
     319 8.211879055e-04 0.019366875443 -1.129219587e-07 6.9
     320 8.955117609e-04 0.019130769454 2.589505672e-07 6.6
     321 9.765625000e-04 0.018894663466 -6.508972539e-08 7.2
     322 1.064948958e-03 0.018658557478 4.241218476e-08 7.4
     323 1.161335073e-03 0.018422451490 2.453652983e-07 6.6
     324 1.266444878e-03 0.018186345502 2.296493848e-07 6.6
     325 1.381067932e-03 0.017950239514 4.120031982e-08 7.4
     326 1.506065259e-03 0.017714133526 5.827584493e-08 7.2
     327 1.642375811e-03 0.017478027538 -1.735446653e-08 7.8
     328 1.791023522e-03 0.017241921549 1.809493316e-07 6.7
     329 1.953125000e-03 0.017005815561 4.947543719e-08 7.3
     330 2.129897915e-03 0.016769709573 -4.023365396e-08 7.4
     331 2.322670146e-03 0.016533603585 -4.402659903e-08 7.4
     332 2.532889755e-03 0.016297497597 8.187645184e-08 7.1
     333 2.762135864e-03 0.016061391609 -2.069273115e-07 6.7
     334 3.012130518e-03 0.015825285621 3.071764043e-08 7.5
     335 3.284751622e-03 0.015589179632 -2.785656417e-08 7.6
     336 3.582047044e-03 0.015353073644 -3.027688034e-08 7.5
     337 3.906250000e-03 0.015116967656 -1.904897100e-07 6.7
     338 4.259795831e-03 0.014880861668 -1.683114372e-07 6.8
     339 4.645340293e-03 0.014644755680 9.292967829e-08 7.0
     340 5.065779510e-03 0.014408649692 1.383866128e-07 6.9
     341 5.524271728e-03 0.014172543704 5.614294096e-08 7.3
     342 6.024261037e-03 0.013936437716 -6.658807195e-08 7.2
     343 6.569503244e-03 0.013700331727 -1.435434784e-07 6.8
     344 7.164094088e-03 0.013464225739 -8.948700758e-08 7.0
     345 7.812500000e-03 0.013228119751 -4.838022560e-08 7.3
     346 8.519591661e-03 0.012992013763 -1.436152919e-07 6.8
     347 9.290680586e-03 0.012755907775 -4.376058627e-08 7.4
     348 1.013155902e-02 0.012519801787 1.358225326e-07 6.9
     349 1.104854346e-02 0.012283695799 9.346486218e-08 7.0
     350 1.204852207e-02 0.012047589810 -2.908217578e-08 7.5
     351 1.313900649e-02 0.011811483822 -9.213229890e-08 7.0
     352 1.432818818e-02 0.011575377834 4.174524038e-08 7.4
     353 1.562500000e-02 0.011339271846 1.379139563e-07 6.9
     354 1.703918332e-02 0.011103165858 1.954543860e-09 8.7
     355 1.858136117e-02 0.010867059870 1.447486908e-09 8.8
     356 2.026311804e-02 0.010630953882 -2.719908920e-08 7.6
     357 2.209708691e-02 0.010394847894 1.180522117e-07 6.9
     358 2.409704415e-02 0.010158741905 2.375664332e-09 8.6
     359 2.627801298e-02 0.009922635917 3.493183931e-08 7.5
     360 2.865637635e-02 0.009686529929 7.828329762e-09 8.1
     361 3.125000000e-02 0.009450423941 5.459876207e-08 7.3
     362 3.407836665e-02 0.009214317953 4.834763923e-08 7.3
     363 3.716272234e-02 0.008978211965 -8.341425239e-08 7.1
     364 4.052623608e-02 0.008742105977 1.954604523e-08 7.7
     365 4.419417382e-02 0.008505999988 -2.239338204e-08 7.6
     366 4.819408829e-02 0.008269894000 -1.243319980e-08 7.9
     367 5.255602595e-02 0.008033788012 4.993745217e-08 7.3
     368 5.731275270e-02 0.007797682024 1.741400035e-08 7.8
     369 6.250000000e-02 0.007561576036 4.626561989e-08 7.3
     370 6.815673329e-02 0.007325470048 -3.471772847e-08 7.5
     371 7.432544469e-02 0.007089364060 8.226451520e-09 8.1
     372 8.105247217e-02 0.006853258072 -4.223112349e-08 7.4
     373 8.838834765e-02 0.006617152083 1.657740989e-08 7.8
     374 9.638817659e-02 0.006381046095 3.659368508e-08 7.4
     375 1.051120519e-01 0.006144940107 9.755834918e-09 8.0
     376 1.146255054e-01 0.005908834119 8.519715378e-09 8.1
     377 1.250000000e-01 0.005672728131 2.632288454e-08 7.6
     378 1.363134666e-01 0.005436622143 1.039711917e-08 8.0
     379 1.486508894e-01 0.005200516155 -3.554264660e-08 7.4
     380 1.621049443e-01 0.004964410166 -1.102801761e-08 8.0
     381 1.767766953e-01 0.004728304178 -2.268685084e-08 7.6
     382 1.927763532e-01 0.004492198190 -1.758708290e-09 8.8
     383 2.102241038e-01 0.004256092202 -1.118194382e-08 8.0
     384 2.292510108e-01 0.004019986214 -1.232663527e-09 8.9
     385 2.500000000e-01 0.003783880226 -3.901555301e-09 8.4
     386 2.726269332e-01 0.003547774238 2.407146371e-09 8.6
     387 2.973017788e-01 0.003311668250 3.332240928e-09 8.5
     388 3.242098887e-01 0.003075562261 1.847989872e-09 8.7
     389 3.535533906e-01 0.002839456273 5.534314118e-10 9.3
     390 3.855527064e-01 0.002603350285 -4.492217709e-09 8.3
     391 4.204482076e-01 0.002367244297 1.137717809e-09 8.9
     392 4.585020216e-01 0.002131138309 -5.251989954e-10 9.3
     393 5.000000000e-01 0.001895032321 -1.494597335e-09 8.8
     >
     > ## } ## only when we find inaccurate regions
     > showProc.time()
     Time elapsed: 0.11 0.01 0.17
     >
     >
     > ## Oops: another qgamma() / qchisq() problem: mostly NaN's == all solved now
     > curve(qgamma(x, 20), 1e-16, 1e-10, log='x')
     > curve(qgamma(x, 20), 1e-300, .99 , log='xy') # and add the critical region from above:
     > abline(v=c(1e-16, 1e-10), col="light blue")
     > curve(qgamma(x, 20), 1e-26, 1e-07, log='x')
     > ##-> now using log=TRUE in same region:
     > curve(qgamma(x, 20, log=TRUE), -38, -16)## no problem!!
     > curve(qgamma(exp(x), 20), add=TRUE, col="green3", n=2001)
     > ## had problem here, but no longer !
     >
     > ##--> Further fix for qgamma: when 'x' is very small: use "log=TRUE of log(x)"!
     >
     > ## had bug (gave NaN), but no longer:
     > (q_12 <- qgamma(1e-12, 20))
     [1] 2.330042
     > all.equal(1e-12, pgamma(q_12, 20), tol=0)# show rel.err (Lnx 64-bit: 4.04e-16)
     [1] "Mean relative difference: 4.038968e-16"
     > stopifnot(
     + all.equal(1e-12, pgamma(q_12, 20), tolerance = 1e-14)
     + )
     >
     >
     > ## --- Nice graphic : --- but amazingly *S..L..O..W*
     >
     > p.qgammaSml <- function(from= 1e-110, to = 1e-5, ylim = c(0.4, 1000),
     + n = 201, k.lab = 3,
     + a1 = c(10, seq(10.1,20, by=.2), 21:105),
     + a2 = seq(110,330, by=10),
     + a3 = seq(350,1600, by=50))
     + {
     + ## Purpose: nice qgamma() lines ``for small x'' aka p
     + ## ----------------------------------------------------------------------
     + ## Author: Martin Maechler, Date: 22 Mar 2004, 14:23
     + x <- exp(seq(log(from), log(to), length = n))
     +
     + op <- par(las=1, lab = c(10,10, 7), xaxs = "i", mex = 0.8)
     + on.exit(par(op))
     + plot(x, qgamma(x, a1[1]), log="xy", ylim=ylim, type='l', xaxt = "n",
     + main = paste("qgamma(x, a) for very small x, a in [",
     + formatC(a1[1]),", ",formatC(max(a1,a2,a3)),"] - log-log", sep=''),
     + sub = R.version.string)
     + lab.x <- pretty(log10(c(from,to)), 20)
     + axis(1, at=10^lab.x, lab = paste("10^",formatC(lab.x),sep=''))
     + if(is.nan(qgamma(1e-12, 20)))
     + text(1e-60, 20, "all NaN", cex = 2)
     + if(!is.finite(qgamma(1e-140, 155)))
     + text(1e-240, 5, "all +Inf", cex = 2)
     +
     + lines.txt <- function(a.s, col = par("col")) {
     + col <- rep(col, length=length(a.s))
     + for(i in seq(along=a.s)) {
     + qx <- qgamma(x, (a <- a.s[i]))
     + if(i %% k.lab == 0 &&
     + any(ifi <- is.finite(qx) & qx >= ylim[1])) {
     + ik <- (i%%(2*k.lab))/k.lab # = 0 or 1
     + j <- quantile(which(ifi), c(.02,(1:3)/4+ ik/10, .98))
     + ## "segments" around the labels :
     + i0 <- 1
     + for(jj in j) {
     + ii <- i0:(jj-1)
     + i2 <- jj + -1:1
     + lines(x[ii], qx[ii], col=col[i])
     + lines(x[i2], qx[i2], col=col[i], type = 'c')
     + i0 <- jj+1
     + }
     + text(x[j], qx[j], formatC(a), col= "gray40", cex = 0.8)
     + }
     + else
     + lines(x, qx, col=col[i])
     +
     + }
     + }
     + oo <- options(warn = -1)
     + lines.txt(a1[-1])
     + lines.txt(a2, col= 2)
     + lines.txt(a3, col= rainbow(length(a3), .8, .8,
     + start = (max(a3)-min(a3))/(1+max(a3))))
     + invisible(options(oo))
     + }
     >
     > showProc.time()
     Time elapsed: 0.03 0 0.05
     >
     > p.qgammaSml()
     > p.qgammaSml(1e-300)
     > p.qgammaSml(1e-300,1e-50, a2= seq(100,360, by=4), a3=seq(350,1500, by=10))
     >
     > showProc.time()
     Time elapsed: 1.54 0.02 1.59
     >
     > ## The "upper" problematic corner:
     > p.qgammaSml(1e-19, 1e-3, a2=NULL,a3=NULL, ylim=c(.1,20))
     > p.qgammaSml(1e-19, 1e-3, a2=seq(1,12, by=.04), ylim=c(.1,20),a3=NULL,k.lab=10)
     > ## now shows the problem (quite well):
     > ## could it be in pgamma()'s inaccuracy, leading to qgamma() bias ?
     > aa <- c(seq(9, 22, by=0.25),seq(22.3,40,by=0.4))
     > caa <- formatC(range(aa))
     > sfsmisc::mult.fig(2)
     > curve(pgamma(x, sh=aa[1]), 0.5, 20, log = 'xy', ylim = c(1e-60, .2),
     + main = sprintf("pgamma(x, a) for a in [%s,%s]", caa[1],caa[2]))
     > for(sh in aa) curve(pgamma(x, sh), add = TRUE, col=2)
     > abline(h=c(1e-15), col="light blue", lty=2)
     >
     > curve(pgamma(x, sh=aa[1]), 0.5, 20, log = 'xy', ylim = c(1e-15, .8),
     + main = sprintf("pgamma(x, a) for a in [%s,%s]", caa[1],caa[2]))
     > for(sh in aa) curve(pgamma(x, sh), add = TRUE, col=2)
     > ## the "border curve" between "Pearson" and "Continued fraction (upper tail)"
     > ## in pgamma.c :
     > curve(pgamma(max(1,x), x), add = TRUE, col=4)
     > ## ==> pgamma() is perfect here {series expansion up to eps_C accuracy}!
     >
     > aa <- c(seq(9, 22, by=0.25),seq(22.3,40.4,by=0.4))
     > p.qgammaSml(1e-24, 1e-5, a1=aa, a2=NULL,a3=NULL, ylim=c(.8,8))
     > ## -------- save the above?
     > aa1 <- c(aa,seq(40.5,90, by=0.5))
     > p.qgammaSml(1e-60, 1e-5, a1=aa1, a2=NULL,a3=NULL, ylim=c(.9, 16))
     > aa2 <- c(aa1, seq(91,150, by= 1))
     > p.qgammaSml(1e-90, 1e-5, a1=aa2, a2=NULL,a3=NULL, ylim=c(.9, 35))
     > aa3 <- c(aa2, seq(150,250, by= 2), seq(253, 400, by=5))
     > p.qgammaSml(1e-200, 1e-5, a1=aa3, a2=NULL,a3=NULL, ylim=c(.9, 100))
     > p.qgammaSml(1e-200, 1e-5, a1=aa3, a2=NULL,a3=NULL, ylim=c(.9, 200),k.lab=9e9)
     > p.qgammaSml(1e-60, 1e-5, a1=aa3, a2=NULL,a3=NULL, ylim=c(.9, 200),k.lab=9e9)
     >
     > showProc.time()
     Time elapsed: 4.14 0.05 4.23
     >
     > ## lower a \> 10
     >
     > curve(qgamma(x, 19), 1e-14, 1e-9, log='x')
     > curve(qgamma(x, 18), 1e-14, 1e-9, log='x')
     > curve(qgamma(x, 15), 1e-11, 5e-9, log='x')
     > curve(qgamma(x, 13), 5e-10, 1e-8, log='x')
     > curve(qgamma(x, 11), 1e-8, 5e-8, log='x')
     > curve(qgamma(x, 10.5), 4.2e-8, 6e-8, log='x')
     > curve(qgamma(x, 10.3), 6e-8, 7e-8, log='x')
     > curve(qgamma(x, 10.2), 7.1e-8, 7.6e-8, log='x')
     > curve(qgamma(x, 10.15),7.7e-8, 7.9e-8, log='x')
     > curve(qgamma(x, 10.14),7.88e-8,7.92e-8, log='x',n=10001)
     >
     > ## no more problems for smaller a!! here:
     > curve(qgamma(x, 10.13), 1e-10, 5e-4, log='x',n=20001)
     > curve(qgamma(x, 10.12), 1e-10, 5e-4, log='x',n=20001)
     > curve(qgamma(x, 10.1), 1e-10, 5e-4, log='x',n=20001)
     >
     > showProc.time()
     Time elapsed: 0.73 0.01 0.81
     >
     > ##--- the "+Inf" / premature "0" case:
     > curve(qgamma(x, 155, log=TRUE), -1500, 0, log='y', n=2001,col=2)
     > curve(qgamma(x, 1e3, log=TRUE), -1500, 0, log='y', n=2001,col=2)
     > ## now works, but slowly and with kink
     > curve(qgamma (x, 1e5, log=TRUE), -3e5, 0, log='y', n=2001,col=2,lwd=3)
     > curve(qgammaAppr(x, 1e5, log=TRUE), add = TRUE, n=2001, col="blue",lwd=.4)
     > ## --- curves are almost "identical"
     > ## ===> the kink *does* come from the initial approx... hmm
     >
     > ## still "identical"
     > curve(qgamma (x, 1e4, log=TRUE), -3e4, 0, log='y', n=2001,col=2)
     > curve(qgammaAppr(x, 1e4, log=TRUE), add = TRUE, n=2001, col="tomato3")
     >
     > ## now see some difference (approx. has kink at ~ -165)
     > curve(qgamma (x, 100, log=TRUE), -200, 0, log='y', n=2001,col=2)
     > curve(qgammaAppr(x, 100, log=TRUE), add = TRUE, n=2001, col="tomato3")
     > ##
     > (kk <- 100 * 2/1.24)# 161.29
     [1] 161.2903
     > curve(qgamma (x, 100, log=TRUE), -1.1*kk, -.95*kk, log='y', n=2001,col=2)
     > curve(qgammaAppr(x, 100, log=TRUE), add = TRUE, n=2001, col="tomato3")
     > abline(v = -kk, col='blue', lty=2)# exactly: kink is at a * 2 / 1.24 = a / .62
     > curve(qgammaAppr(x - 100/.62, 100,log=TRUE), -1e-3, +1e-3)
     >
     > showProc.time()
     Time elapsed: 0.2 0 0.2
     >
     > p.qgammaLog <- function(alpha, xl.f = 1.5, xr.f = 0.4, n = 2001)
     + {
     + ## Purpose:
     + ## ----------------------------------------------------------------------
     + ## Arguments:
     + ## ----------------------------------------------------------------------
     + ## Author: Martin Maechler, Date: 30 Mar 2004, 18:44
     + kk <- -alpha / .62 # = (alpha * 2) / (-1.24)
     + curve(qgamma(x, alpha, log=TRUE), xl.f*kk, xr.f*kk, log='y',
     + n=n, col=2, lwd=3.6, lty = 4,
     + main= paste("qgamma(x, alpha=",formatC(alpha,digits=10),", log = TRUE)"))
     + lines(kk, qgamma(kk, alpha, log=TRUE), type = 'h', lty = 3)
     + curve(qgamma (exp(x), alpha), add = TRUE, col="orange", n=n, lwd= 2)
     + curve(qgammaAppr(x, alpha, log=TRUE), add = TRUE, col=3, n=n,lwd = .4)
     + }
     >
     > showProc.time()
     Time elapsed: 0 0 0
     >
     > p.qgammaLog(25)
     > p.qgammaLog(16)# ~ [-25, -20]
     > p.qgammaLog(12, 1.2, 0.8)# small problem remaining
     > p.qgammaLog(11, 1.2, 0.8)# even smaller
     > p.qgammaLog(10.5, 1.1, 0.9)# even smaller
     > p.qgammaLog(10.25, 1.1, 0.9)# even smaller
     > ## 2019-08: __nothing__ visible from here on:
     > p.qgammaLog(10.18, 1.02, 0.98)# even smaller
     > p.qgammaLog(10.15, 1.02, 0.98)# even smaller
     > p.qgammaLog(10.14, 1.001, 0.999)# even smaller
     > p.qgammaLog(10.139, 1.0002, 0.9998)#
     > p.qgammaLog(10.138, 1.0002, 0.9998)#
     > p.qgammaLog(10.137, 1.00001, 0.99999)#
     > p.qgammaLog(10.13699, 1.0000001, 0.9999999)#
     > p.qgammaLog(10.1369899, 1.0000001, 0.9999999)#
     > p.qgammaLog(10.1369894, 1.0000001, 0.9999999)#
     > p.qgammaLog(10.1369893, 1.0000001, 0.9999999)# even smaller at -16.34998
     >
     > showProc.time()
     Time elapsed: 0.6 0.02 0.69
     >
     > ##-- here is the boundary --- for 64-bit AMD Opteron ---
     > ## and for 32-bit AMD Athlon
     >
     > p.qgammaLog(10.1369892, 1.0000001, 0.9999999)# no more
     > p.qgammaLog(10.136989, 1.0000001, 0.9999999)#
     > p.qgammaLog(10.136988, 1.0000001, 0.9999999)#
     > p.qgammaLog(10.136985, 1.0000001, 0.9999999)#
     > p.qgammaLog(10.13698, 1.0000001, 0.9999999)#
     > p.qgammaLog(10.13697, 1.0000001, 0.9999999)#
     > p.qgammaLog(10.13695, 1.0000001, 0.9999999)#
     > p.qgammaLog(10.1368, 1.000001, 0.999999)#
     > p.qgammaLog(10.1365, 1.000001, 0.999999)#
     > p.qgammaLog(10.136, 1.000001, 0.999999)#
     > p.qgammaLog(10.125, 1.1, 0.9)# --- see it now
     > p.qgammaLog(10, 1.2, 0.8)
     > p.qgammaLog(9)
     >
     > showProc.time()
     Time elapsed: 0.51 0.01 0.56
     >
     > ## For large alpha: show difference to see problem better
     > ## ---> for alpha >= 10, the x problem starts *roughly* at x = -0.8*alpha
     > ##
     >
     > sfsmisc::mult.fig(2)
     > curve(qgammaAppr(x, 5, log=TRUE), - 8.1, -8, n=2001)
     > curve(qgammaAppr(x- 5/.62, 5, log=TRUE), -1e-15, 0)
     >
     > ## is the kink from pgamma() ? : no: this looks fine,
     > curve(pgamma(x, 1e5, log=TRUE), 1, 2e5, log='x', n=2001,col=2)
     > ## and this does too:
     > curve( dgamma(x, 1e5), .5e5, 2e5); par(new=TRUE)
     > curve( dgamma(x, 1e5, log=TRUE), .5e5, 2e5, col=2, yaxt="n")
     > axis(4,col.axis=2); par(new=TRUE)
     > curve( pgamma(x, 1e5), .5e5, 2e5, n=2001, col=3); par(new=TRUE)
     > curve( pgamma(x, 1e5, log=TRUE), .5e5, 2e5, n=2001, col=4); par(new=TRUE)
     > curve(-pgamma(x, 1e5, log=TRUE,lower=FALSE), .5e5, 2e5, n=2001, col=4)
     > ## all looking nice
     >
     >
     > x <- 10^seq(2,6, length=4001)
     > qx <- qgamma(pgamma(x, 1e5, log=TRUE), 1e5, log=TRUE)
     > plot(x, qx, type ='l', col=2, asp = 1); abline(0,1, lty=3)
     >
     > showProc.time()
     Time elapsed: 0.08 0 0.09
     > <0c>
     > ###------------- Approximations of qgamma() ------
     > ##
     >
     > ## source("/u/maechler/R/MM/NUMERICS/dpq-functions/qchisqAppr.R")
     > ##--> qchisqAppr()
     > ##--> qchisqWH [ = Wilson Hilferty ]
     > ##--> qchisqKG [ = Kennedy & Gentle's improvements "a la AS 91" ]
     > ## dyn.load("/u/maechler/R/MM/NUMERICS/dpq-functions/qchisq_appr.so")
     >
     > ## Consider the two different implementations of
     > ## lgamma1p(a) := lgamma(1+a) == log(gamma(1+a) == log(a*gamma(a)) "stable":
     >
     > if(!exists("lseq", mode="function"))
     + lseq <- if(requireNamespace("sfsmisc")) sfsmisc::lseq else
     + function(from, to, length) exp(seq(log(from), log(to), length.out = length))
     >
     > if(require("Rmpfr")) { ##---------------- MPFR numbers -------------------------
     +
     + .mpfr.all.eq <- Rmpfr::all.equal
     + AllEq <- function(target, current, ...)
     + .mpfr.all.eq(target, current, ...,
     + formatFUN = function(x, ...) Rmpfr::format(x, digits = 9))
     +
     + print(gammaE <- Const("gamma",200)); pi. <- Const("pi",200)
     + print(a0 <- (gammaE^2 + pi.^2/6)/2)
     + print(psi2.1 <- -2*zeta(mpfr(3,200)))# == psigamma(1,2) =~ -2.4041138
     + print(a1 <- (psi2.1 - gammaE*(pi.^2/2 + gammaE^2))/6)
     +
     + x <- lseq(1e-30, 0.8, length = if(doExtras) 1000 else 125)
     + x. <- mpfr(x, 200)
     + xct. <- log(x. * gamma(x.)) ## using MPFR arithmetic .. no overflow ...
     + xc2. <- log(x.) + lgamma(x.)## (ditto)
     + print(AllEq(xct., xc2., tol = 0)) # 3.15779......e-57
     + xct <- as.numeric(xct.)
     + stopifnot(exprs = {
     + AllEq(xct., xc2., tol = 1e-45)
     + AllEq(xct , xc2., tol = 1e-15)
     + ##
     + all.equal(lgamma1p(x), lgamma1p(x, tol= 1e-16), tol=0)
     + ## -> no difference; i.e., default tol = 1e-14 seems fine enough!
     + })
     + showProc.time()
     +
     + m.appr <- cbind(log(x*gamma(x)), lgamma(1+x), log(x) + lgamma(x),
     + lgamma1p.(x, k=1, cut=3e-6),
     + lgamma1p.(x, k=2, cut=1e-4),
     + lgamma1p.(x, k=3, cut=8e-4),
     + lgamma1p(x))#, tol= 1e-14), # = default
     +
     + eMat <- m.appr - xct # absolute error
     + ## Relative errors:
     + str(reMat. <- m.appr/xct. - 1)
     + str(reMat <- as(reMat., "array")) # as(., "matrix") fails in older versions
     +
     + matplot(x, eMat , log="x", type="l", lty=1) #-> problematic log(x) + lgamma(x) for "large"
     + matplot(x, abs( eMat), log="xy", type="l", lty=1) #-> but good for small; lgamma1p is much better
     + matplot(x, abs(reMat), log="xy", type="l", lty=1)
     + abline(v= 3.47548562941137e-08, col = "gray80", lwd=3)#<- the cutoff value of lgamma1p()
     + ##---> should use earlier cutoff!
     + ## zoom in:
     + matplot(x, abs(reMat), log="xy", type="l", lty=1, xlim=c(8e-9, 1e-3))
     + abline(v= 3.47548562941137e-08, col = "gray80", lwd=3)#<- the cutoff value of lgamma1p()
     +
     + ## rm(x., xct., xc2., reMat., eMat, AllEq)
     + detach("package:Rmpfr")
     + showProc.time()
     +
     + } ## if( MPFR ) ----------------------------------------------------------------
     Loading required package: Rmpfr
     Loading required package: gmp
    
     Attaching package: 'gmp'
    
     The following objects are masked from 'package:base':
    
     %*%, apply, crossprod, matrix, tcrossprod
    
     C code of R package 'Rmpfr': GMP using 64 bits per limb
    
    
     Attaching package: 'Rmpfr'
    
     The following object is masked from 'package:gmp':
    
     outer
    
     The following objects are masked from 'package:stats':
    
     dbinom, dnorm, dpois, pnorm
    
     The following objects are masked from 'package:base':
    
     cbind, pmax, pmin, rbind
    
     1 'mpfr' number of precision 200 bits
     [1] 0.5772156649015328606065120900824024310421593359399235988057672
     1 'mpfr' number of precision 200 bits
     [1] 0.9890559953279725553953956515006347079391835207282140904431957
     1 'mpfr' number of precision 200 bits
     [1] -2.404113806319188570799476323022899981529972584680997763584544
     1 'mpfr' number of precision 200 bits
     [1] -0.9074790760808862890165601673562751149286114490725637609413306
     [1] "Mean relative difference: 3.181053350957963600882729201783933343641695627742022642032500e-57"
     Time elapsed: 1.42 0 1.44
     Class 'mpfrArray' [package "Rmpfr"] of dimension c(125, 7) and precision 200
     -1 -1 ...
     num [1:125, 1:7] -1.00 -1.00 -1.00 3.64e+13 -1.00 ...
     - attr(*, "dimnames")=List of 2
     ..$ : NULL
     ..$ : NULL
     Time elapsed: 0.08 0.02 0.09
     Warning messages:
     1: In xy.coords(x, y, xlabel, ylabel, log = log) :
     348 y values <= 0 omitted from logarithmic plot
     2: In xy.coords(x, y, xlabel, ylabel, log) :
     1 y value <= 0 omitted from logarithmic plot
     >
     >
     > ## ../R/qchisqAppr.R -- talks about the "small shape" qgamma() approxmation
     > ## ----------------- --> .qgammaApprBnd() :
     > curve(.qgammaApprBnd, 1e-18, 1e-15, col=2)
     > abline(h=0, col="gray70", lty=2)
     > eps.c <- .Machine$double.eps
     > axis(3,at=(1:3)* eps.c,
     + label=expression(epsilon[c], 2*epsilon[c], 3*epsilon[c]))
     > (rt.b <- uniroot(.qgammaApprBnd, c(1,3)*eps.c, tol=1e-12))
     $root
     [1] 2.220446e-16
    
     $f.root
     [1] 1.281676e-16
    
     $iter
     [1] 0
    
     $init.it
     [1] NA
    
     $estim.prec
     [1] 4.440892e-16
    
     > rt.b$root ## 3.954775e-16
     [1] 2.220446e-16
     > rt.b$root / eps.c ## 1.781072
     [1] 1
     > ##==> for a < 1.781*eps, bnd > 0 ==> we have log(p) < bnd for all p
     > ## otherwise, we should effectively 'test'
     > curve(.qgammaApprBnd, 1e-16, 1e-10, log="x", col=2)
     > showProc.time()
     Time elapsed: 0 0 0
     >
     >
     > ## source ("/u/maechler/R/MM/NUMERICS/dpq-functions/beta-gamma-etc/qgamma-fn.R")
     > ## ##--> qchisqAppr.R() -- which has 'kind = ' argument!
     > ## ##--> qgamma.R()
     >
     > p.qchi.appr <-
     + function(x, qm= { m <- cbind(qchisq(x, df, log=TRUE),
     + sapply(knds, function(kind)
     + qchisqAppr.R(x,df,log=TRUE,kind=kind)))
     + colnames(m) <- c("True", "default", knds[-1])
     + m },
     + df,
     + knds = list(NULL,"chi.small", "WH", "p1WH", "df.small"),
     + call = match.call(), main = deparse(call), log = "y", ...)
     + {
     + ## Purpose:
     + ## ----------------------------------------------------------------------
     + ## Arguments:
     + ## ----------------------------------------------------------------------
     + ## Author: Martin Maechler, Date: 25 Mar 2004, 22:08
     +
     + col <- c(2,1,3:6)
     + lty <- c(1,3,1,1,1,1)
     + lwd <- c(2,2,1,1,1,1)
     + matplot(x, qm, col=col, lty=lty, lwd=lwd, log = log,
     + type = 'l', main = main, ...)
     + y0 <- c( .02, .98) %*% par('usr')[3:4]
     + if(par("ylog")) y0 <- 10^y0
     + legend(min(x), y0, colnames(qm), col=col, lty=lty, lwd=lwd)
     + invisible(list(x=x, qm=qm, call=match.call()))
     + }
     >
     >
     > pD.chi.appr <- function(pqr, err.kind=c("relative", "absolute"),
     + type = "l", log = "y",
     + lwds = c(2, rep(1, k-1)),
     + cols = seq(along=lwds),
     + ltys = rep(1,k),
     + ...)
     + {
     + ## Purpose: Plot Difference from "True" qchisq()
     + ## ----------------------------------------------------------------------
     + ## Arguments: pqr: a list as resulting from p.chi.appr()
     + ## ----------------------------------------------------------------------
     + ## Author: Martin Maechler, Date: 31 Mar 2004, 09:38
     + err.kind <- match.arg(err.kind)
     + if(!is.list(pqr) || !is.numeric(k <- ncol(pqr$qm)-1) || k <= 0)
     + stop("invalid first argument 'pqr'")
     + with(pqr, {
     + err <- abs(if(err.kind == "relative")
     + (1- qm[,-1] / qm[,1]) else (qm[,-1] - qm[,1]))
     + matplot(x, err, ylab = "",
     + main = paste(err.kind,"Error from", deparse(call)),
     + type= type, log= log, lty=ltys, col=cols, lwd=lwds, ...)
     + legend(par("xaxp")[1], par("yaxp")[2], colnames(qm)[-1],
     + lty=ltys, col=cols, lwd=lwds)
     + })
     + }
     >
     > ## if(FALSE) # you can manually set
     > ## do.pqchi <- TRUE # before source()ing this file
     > ## if(!exists("do.pqchi") || !is.logical(do.pqchi))
     > ## do.pqchi <- interactive()
     >
     > ## if(do.pqchi) { #------------- FIXME look at speed up or cache indeed !?? <<<<<
     >
     > ## pqFile <- "/u/maechler/R/MM/NUMERICS/dpq-functions/pqchi.rda"
     > ## ## ls -l ................... 1325446 Nov 2 2009 pqchi.rda
     > ## if(file.exists(pqFile)) {
     > ## attach(pqFile) ## it loads more than we create here __FIXME__
     > ## print(ls(2, all.names=TRUE))
     > ## ## [1] "pq.1" "pq.25" "pq.25." "pq.31" "pq.33" "pq.33." "pq.33.2"
     > ## ## [8] "pq.33.3" "pq.33.4" "pq.5" "pq.5." "pq1" "pq1." "pq1.95"
     > ## ## [15] "pq1.95." "pq1.95.2" "pq10" "pq10." "pq10.2" "pq100" "pq2"
     > ## ## [22] "pq2." "pq2.05" "pq2.05." "pq2.05.2" "pq2.5" "pq2.5." "pq2.5.2"
     > ## ## [29] "pq200" "pq2L" "pq4" "pq4." "pq4.2" "pqFile"
     > ## }
     >
     > showProc.time()
     Time elapsed: 0.01 0 0.02
     > x <- seq(-300, -.01, length=501)
     >
     > all.equal(qchisqAppr.R(x, 200, log=TRUE),
     + qchisqAppr (x, 200, log=TRUE),tol=0)
     [1] "Mean relative difference: 1.968318e-16"
     > ## 4.48 e-16 / TRUE (Opteron)
     >
     > all.equal(qchisqAppr.R(x, 2, log=TRUE),
     + qchisqAppr (x, 2, log=TRUE),tol=0)
     [1] "Mean relative difference: 1.535551e-16"
     > ## 3.90 e-16 / TRUE (Opteron)
     >
     > all.equal(qchisqAppr.R(x, 0.1, log=TRUE),
     + qchisqAppr (x, 0.1, log=TRUE),tol=0)
     [1] TRUE
     > ## 7.15 e-15 / 1.179e-8 !!!!! (Opteron)
     >
     > pq200 <- p.qchi.appr(x = seq(-300, -.01, length=501), df = 200)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     939 y values <= 0 omitted from logarithmic plot
     > pq100 <- p.qchi.appr(x = seq(-160, -.01, length=501), df = 100)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     941 y values <= 0 omitted from logarithmic plot
     > ## after (slow) computing, quickly repeat:
     > with(pq200, p.qchi.appr(x=x, qm=qm, call=call))
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     939 y values <= 0 omitted from logarithmic plot
     > with(pq100, p.qchi.appr(x=x, qm=qm, call=call))
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     941 y values <= 0 omitted from logarithmic plot
     >
     > ## this "hangs forever" -- before I introduced 'maxit' (for 'nu.small'):
     > pq10 <- p.qchi.appr(x = seq(-12, -.005, length=501), df = 10)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     891 y values <= 0 omitted from logarithmic plot
     > ## want to see the jump:
     > pq10. <- p.qchi.appr(x = seq(-10, -4, length=501), df = 10)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     1002 y values <= 0 omitted from logarithmic plot
     > pq10.2<- p.qchi.appr(x = seq(-8.5,-7.5, length=501), df = 10)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     1002 y values <= 0 omitted from logarithmic plot
     > with(pq10.2, p.qchi.appr(x=x, qm=qm, call=call))
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     1002 y values <= 0 omitted from logarithmic plot
     >
     >
     > pq2.5 <- p.qchi.appr(x = seq(-3.4, -.01, length=501), df = 2.5)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     244 y values <= 0 omitted from logarithmic plot
     > pq2.5. <- p.qchi.appr(x = seq(-2.1, -1.8, length=901), df = 2.5)#the jump
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     901 y values <= 0 omitted from logarithmic plot
     > ## what about p1WH (which is fantastic for df=2)?
     > pq2.5.2<- p.qchi.appr(x = seq(-0.5, -1e-3, length=901), df = 2.5)
     > with(pq2.5, p.qchi.appr(x=x, qm=qm, call=call))
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     244 y values <= 0 omitted from logarithmic plot
     > with(pq2.5.2, p.qchi.appr(x=x, qm=qm, call=call))
     > pD.chi.appr(pq2.5.2)# nothing special
     >
     > pq2.05 <- p.qchi.appr(x = seq(-3.4, -.01, length=501), df = 2.05)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     81 y values <= 0 omitted from logarithmic plot
     > pq2.05. <- p.qchi.appr(x = seq(-2.5, -1.5, length=901), df = 2.05)#the jump
     > ## ^^ the jump from chi.small to WH is much too late here
     > ## what about p1WH (which is fantastic for df=2)?
     > pq2.05.2<- p.qchi.appr(x = seq(-0.4, -1e-5, length=901), df = 2.05)
     > pD.chi.appr(pq2.05.2) # p1WH is starting to become better
     > # and the jump (WH -> p1WH) is too late
     >
     > with(pq2.05, p.qchi.appr(x=x, qm=qm, call=call))
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     81 y values <= 0 omitted from logarithmic plot
     > with(pq2.05.2, p.qchi.appr(x=x, qm=qm, call=call))
     >
     >
     > ## Here, all are 'ok' (but "nu.small"):
     > pq2L <- p.qchi.appr(seq(-300, -.01, length=201), df = 2)
     There were 50 or more warnings (use warnings() to see the first 50)
     >
     > pq2 <- p.qchi.appr(x = seq(-5, -.001, length=501), df = 2)
     > pq2. <- p.qchi.appr(x = seq(-2.5, -1, length=901), df = 2)
     > with(pq2., p.qchi.appr(x=x, qm=qm, call=call))
     >
     > pq4 <- p.qchi.appr(x = seq(-8, -0.01, length = 501), df = 4)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     845 y values <= 0 omitted from logarithmic plot
     > summary(warnings())
     1 identical warnings:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     845 y values <= 0 omitted from logarithmic plot
     > pq4.2 <- p.qchi.appr(x = seq(-0.1, -1e-04, length = 901), df = 4)
     >
     > pq1.95 <- p.qchi.appr(x = seq(-3., -.01, length=501), df = 1.95)
     > pq1.95. <- p.qchi.appr(x = seq(-2.2, -1.5, length=901), df = 1.95)#the jump -1.57
     > ## ^^ the jump from chi.small to WH is *much* too late here
     > ## what about p1WH (which is fantastic for df=2)?
     > pq1.95.2<- p.qchi.appr(x = seq(-0.4, -1e-7, length=901), df = 1.95)
     > pD.chi.appr(pq1.95.2) # p1WH is starting to become better
     > # and the jump (WH -> p1WH) is too late
     >
     > with(pq1.95, p.qchi.appr(x=x, qm=qm, call=call))
     > with(pq1.95.2, p.qchi.appr(x=x, qm=qm, call=call))
     >
     >
     > pq1 <- p.qchi.appr(x = seq(-4, -.001, length=501), df = 1)
     There were 50 or more warnings (use warnings() to see the first 50)
     > pq1. <- p.qchi.appr(x = seq(-1.8, -.5, length=901), df = 1)
     > with(pq1., p.qchi.appr(x=x, qm=qm, call=call))
     >
     > pq.5 <- p.qchi.appr(x = seq(-1.5, -.001, length=501), df = 0.5)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     243 y values <= 0 omitted from logarithmic plot
     > pq.5. <- p.qchi.appr(x = seq(-0.8, -0.2, length=901), df = 0.5, ylim=c(.04,.6))
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     40 y values <= 0 omitted from logarithmic plot
     >
     > pq.33 <- p.qchi.appr(x= seq(-0.9, -.001,length=501), df= 0.33)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     246 y values <= 0 omitted from logarithmic plot
     > pq.33. <- p.qchi.appr(x= seq(-0.4, -0.02,length=901), df= 0.33)
     > pq.33.2<- p.qchi.appr(x= seq(-0.4, -0.2, length=901), df= 0.33, ylim=c(.15,.60))
     > with(pq.33.2, p.qchi.appr(x=x, qm=qm, call=call,ylim=c(.15,.45)))
     >
     > pq.33.3<- p.qchi.appr(x= seq(-0.4, -0.005, length=901), df= 0.33, ylim=c(.15, 4.00))
     > with(pq.33.3, p.qchi.appr(x=x, qm=qm, call=call))#,ylim=c(.15, 8)))
     >
     > pq.33.4<- p.qchi.appr(x= seq(-0.003, -1e-6, length=901), df= 0.33,ylim=c(5,25))
     > with(pq.33.4, p.qchi.appr(x=x, qm=qm, call=call,ylim=c(5,25)))
     >
     > ## nu <= 0.32 is the "magic" border of Best & Roberts
     >
     > pq.31 <- p.qchi.appr(x = seq(-0.45,-.010, length=501), df = 0.31)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     27 y values <= 0 omitted from logarithmic plot
     > with(pq.31, p.qchi.appr(x=x, qm=qm, call=call))
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     27 y values <= 0 omitted from logarithmic plot
     >
     > pq.25 <- p.qchi.appr(x = seq(-0.3, -0.02, length=901), df = 0.25)
     >
     > pq.1 <- p.qchi.appr(x = seq(-0.16, -0.01, length=901), df = 0.1)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     226 y values <= 0 omitted from logarithmic plot
     > with(pq.1, p.qchi.appr(x=x, qm=qm, call=call))
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log = log) :
     226 y values <= 0 omitted from logarithmic plot
     > showProc.time()
     Time elapsed: 9.36 0.08 9.56
     >
     > ## if(!file.exists(pqFile)) # don't overwrite for now (as it contains pq2L ,
     > ## save(list=ls(pat="^pq"), file = pqFile)
     >
     > ##}## end if(do.pqchi){ only if interactive } ======================================
     >
     > pD.chi.appr(pq2L, "abs")
     Warning messages:
     1: In xy.coords(x, y, xlabel, ylabel, log = log) :
     363 y values <= 0 omitted from logarithmic plot
     2: In xy.coords(x, y, xlabel, ylabel, log) :
     180 y values <= 0 omitted from logarithmic plot
     > pD.chi.appr(pq2L, "rel")
     Warning messages:
     1: In xy.coords(x, y, xlabel, ylabel, log = log) :
     363 y values <= 0 omitted from logarithmic plot
     2: In xy.coords(x, y, xlabel, ylabel, log) :
     180 y values <= 0 omitted from logarithmic plot
     > ## --> want only much smaller x-range:
     > pD.chi.appr(pq2,"abs")#--> fantastic p1WH
     Warning messages:
     1: In xy.coords(x, y, xlabel, ylabel, log = log) :
     340 y values <= 0 omitted from logarithmic plot
     2: In xy.coords(x, y, xlabel, ylabel, log) :
     1 y value <= 0 omitted from logarithmic plot
     > pD.chi.appr(pq2) # (ditto)
     Warning messages:
     1: In xy.coords(x, y, xlabel, ylabel, log = log) :
     340 y values <= 0 omitted from logarithmic plot
     2: In xy.coords(x, y, xlabel, ylabel, log) :
     1 y value <= 0 omitted from logarithmic plot
     >
     > pD.chi.appr(pq4.2)# p1WH: only at very right
     > pD.chi.appr(pq4.2, xlim=c(-.016,0))# p1WH: only at very right
     >
     >
     > ## no Newton step here, eg:
     > (qgA100 <- qgammaAppr(1e-100, 100))
     [1] 3.799269
     > (qg.100 <- qgamma (1e-100, 100))
     [1] 3.950799
     > all.equal(qgA100, qg.100)
     [1] "Mean relative difference: 0.03988396"
     > ## too much different
     > dgamma(1e-100, 100, log = TRUE)# -23154.7 i.e., "non-log" is 0
     [1] -23154.73
     >
     > qgamma.R(1e-100, 100, verbose = TRUE)#-> final Newton fails!
     qgamma(p= 1e-100, alpha= 100, scale= 1, l.t.= 1, log.p= 0): (p1,df)=(-230.259,200) ==> small chi-sq., ch0 = 7.59854
    
     ==> ch = 7.59854: Ph.II iter; ch=7.59854, p2=9.76738e-101
     it=2, ch=2.45916e+09, p2=-1
     --> Del became NA
    
     it=1: p=-230.259, x = 3.79927, p.=-234.019; p1:=D{p}=-3.76094
     new t= 3.94774005, p.= -230.332933 it=2, d{p}=-0.074424
     new t= 3.95079758, p.= -230.258539 it=3, d{p}=-2.99766e-05
     new t= 3.95079881, p.= -230.258509 it=4, d{p}=-4.88853e-12
     new t= 3.95079881, p.= -230.258509 it=5, d{p}=0
     [1] 3.950799
     >
     > ## but here, the final Newton iterations do work :
     > x <- qgamma.R(log(1e-100), 100, log = TRUE, verbose = TRUE)
     qgamma(p=-230.259, alpha= 100, scale= 1, l.t.= 1, log.p= 1): (p1,df)=(-230.259,200) ==> small chi-sq., ch0 = 7.59854
    
     it=1: p=-230.259, x = 3.79927, p.=-234.019; p1:=D{p}=-3.76094
     new t= 3.94774005, p.= -230.332933 it=2, d{p}=-0.074424
     new t= 3.95079758, p.= -230.258539 it=3, d{p}=-2.99766e-05
     new t= 3.95079881, p.= -230.258509 it=4, d{p}=-4.88853e-12
     new t= 3.95079881, p.= -230.258509 it=5, d{p}=0
     > pgamma(x, 100) # = 1e-100 ! perfect
     [1] 1e-100
     > showProc.time()
     Time elapsed: 0.11 0 0.11
     >
     > ###--> Use this to devise an improved final Newton algorithm !!!
     >
     > <0c>
     > ## From: Prof Brian Ripley <ripley@stats.ox.ac.uk>
     > ## To: skylab.gupta@gmail.com
     > ## Cc: R-bugs@biostat.ku.dk, r-devel@stat.math.ethz.ch
     > ## Subject: Re: [Rd] qgamma inaccuracy (PR#12324)
     > ## Date: Tue, 12 Aug 2008 20:50:50 +0100 (BST)
     >
     > ## This is a really extreme usage. AFAICS the code works well enough down to
     > ## shape=1e-10 or so, e.g.
     >
     > qgamma(1e-10, 5e-11, lower.tail=FALSE)
     [1] 0.08237203
     > ## [1] 0.08237203
     > ## in R 2.9.. .. 2.10.0 -- with an accuracy warning {which is *wrong*!}
     >
     > ## I would be interested to know what substantive problem you were trying to
     > ## solve here that required such values.
     >
     > ## I am pretty sure that a completely different algorithm will be required.
     >
     > ## MM: It looks like this is basis for a new algo:
     > a <- 1e-14
     > gammaE <- 0.57721566490153286060651209008240243079 # Euler's gamma
     > curve(pgamma(x, a, lower.tail=FALSE)/a + log(x) + gammaE, 1e-300, 1e-1, log="",n=1000)
     > curve(pgamma(x, a, lower.tail=FALSE)/a + log(x) + gammaE - x, 1e-300, 1e-1, log="",n=1000)
     > ## ==> Q = 1 - P = pgamma(x,a, lower=FALSE) ~= a*(-log(x) - gammaE + x - x^2/4)
     > ## i.e., Q ~= -a(log(x) + gammaE { -x + x^2/4 }
     > ## -Q/a - gammaE ~= log(x) { -x + x^2/4 }
     > ## ==> x ~= exp(- (Q/a + gammaE))
     > ## e.g., example below:
     > Q <- 1e-100; a <- 5e-101
     > ## MM: Find inverse :
     > str(r.a <- uniroot(function(x) pgamma(x,a,lower.tail=FALSE) - Q,
     + int = c(0.01, 0.1), tol=1e-20))
     List of 5
     $ root : num 0.0824
     $ f.root : num 0
     $ iter : int 8
     $ init.it : int NA
     $ estim.prec: num 5.95e-11
     > dput(x0 <- r.a$root) ## 0.0823720296207203
     0.0823720296207203
     > (x1 <- exp(- (Q/a + gammaE)))## 0.07598528 .. not so good
     [1] 0.07598528
     > qgammaApprSmallP(Q, a, lower.tail=FALSE)## ~= 0.07598528 -- the same!!
     [1] 0.07598528
     >
     > pgamma(x0, a, lower.tail=FALSE) ## 1.00000e-100.
     [1] 1e-100
     > pgamma(x1, a, lower.tail=FALSE) ## 1.03728e-100 ... hmm "close"
     [1] 1.037283e-100
     > ##
     >
     > ## MM: -- now look at the bigger picture
     > p.qg.2a <- function(l2x.min= -15, l2x.max = -100, n = 501,
     + do.offset = FALSE,
     + type = "o", log = "x", cex = 0.6, ...) {
     + x.log <- any("x" == strsplit(log,"")[[1]])
     + x <- if(x.log) 2^seq(l2x.min, l2x.max, length=n)
     + else seq(2^l2x.min,2^l2x.max, length=n)
     + if(do.offset)
     + plot(x, qgamma(2*x, x, lower.tail=FALSE) - 0.0823720296206873,
     + type=type, cex=cex, log=log, ...)
     + else plot(x, qgamma(2*x, x, lower.tail=FALSE),
     + type=type, cex=cex, log=log, ...)
     + }
     > p.qg.2a() # was "very bad" in R <= 2.10.0, now --> 0.082372... "perfect" smooth
     > ## still a little ---but acceptable--- remaining inaccuracy ...zooming in:
     > p.qg.2a(-43,-55, do.offset=TRUE)
     > p.qg.2a(,-1024)
     > p.qg.2a(,-1024, log="", pch=".")## linear in x !!
     > ## zoom in at the limit
     > p.qg.2a(-30,-1024, do.offset=TRUE, ylim = 1e-11*c(-1,1))
     > p.qg.2a(-33,-1024, do.offset=TRUE, ylim = 1e-12*c(-1,1))
     > p.qg.2a(-33,-1024, do.offset=TRUE, ylim = 1e-13*c(-1,1))
     >
     > a <- 2^-(7:900)
     > qg <- qgamma(2*a, a, lower.tail=FALSE)
     > re <- 1-pgamma(qg, a, lower.tail=FALSE)/(2*a)
     > plot(a, re, log="x", type="b", col=2)
     > stopifnot(abs(re) < 2e-12) # but really, *should be a bit better
     >
     > showProc.time()
     Time elapsed: 0.14 0 0.15
     > ## For completeness we may write that in due course, but for now (R 2.7.2) I
     > ## suggest just issuing a warning for miniscule 'shape'.
     >
     > ## On Thu, 7 Aug 2008, skylab.gupta@gmail.com wrote:
     >
     > ## > Full_Name:
     > ## > Version: 2.7.1 (2008-06-23)
     > ## > OS: windows vista
     > ## > Submission from: (NULL) (216.82.144.137)
     > ## >
     > ## >
     > ## > Hello,
     > ## >
     > ## > I have been working with various probability distributions in R, and it seems
     > ## > the gamma distribution is inaccurate for some inputs.
     >
     > ## > For example, qgamma(1e-100, 5e-101, lower.tail=FALSE) gives: 1.0. However, it
     > ## > seems this is incorrect; I think the correct answer should be
     > ## > 0.082372029620717283. When I check these numbers using pgamma, I get:
     >
     > (qg <- qgamma(1e-100, 5e-101, lower.tail=FALSE))# 1 (wrong, originally)
     [1] 0.08237203
     > ## 0.08237203 now (2009-11-04), i.e. ok
     > pgamma(qg, 5e-101, lower.tail=FALSE)# now -> 1e-100 : ok
     [1] 1e-100
     >
     > pgamma(0.082372029620717283,5e-101, lower.tail=FALSE)
     [1] 1e-100
     > ## 1.0000000000000166e-100.
     >
     > RE.pqgamma <- function(p, shape, lower.tail = TRUE, log.p = FALSE) {
     + ## Relative Error of pgamma(qgamma(*), ..):
     + 1 - pgamma(qgamma(p, shape, lower.tail=lower.tail, log.p=log.p),
     + shape=shape, lower.tail=lower.tail, log.p=log.p) / p
     + }
     > RE.qpgamma <- function(q, shape, lower.tail = TRUE, log.p = FALSE) {
     + ## Relative Error of qgamma(pgamma(*), ..):
     + 1 - qgamma(pgamma(q, shape, lower.tail=lower.tail, log.p=log.p),
     + shape=shape, lower.tail=lower.tail, log.p=log.p) / q
     + }
     >
     > ## Ok, how extreme can we get -- let a := alpha := shape --> 0 :
     > x <- 1e-100
     > a <- 2^-(7:300)# is still "ok":
     > plot(a, (re <- RE.pqgamma(x, a, lower.tail=FALSE)), type="b", col=2, log="x")# oops!
     > a <- 2^-(7:400)
     > plot(a, (re <- RE.pqgamma(x, a, lower.tail=FALSE)), type="b", col=2, log="x")# oops!
     > ## Oops!
     > ## but, it is clear
     > qgamma(x, 2^-400, lower.tail=FALSE)## is exactly 0
     [1] 0
     >
     > ## -> it goes to 0 quickly . .. zooming in:
     > curve(qgamma(1e-100, x, lower.tail=FALSE), 1e-110, 1e-70, log="xy", col=2, axes=FALSE)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     18 y values <= 0 omitted from logarithmic plot
     > eaxis(1);eaxis(2)
     >
     > ## from when on is it exactly 0:
     > uniroot(function(u) qgamma(1e-100, 2^u,lower.tail=FALSE)-1e-315, c(-400, -300))$root
     [1] -341.6941
     > ## -341.6941
     > ## => use
     > a <- 2^-(7:341)
     > plot(a, (re <- RE.pqgamma(x, a, lower.tail=FALSE)),
     + type="l", col=2, log="x")# small glips
     > ## zoom in:
     > curve(abs(RE.pqgamma(1e-100, x, lower.tail=FALSE)), 2^-341, 1e-90, log="xy", n=2000)
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     6 y values <= 0 omitted from logarithmic plot
     > curve(RE.pqgamma(1e-100, x, lower.tail=FALSE), 1e-100, 10e-100, n=2000)
     > curve(RE.pqgamma(1e-100, x, lower.tail=FALSE), 4e-100, 6e-100, n=2000)
     > ## Ok: at least here is a problem
     > RE.pqgamma(1e-100, 5e-100, lower.tail=FALSE)# -0.1538
     [1] 3.874678e-14
     >
     > ## more general
     > curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 1e-100, 1e-20, n=10000, log="x")
     > ## problem *everywhere* , starting quite early: (lesser problem at ~ 1e-16 !)
     > curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 1e-25, 1e-15, n=1000, log="x")
     > curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 1e-21, 10e-21,n=1000)
     > curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 2e-21, 6e-21, n=1000)
     > curve(RE.pqgamma(x, 5*x, lower.tail=FALSE), 4e-21, 4.5e-21, n=1000)
     > ## and indeed, it's qgamma() that jumps here
     > curve(qgamma(x, 5*x, lower.tail=FALSE), 4e-21, 4.5e-21, ylim=c(.97, 1.1))
     >
     > ## well, looking at pgamma(), finally reveals the buglet is there first:
     > ## There's a jump at x = 1 !!
     > curve(pgamma(x, 1e-30, lower=FALSE), .9999, 1.0001)
     > curve(pgamma(x, 1e-20, lower=FALSE), .9999, 1.0001)
     > curve(pgamma(x, 1e-17, lower=FALSE), .9999, 1.0001)
     > curve(pgamma(x, 1e-15, lower=FALSE), .9999, 1.0001)
     > curve(pgamma(x, 1e-13, lower=FALSE), .9999, 1.0001)
     > curve(pgamma(x, 1e-12, lower=FALSE), .9999, 1.0001)# still clearly visible
     > curve(pgamma(x, 1e-11, lower=FALSE), .9999, 1.0001)# barely visible
     > ## for larger alpha == shape, must zoom in more and more:
     > curve(pgamma(x, 1e-11, lower=FALSE), .99999, 1.00001)#
     > curve(pgamma(x, 1e-10, lower=FALSE), .999999, 1.000001)
     > curve(pgamma(x, 1e-9, lower=FALSE), .9999999, 1.0000001)
     > curve(pgamma(x, 1e-8, lower=FALSE), .99999999, 1.00000001)
     > curve(pgamma(x, 1e-7, lower=FALSE), .999999999, 1.000000001)
     > curve(pgamma(x, 1e-6, lower=FALSE), .9999999999, 1.0000000001)
     > curve(pgamma(x, 1e-5, lower=FALSE), .99999999999, 1.00000000001)
     > curve(pgamma(x, 1e-4, lower=FALSE), .999999999999, 1.000000000001)
     > ## now we get close to noise level:
     > curve(pgamma(x, 1e-3, lower=FALSE), .9999999999999, 1.0000000000001)
     > curve(pgamma(x, 1e-2, lower=FALSE), .99999999999999, 1.00000000000001)
     Warning message:
     In plot.window(...) :
     relative range of values = 67 * EPS, is small (axis 1)
     >
     > showProc.time()
     Time elapsed: 0.31 0.07 0.44
     >
     > del.pgamma <- function(a, eps = 1e-13)
     + {
     + ## Purpose: *relative* jump size at x = 1 of pgamma(x, a, lower=FALSE)
     + ## ----------------------------------------------------------------------
     + ## Author: Martin Maechler, Date: 5 Nov 2009, 16:08
     + stopifnot(a > 0, length(a) == 1, eps > 0, length(eps) == 1)
     + pp <- pgamma(1+c(-eps,eps), a, lower.tail = FALSE)
     + (pp[2] - pp[1])*2/(pp[2] + pp[1])
     + }
     >
     > a <- lseq(1e-300, 1e-3, length=400)
     > dpa <- sapply(a, del.pgamma)
     > plot(a, -dpa, log="xy", type="l", col=2, yaxt="n");eaxis(2)
     > ## ok, it remains constant all the way to 1e-300
     > ## --> focus
     > a <- lseq(1e-40, 1e-5, length=400)
     > dpa <- sapply(a, del.pgamma)
     > plot(a, -dpa, log="xy", type="l", col=2, axes=FALSE)
     > eaxis(1, at = 10^-seq(5,40, by=5))
     > eaxis(2)
     >
     >
     > xm <- .Machine$double.xmin
     > pgamma(xm, shape= 1e-20)# is "practically 1" --> *most* qgamma() will be exactly 0
     [1] 1
     > ## how close to 1 ? ---> use upper_tail [possibly on log scale:]
     > pgamma(xm, shape= 1e-20, lower.tail=FALSE) # 7.078192e-18
     [1] 7.078192e-18
     > pgamma(xm, shape= 1e-20, lower.tail=FALSE, log=TRUE) # -39.48951
     [1] -39.48951
     >
     > ## Where is the 'boundary' (from which on qgamma() must return 0, since it can't give
     > ## xm = 2.2e-308 ) ?
     > a <- 2^-(7:1030)
     > plot(a, (p <- pgamma(xm, a, lower.tail=FALSE, log=TRUE)),
     + cex=.5,type ="l", col=2, log="x")
     > summary(lm. <- lm(p ~ log(a) + a + I(a^2))) ## coeff. of log(a) *is* 1
    
     Call:
     lm(formula = p ~ log(a) + a + I(a^2))
    
     Residuals:
     Min 1Q Median 3Q Max
     -0.0081566 -0.0000098 0.0000171 0.0000440 0.0107065
    
     Coefficients:
     Estimate Std. Error t value Pr(>|t|)
     (Intercept) 6.562e+00 3.327e-05 197240 <2e-16 ***
     log(a) 1.000e+00 8.024e-08 12463221 <2e-16 ***
     a -3.403e+02 2.039e-01 -1669 <2e-16 ***
     I(a^2) 1.551e+04 2.911e+01 533 <2e-16 ***
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     Residual standard error: 0.0005225 on 1020 degrees of freedom
     Multiple R-squared: 1, Adjusted R-squared: 1
     F-statistic: 5.249e+13 on 3 and 1020 DF, p-value: < 2.2e-16
    
     > summary(lm. <- lm(p ~ offset(log(a)) + a + I(a^2)))
    
     Call:
     lm(formula = p ~ offset(log(a)) + a + I(a^2))
    
     Residuals:
     Min 1Q Median 3Q Max
     -0.0081788 0.0000179 0.0000179 0.0000179 0.0107351
    
     Coefficients:
     Estimate Std. Error t value Pr(>|t|)
     (Intercept) 6.562e+00 1.639e-05 400476.3 <2e-16 ***
     a -3.404e+02 2.033e-01 -1674.4 <2e-16 ***
     I(a^2) 1.552e+04 2.907e+01 533.7 <2e-16 ***
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     Residual standard error: 0.0005232 on 1021 degrees of freedom
     Multiple R-squared: 1, Adjusted R-squared: 1
     F-statistic: 7.853e+13 on 2 and 1021 DF, p-value: < 2.2e-16
    
     >
     > p.l <- pgamma(xm, a, lower.tail=FALSE, log=TRUE) - log(a)
     > dput(mean(tail(p.l))) ## 6.56218869790132
     6.56218869790132
     > ##=> pgamma(xm, a) ~= log(a) + 6.5621886979 +
     > ## ok, to get this better, now need different a:
     > al <- seq(1e-300, 1e-3, length=200)
     > plot(al, (pl <- pgamma(xm, al, lower.tail=FALSE, log=TRUE)) - (log(al)+6.56218869790132),
     + cex=.5,type ="l", col=2)
     > summary(lm.. <- lm(pl ~ offset(log(al) + 6.56218869790132) + 0 + al + I(al^2)))
    
     Call:
     lm(formula = pl ~ offset(log(al) + 6.56218869790132) + 0 + al +
     I(al^2))
    
     Residuals:
     Min 1Q Median 3Q Max
     -1.201e-05 -4.268e-06 -1.336e-06 3.075e-06 5.247e-06
    
     Coefficients:
     Estimate Std. Error t value Pr(>|t|)
     al -3.539e+02 2.032e-03 -174123 <2e-16 ***
     I(al^2) 2.075e+04 2.617e+00 7929 <2e-16 ***
     ---
     Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
     Residual standard error: 4.153e-06 on 198 degrees of freedom
     Multiple R-squared: 1, Adjusted R-squared: 1
     F-statistic: 1.359e+16 on 2 and 198 DF, p-value: < 2.2e-16
    
     > confint(lm..)
     2.5 % 97.5 %
     al -353.8628 -353.8547
     I(al^2) 20745.6597 20755.9814
     > coef(lm..)
     al I(al^2)
     -353.8587 20750.8205
     > ## al I(al^2)
     > ## -353.8587 20750.8205
     > ##=> pgamma(xm, a) ~= log(a) + 6.5621886979 - 353.858745 * a
     >
     > ## > Similarly, for example: -- now (2009-11-04) ok
     > qgamma(1e-100,0.005,lower.tail=FALSE) # = 109.36757177 instead of 219.5937..
     [1] 219.5937
     > pgamma(109.36757177007101, 0.005, lower.tail=FALSE)# = 1.4787306506694758e-52.
     [1] 1.478731e-52
     >
     > ## > This looks completely wrong. The correct value, I think, should be
     > ## > 219.59373661415756. In fact,
     > pgamma(219.59373661415756, 0.005, lower.tail=FALSE)# = 9.9999999999999558e-101.
     [1] 1e-100
     > ## >
     > ## > In fact, when I do the following in R, the results are completely wrong,
     > ## >
     > a <- 5*10^-(1:40)
     > z1 <- qgamma(1e-100,a,lower.tail=FALSE)
     > (y <- pgamma(z1,a,lower.tail=FALSE))
     [1] 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100
     [11] 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100
     [21] 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100
     [31] 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100
     > ## The value of y that I get should be close to 1e-100, but they are not:
     > ## [1] 1.000000e-100 1.871683e-51 1.478731e-52 1.444034e-53 1.440606e-54
     > ## [6] 1.440264e-55 1.440230e-56 1.440226e-57 1.440226e-58 1.440226e-59
     > summary(abs(1 - y/1e-100))
     Min. 1st Qu. Median Mean 3rd Qu. Max.
     1.332e-15 5.579e-15 9.992e-15 1.623e-14 2.104e-14 6.717e-14
     > plot(a, abs(1 - y/1e-100), log="xy", type="b")
     > stopifnot(abs(1 - y/1e-100) < 2e-13)# max (32b Linux, P.M) = 4.186e-14
     >
     > ## > The correct values of z1 should be:
     > z1true <- c(226.97154111939946, 222.15218724493326, 219.59373661415756,
     + 217.27485383840451, 214.98015408183574, 212.68797118872064,
     + 210.39614286838227, 208.10445550564617, 205.81289009100664,
     + 203.52144711679352)
     > all.equal(z1[1:10], z1true, tol=0)# 1.307e-16 on 32-bit (Pentium M); 1.686e-16 (64b Lnx, 2019)
     [1] "Mean relative difference: 1.68599e-16"
     > stopifnot(all.equal(z1[1:10], z1true, tol = 1e-15))
     > showProc.time()
     Time elapsed: 0.24 0.02 0.27
     > ##>
     > ##> With these values of z1true, we get the expected values:
     >
     > (y <- pgamma(z1true, a, lower.tail=FALSE))
     [1] 1.000000e-100 1.000000e-100 1.000000e-100 1.000000e-100 1.000000e-100
     [6] 1.000000e-100 1.000000e-100 1.000000e-100 1.000000e-100 1.000000e-100
     [11] 5.869617e-112 7.428625e-111 9.705940e-111 9.970232e-111 9.997024e-111
     [16] 9.999703e-111 9.999970e-111 9.999997e-111 1.000000e-110 1.000000e-110
     [21] 5.869617e-122 7.428625e-121 9.705940e-121 9.970232e-121 9.997024e-121
     [26] 9.999703e-121 9.999970e-121 9.999997e-121 1.000000e-120 1.000000e-120
     [31] 5.869617e-132 7.428625e-131 9.705940e-131 9.970232e-131 9.997024e-131
     [36] 9.999703e-131 9.999970e-131 9.999997e-131 1.000000e-130 1.000000e-130
     > ## [1] 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100 1e-100
     >
     > ## > I am using the precompiled binary version of R, under Windows Vista.
     > ## > -----------
     > ## >> version
     > ## > _
     > ## > platform i386-pc-mingw32
     > ## > arch i386
     > ## > os mingw32
     > ## > system i386, mingw32
     > ## > status
     > ## > major 2
     > ## > minor 7.1
     > ## > year 2008
     > ## > month 06
     > ## > day 23
     > ## > svn rev 45970
     > ## > language R
     > ## > version.string R version 2.7.1 (2008-06-23)
     > ## > ------------
     > ## >
     > ## > So, it seems qgamma is inaccurate for small probability values in the upper
     > ## > tail, when the shape parameter is also small.
     >
     >
     > ###_-- MM: Still wrong:
     >
     > (xm <- 2^-1074.9999) # is less than .Machine $ double.xmin == the really x > 0
     [1] 4.940656e-324
     >
     > pgamma(xm, .00001)# 0.992589
     [1] 0.992589
     > qgamma(.99, .00001)##--> NaN -- should give 0 or "xmin" or so
     [1] 0
     > ## FIXME -- ok, now
     >
     > ## but
     > curve(qgamma(x, .001, lower=FALSE), .4, .8, n=1001, log="y")
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     688 y values <= 0 omitted from logarithmic plot
     > curve(qgamma(x, 1e-5, lower=FALSE), .002, .2, n=1001, log="xy")
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     716 y values <= 0 omitted from logarithmic plot
     > curve(qgamma(x, 1e-7, lower=FALSE), 1e-5, .04, n=1001, log="xy")
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     759 y values <= 0 omitted from logarithmic plot
     > curve(qgamma(x, 1e-12, lower=FALSE), 1e-12, 1e-2, n=1001, log="xy")
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     713 y values <= 0 omitted from logarithmic plot
     >
     > ## or
     > curve(qgamma(x, 1e-121, lower=FALSE), 7e-119, 8e-119,
     + n=2001, log="y", yaxt="n")
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     1117 y values <= 0 omitted from logarithmic plot
     > try( # reveals eaxis() bug ? -- for the *subnormal* numbers
     + eaxis(2, at = 10^-seq(304,324, by=2))
     + )
     Error in eaxis(2, at = 10^-seq(304, 324, by = 2)) :
     invalid 'log=TRUE' for at <= 0: not a true log scale plot?
     >
     > curve(qgamma(x, .001, lower=FALSE), .4, .6, n=1001, log="y")
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     376 y values <= 0 omitted from logarithmic plot
     > curve(qgamma(x, .001, lower=FALSE), .5, .55, n=1001, log="y")
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     503 y values <= 0 omitted from logarithmic plot
     > ## gives a *warning* from axis() because of subnormal y-range {was error, fixed in R-devel (2018-08)}
     > curve(qgamma(x, .001, lower=FALSE), .52, .53, n=1001, log="y")
     Error in axis(side = side, at = at, labels = labels, ...) :
     CreateAtVector [log-axis()]: axp[0] = 0 < 0!
     Calls: curve ... plot.default -> localAxis -> Axis -> Axis.default -> axis
     In addition: Warning messages:
     1: In xy.coords(x, y, xlabel, ylabel, log) :
     514 y values <= 0 omitted from logarithmic plot
     2: In plot.window(...) : Internal(pretty()): very small range.. corrected
     3: In axis(side = side, at = at, labels = labels, ...) :
     CreateAtVector "log"(from axis()): axp[0] = 0 !
     Execution halted
Flavor: r-oldrel-windows-ix86+x86_64

Version: 0.3-5
Check: tests
Result: ERROR
     Running ‘chisq-nonc-ex.R’ [22s/22s]
     Running ‘dnchisq-tst.R’ [0s/0s]
     Running ‘pnbeta-tst.R’ [0s/0s]
     Running ‘pnt-precision-problem.R’ [11s/11s]
     Running ‘ppois-ex.R’ [1s/1s]
     Running ‘qPoisBinom-ex.R’ [0s/0s]
     Running ‘qbeta-dist.R’ [11s/11s]
     Running ‘qbeta-tst.R’ [0s/0s]
     Running ‘qgamma-ex.R’ [14s/15s]
    Running the tests in ‘tests/qgamma-ex.R’ failed.
    Last 13 lines of output:
     Warning message:
     In xy.coords(x, y, xlabel, ylabel, log) :
     503 y values <= 0 omitted from logarithmic plot
     > ## gives a *warning* from axis() because of subnormal y-range {was error, fixed in R-devel (2018-08)}
     > curve(qgamma(x, .001, lower=FALSE), .52, .53, n=1001, log="y")
     Error in axis(side = side, at = at, labels = labels, ...) :
     CreateAtVector [log-axis()]: axp[0] = 0 < 0!
     Calls: curve ... plot.default -> localAxis -> Axis -> Axis.default -> axis
     In addition: Warning messages:
     1: In xy.coords(x, y, xlabel, ylabel, log) :
     514 y values <= 0 omitted from logarithmic plot
     2: In plot.window(...) : Internal(pretty()): very small range.. corrected
     3: In axis(side = side, at = at, labels = labels, ...) :
     CreateAtVector "log"(from axis()): axp[0] = 0 !
     Execution halted
Flavor: r-oldrel-osx-x86_64