R Under development (unstable) (2024-11-15 r87338 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ### dgamma(): ----------------------- was part of  ./special-fun-ex.R -------------------
> stopifnot(require("Rmpfr"))
Loading required package: Rmpfr
Loading required package: gmp

Attaching package: 'gmp'

The following objects are masked from 'package:base':

    %*%, apply, crossprod, matrix, tcrossprod

C code of R package 'Rmpfr': GMP using 64 bits per limb


Attaching package: 'Rmpfr'

The following object is masked from 'package:gmp':

    outer

The following objects are masked from 'package:stats':

    dbinom, dgamma, dnbinom, dnorm, dpois, dt, pnorm

The following objects are masked from 'package:base':

    cbind, pmax, pmin, rbind

> require(sfsmisc)# -> eaxis(); relErrV()
Loading required package: sfsmisc

Attaching package: 'sfsmisc'

The following objects are masked from 'package:gmp':

    factorize, is.whole

> 
> (doExtras <- Rmpfr:::doExtras())
[1] FALSE
> options(nwarnings = 50000, width = 99)
> 
> ##                                      vvvvvvvvvvvvvvvv
> ## to enhance  |rel.Err| plots:  from ./special-fun-ex.R, also in ~/R/Pkgs/DPQ/tests/pow-tst.R }
> drawEps.h <- function(p2 = -(53:51), side = 4, lty=3, lwd=2, col=adjustcolor(2, 1/2)) {
+     abline(h = 2^p2, lty=lty, lwd=lwd, col=col)
+     axis(side, las=2, line=-1, at = 2^p2,
+          labels = as.expression(lapply(p2, function(p) substitute(2^E, list(E=p)))),
+          col.axis = col, col=NA, col.ticks=NA)
+ }
> 
> (do.pdf <- !dev.interactive(orNone = TRUE))
[1] TRUE
> 
> if(do.pdf) pdf("special-fun-dgamma.pdf")
> 
> xe <- c(-2e5, -1e5, -2e4, -1e4, -2000, -1000, -500, -200, -100, -50, -20, -10)
> (xe <- c(xe, -8:8, -rev(xe)))
 [1] -2e+05 -1e+05 -2e+04 -1e+04 -2e+03 -1e+03 -5e+02 -2e+02 -1e+02 -5e+01 -2e+01 -1e+01 -8e+00
[14] -7e+00 -6e+00 -5e+00 -4e+00 -3e+00 -2e+00 -1e+00  0e+00  1e+00  2e+00  3e+00  4e+00  5e+00
[27]  6e+00  7e+00  8e+00  1e+01  2e+01  5e+01  1e+02  2e+02  5e+02  1e+03  2e+03  1e+04  2e+04
[40]  1e+05  2e+05
> two <- mpfr(2, 256)
> ## For centering at E[.], will use xP(x, shp) :
> xP <- function(x, d) {
+     ## cannot eliminate them, as for <mpfr> they are all finite ..
+     ## x <- x[is.finite(x)]
+     x - d*(x > d)
+ }
> aEQformat <- function(xy, ...) format(xy, digits = 7, ...)
> allEQ_0 <- function (target, current, ...)
+     all.equal(target, current, tolerance = 0, formatFUN = aEQformat, ...)
> stopIfNot <-
+     if("allow.logical0" %in% names(formals(stopifnot))) { # experimental (MM only)
+         stopifnot
+     } else function(exprs, allow.logical0) stopifnot(exprs=exprs)
> abs19 <- function(r) pmax(abs(r), 1e-19) # cut  |err| to positive {for log-plots}
> 
> for(shp in 2^c(-20, -3, -1:1, 4, 10, 14, 20, 50)) {
+     cat("shape = 2^", log2(shp), ":\n-------------\n")
+     d.dg  <- dgamma(xP(2 ^ xe, shp) -> x, shape=shp)
+     m.dg  <- dgamma(xP(two^xe, shp), shape=shp)
+     m.ldg <- dgamma(xP(two^xe, shp), shape=shp, log=TRUE)
+     relE <- asNumeric(relErrV(m.dg, d.dg))
+     ## Plots: do *not* observe any problems yet
+     plot(x, relE, log="x", type="l",
+          main = paste0("rel.Errors dgamma(., shape = 2^", log2(shp),")"))
+     abline(h=0, col=adjustcolor("gray10", 1/2), lty=3, lwd=2)
+     plot(x, abs19(relE), log="xy", type="l", ylim = pmax(4e-17, range(abs19(relE), finite=TRUE)))
+     abline(h = 2^-(52:50), col=adjustcolor("red4",1/2), lty=3)
+     ##
+     stopIfNot(exprs = {
+         !is.unsorted(xe)
+         is.finite(m.dg)
+         m.dg >= 0
+         shp > 1  || all(diff(m.dg) <= 0)
+         shp > 100|| all((m.dg > 0) >= (d.dg > 0))
+         any(fin.d <- is.finite(d.dg))
+         m.dg[!fin.d] > 1e300
+         { cat("all.EQ(<mpfr>, <doubl>):", allEQ_0  (m.dg[fin.d], d.dg[fin.d]), "\n")
+           shp > 100  ||                   all.equal(m.dg[fin.d], d.dg[fin.d],
+                                                     tol = 1e-13) # 2.063241e-14
+         }
+         ## compare with log scale :
+         if(any(pos.d <- m.dg > 0)) {
+             cat("For non-0 <mpfr>-values;  all.EQ(log(d), d*(log)):",
+                 allEQ_0  (log(m.dg[pos.d]), m.ldg[pos.d]),"\n")
+             ##
+             all.equal(log(m.dg[pos.d]), m.ldg[pos.d], tol = 1e-14)
+         } else TRUE
+     })#, allow.logical0 = TRUE)
+ }
shape = 2^ -20 :
-------------
all.EQ(<mpfr>, <doubl>): Mean relative difference: 5.883928e-16 
For non-0 <mpfr>-values;  all.EQ(log(d), d*(log)): Mean relative difference: 8.635607e-74 
shape = 2^ -3 :
-------------
all.EQ(<mpfr>, <doubl>): Mean relative difference: 2.037429e-16 
For non-0 <mpfr>-values;  all.EQ(log(d), d*(log)): Mean relative difference: 1.058572e-77 
shape = 2^ -1 :
-------------
all.EQ(<mpfr>, <doubl>): Mean relative difference: 1.359070e-17 
For non-0 <mpfr>-values;  all.EQ(log(d), d*(log)): Mean relative difference: 1.149299e-77 
shape = 2^ 0 :
-------------
all.EQ(<mpfr>, <doubl>): Mean relative difference: 1.867682e-17 
For non-0 <mpfr>-values;  all.EQ(log(d), d*(log)): Mean relative difference: 4.108006e-77 
shape = 2^ 1 :
-------------
all.EQ(<mpfr>, <doubl>): Mean relative difference: 1.133307e-16 
For non-0 <mpfr>-values;  all.EQ(log(d), d*(log)): TRUE 
shape = 2^ 4 :
-------------
all.EQ(<mpfr>, <doubl>): Mean relative difference: 1.918428e-16 
For non-0 <mpfr>-values;  all.EQ(log(d), d*(log)): TRUE 
shape = 2^ 10 :
-------------
all.EQ(<mpfr>, <doubl>): Mean relative difference: 1.291409e-16 
For non-0 <mpfr>-values;  all.EQ(log(d), d*(log)): TRUE 
shape = 2^ 14 :
-------------
all.EQ(<mpfr>, <doubl>): Mean relative difference: 1.000000 
For non-0 <mpfr>-values;  all.EQ(log(d), d*(log)): TRUE 
shape = 2^ 20 :
-------------
all.EQ(<mpfr>, <doubl>): Mean relative difference: 7.312341e-17 
For non-0 <mpfr>-values;  all.EQ(log(d), d*(log)): TRUE 
shape = 2^ 50 :
-------------
all.EQ(<mpfr>, <doubl>): Mean relative difference: 2.657575e-18 
For non-0 <mpfr>-values;  all.EQ(log(d), d*(log)): TRUE 
There were 20 warnings (use warnings() to see them)
> 
> ## NB:  dgamma(x, sh)  sh >= 1 calls
> ## --   dpois_raw(sh-1, x)   which then
> ## adds stirlerr(x.) to bd0(x., lambda) ;  where x. <- sh-1; lambda <- x
> ## bd0(x,L) ~= 0 iff x ~= L  <==> sh-1 ~= x  <==>  x+1 ~= sh
> 
> sh2x_gamma <- function(sh, nx, k = 12, f1 = 0.5, f2 = 1.25) {
+     stopifnot(is.numeric(sh), length(sh) == 1, sh >= 0, length(k) == 1, k >= 3,
+               f1 < f2, length(f1) == length(f2), length(f2) == 1)
+     p2 <- 2^-(k:3)
+     1 + sh* unique(sort(c(1-p2, 1, 1+p2, # <- values x very close to sh -- does *not* make any diff (????)
+                           seq(f1, f2, length=nx))))
+ }
> 
> relEgamma <- function(sh, nx = 1001, k = 12, precBits = 256,
+                       x = sh2x_gamma(sh, nx=nx, k=k)) {
+     dg  <- dgamma(x, sh)
+     dgM <- dgamma(mpfr(x,  precBits),
+                   mpfr(sh, precBits))
+     structure(cbind(x, relE = asNumeric(relErrV(dgM, dg))), shape=sh)
+ }
> 
> shs <- 1/32 + seq(6, 16, by = 1/8)
> stopifnot(all(!is.whole(shs*2))) # shs are *not* half-integers
> system.time(
+     LrelE <- lapply(shs, relEgamma)
+ ) # 7.5 sec
   user  system elapsed 
  14.16    0.02   14.18 
> 
> m.relE <- sapply(LrelE, function(m) m[,"relE"])
> qrelE <- t(apply(abs(m.relE), 2, quantile, probs = c(.05, .25, .50, .75, .90, 1)))
> ##               ^^^^^^^^^^^
> 
> ## Heureka! --- this shows quite a difference between R 4.3.3  and R-devel (R 4.4.0) !!
> 
> iS <- sort.list(qrelE[,"50%"], decreasing = TRUE)
> cbind(shs, qrelE)[iS,]
           shs           5%          25%          50%          75%          90%         100%
 [1,]  6.15625 3.406154e-16 5.695360e-16 7.918756e-16 9.980122e-16 1.146722e-15 1.422691e-15
 [2,] 12.40625 3.264512e-17 2.201838e-16 4.556529e-16 6.763920e-16 8.159771e-16 1.237976e-15
 [3,] 13.53125 3.780237e-17 2.189528e-16 4.552480e-16 6.797057e-16 8.226962e-16 1.142693e-15
 [4,] 10.28125 3.861707e-17 2.099343e-16 4.524258e-16 6.717071e-16 8.281026e-16 1.141152e-15
 [5,] 10.03125 4.227963e-17 2.100465e-16 4.516176e-16 6.539794e-16 7.962860e-16 1.259820e-15
 [6,] 12.53125 3.421072e-17 2.088431e-16 4.506959e-16 6.721058e-16 8.413170e-16 1.251471e-15
 [7,] 14.15625 5.953303e-17 2.072735e-16 4.486277e-16 6.723762e-16 8.468715e-16 1.188830e-15
 [8,] 15.90625 3.665735e-17 2.151807e-16 4.471476e-16 6.752052e-16 8.224428e-16 1.135756e-15
 [9,] 11.03125 3.727751e-17 2.150808e-16 4.467589e-16 6.720609e-16 8.059952e-16 1.124140e-15
[10,] 11.78125 3.788590e-17 2.134320e-16 4.452910e-16 6.620367e-16 8.132118e-16 1.138367e-15
[11,] 11.90625 4.240331e-17 2.072419e-16 4.438902e-16 6.553498e-16 8.094856e-16 1.049748e-15
[12,] 14.40625 3.469252e-17 2.256388e-16 4.420190e-16 6.680521e-16 8.113211e-16 1.232621e-15
[13,] 11.28125 4.390362e-17 2.123472e-16 4.404865e-16 6.548729e-16 8.115886e-16 1.131573e-15
[14,] 12.15625 4.285183e-17 2.093994e-16 4.382577e-16 6.633445e-16 8.156693e-16 1.082400e-15
[15,] 10.90625 4.233044e-17 2.121453e-16 4.365036e-16 6.467968e-16 7.956220e-16 1.091077e-15
[16,] 16.03125 3.604725e-17 1.940965e-16 4.363393e-16 6.609130e-16 8.092824e-16 1.182471e-15
[17,] 10.65625 3.410572e-17 1.875870e-16 4.363278e-16 6.857693e-16 8.146012e-16 1.193253e-15
[18,]  9.03125 4.372773e-17 2.086370e-16 4.357133e-16 6.598536e-16 8.076996e-16 1.153380e-15
[19,] 14.90625 3.571820e-17 2.103234e-16 4.356946e-16 6.636308e-16 8.142057e-16 1.139190e-15
[20,] 11.65625 4.147000e-17 2.118409e-16 4.349962e-16 6.574708e-16 8.054285e-16 1.084827e-15
[21,]  9.78125 3.639984e-17 1.945526e-16 4.341260e-16 6.578553e-16 8.118468e-16 1.119916e-15
[22,] 11.15625 3.356890e-17 1.976659e-16 4.327397e-16 6.652691e-16 8.337820e-16 1.202934e-15
[23,] 13.65625 4.800498e-17 2.045067e-16 4.315864e-16 6.424374e-16 8.009232e-16 1.158199e-15
[24,] 15.65625 4.083897e-17 2.081412e-16 4.304377e-16 6.543286e-16 8.186140e-16 1.075301e-15
[25,] 15.53125 3.805408e-17 2.166847e-16 4.292272e-16 6.359984e-16 7.834212e-16 1.071912e-15
[26,] 15.40625 3.878335e-17 1.866844e-16 4.287511e-16 6.699912e-16 7.975596e-16 1.086243e-15
[27,]  9.90625 3.781109e-17 2.040648e-16 4.284148e-16 6.555577e-16 8.131888e-16 1.094031e-15
[28,] 13.15625 3.358864e-17 1.986425e-16 4.280204e-16 6.565787e-16 7.979281e-16 1.100625e-15
[29,]  9.28125 3.665752e-17 1.915547e-16 4.267280e-16 6.576382e-16 8.064930e-16 1.071432e-15
[30,] 10.40625 3.640059e-17 2.052833e-16 4.256691e-16 6.571531e-16 8.305205e-16 1.187500e-15
[31,] 14.53125 2.695020e-17 1.921481e-16 4.254064e-16 6.464207e-16 7.881546e-16 1.089137e-15
[32,]  9.65625 3.548581e-17 2.010661e-16 4.247601e-16 6.614814e-16 8.128390e-16 1.129341e-15
[33,] 14.03125 2.914324e-17 1.920468e-16 4.244952e-16 6.668548e-16 8.017761e-16 1.160369e-15
[34,] 13.78125 4.337613e-17 1.997844e-16 4.237318e-16 6.668899e-16 8.064630e-16 1.015416e-15
[35,] 11.53125 4.261088e-17 2.071314e-16 4.230143e-16 6.572290e-16 8.539094e-16 1.281702e-15
[36,] 15.28125 2.828213e-17 1.914919e-16 4.221052e-16 6.620637e-16 7.928603e-16 1.028134e-15
[37,] 12.90625 3.453294e-17 2.058897e-16 4.214771e-16 6.448743e-16 7.991928e-16 1.203923e-15
[38,] 12.65625 3.327868e-17 2.022907e-16 4.203595e-16 6.541717e-16 7.938786e-16 1.159096e-15
[39,] 13.03125 3.217340e-17 2.002142e-16 4.202348e-16 6.567995e-16 8.017598e-16 1.079404e-15
[40,] 13.90625 3.154934e-17 1.892979e-16 4.197169e-16 6.645068e-16 8.082360e-16 1.123478e-15
[41,]  9.40625 4.816564e-17 1.978657e-16 4.193184e-16 6.477997e-16 8.099672e-16 1.185520e-15
[42,] 12.78125 4.641259e-17 2.010458e-16 4.188343e-16 6.620586e-16 8.268758e-16 1.108666e-15
[43,] 15.78125 3.534388e-17 1.979963e-16 4.163894e-16 6.348631e-16 8.055629e-16 1.037695e-15
[44,] 14.78125 2.973869e-17 1.783082e-16 4.159072e-16 6.600918e-16 8.180505e-16 1.101258e-15
[45,]  9.15625 3.532054e-17 2.001469e-16 4.144972e-16 6.641739e-16 8.333346e-16 1.161535e-15
[46,] 14.28125 3.729690e-17 2.051280e-16 4.141394e-16 6.420828e-16 8.014932e-16 1.061262e-15
[47,] 10.78125 2.793169e-17 1.640889e-16 4.130271e-16 6.354924e-16 7.865265e-16 1.232986e-15
[48,] 12.28125 3.500407e-17 2.057835e-16 4.120828e-16 6.535214e-16 8.014065e-16 1.178198e-15
[49,] 11.40625 2.859480e-17 1.902818e-16 4.091972e-16 6.344905e-16 8.045298e-16 1.064704e-15
[50,] 10.53125 3.203451e-17 1.838540e-16 4.087420e-16 6.574679e-16 8.164804e-16 1.214193e-15
[51,] 13.40625 3.593774e-17 2.025847e-16 4.073746e-16 6.372493e-16 7.880560e-16 1.015534e-15
[52,] 14.65625 2.595704e-17 1.982885e-16 4.045132e-16 6.367027e-16 7.915563e-16 1.075244e-15
[53,] 12.03125 3.031214e-17 1.812565e-16 4.038055e-16 6.278804e-16 7.811765e-16 1.045932e-15
[54,] 15.03125 3.786686e-17 1.870550e-16 4.008973e-16 6.544609e-16 7.940192e-16 1.148569e-15
[55,]  9.53125 2.514509e-17 1.709193e-16 3.972882e-16 6.359508e-16 7.768724e-16 1.096575e-15
[56,] 15.15625 2.784399e-17 1.711396e-16 3.947257e-16 6.364446e-16 7.832298e-16 1.032463e-15
[57,] 10.15625 1.861819e-17 1.137710e-16 3.593983e-16 6.223346e-16 7.745467e-16 1.137758e-15
[58,] 13.28125 2.128289e-17 1.180957e-16 3.315593e-16 5.960430e-16 7.713776e-16 1.066117e-15
[59,]  6.28125 2.373798e-17 1.240821e-16 2.953215e-16 5.183405e-16 6.519983e-16 9.178537e-16
[60,]  6.53125 2.261652e-17 1.134488e-16 2.372971e-16 3.841038e-16 5.216101e-16 9.013567e-16
[61,]  6.03125 2.166066e-17 1.071175e-16 2.285122e-16 4.138395e-16 5.639309e-16 8.601669e-16
[62,]  7.40625 2.018394e-17 1.041420e-16 2.270119e-16 3.676625e-16 4.945741e-16 7.503839e-16
[63,]  7.15625 2.082297e-17 1.123879e-16 2.267798e-16 3.399366e-16 4.388698e-16 7.180046e-16
[64,]  6.90625 2.859889e-17 1.100437e-16 2.223190e-16 3.307946e-16 4.303613e-16 6.866955e-16
[65,]  8.53125 1.816940e-17 1.068312e-16 2.222083e-16 3.370688e-16 4.257588e-16 6.867601e-16
[66,]  8.40625 2.492418e-17 1.057942e-16 2.212861e-16 3.571108e-16 4.654783e-16 7.328423e-16
[67,]  8.78125 2.073643e-17 1.077757e-16 2.204017e-16 3.449525e-16 4.422913e-16 7.872238e-16
[68,]  8.15625 1.602352e-17 9.785526e-17 2.174716e-16 3.376833e-16 4.592346e-16 7.238405e-16
[69,]  7.53125 2.412168e-17 1.016681e-16 2.169689e-16 3.520386e-16 4.414869e-16 8.018646e-16
[70,]  7.78125 2.157755e-17 1.051076e-16 2.147218e-16 3.351922e-16 4.233797e-16 7.166970e-16
[71,]  6.78125 1.940116e-17 9.827782e-17 2.120133e-16 3.533999e-16 4.766353e-16 7.952774e-16
[72,]  7.90625 2.255016e-17 1.039541e-16 2.112024e-16 3.258873e-16 4.197904e-16 6.163939e-16
[73,]  7.28125 1.647475e-17 9.933971e-17 2.099438e-16 3.355351e-16 4.240930e-16 7.467793e-16
[74,]  6.65625 1.797278e-17 9.737020e-17 2.083182e-16 3.363630e-16 4.246787e-16 6.327588e-16
[75,]  8.90625 2.143521e-17 9.959370e-17 2.069356e-16 3.334586e-16 4.216428e-16 6.287733e-16
[76,]  8.28125 2.002590e-17 1.021602e-16 2.055662e-16 3.300786e-16 4.338119e-16 7.138743e-16
[77,]  8.03125 1.870518e-17 1.001423e-16 2.023863e-16 3.317792e-16 4.142875e-16 5.982200e-16
[78,]  8.65625 2.139552e-17 9.931911e-17 2.019536e-16 3.287046e-16 4.252925e-16 7.261358e-16
[79,]  7.65625 1.814795e-17 9.332480e-17 1.988673e-16 3.348150e-16 4.414474e-16 6.842426e-16
[80,]  6.40625 2.187582e-17 9.824436e-17 1.979501e-16 3.321295e-16 4.255058e-16 6.508132e-16
[81,]  7.03125 1.317255e-17 7.770771e-17 1.877866e-16 3.061724e-16 4.287076e-16 6.225555e-16
> ## For R 4.3.3 :
> ##      shs             5%          25%          50%          75%          90%         100%
> ## 14.53125   9.410630e-15 9.815160e-15 1.023178e-14 1.065722e-14 1.092232e-14 1.138372e-14
> ## 15.03125   8.265317e-15 8.702900e-15 9.086072e-15 9.506915e-15 9.756106e-15 1.007928e-14
> ## 15.90625   6.799207e-15 7.137733e-15 7.611360e-15 8.057580e-15 8.343992e-15 8.670817e-15
> ## 13.53125   6.716182e-15 7.103502e-15 7.566360e-15 8.004966e-15 8.276645e-15 8.630780e-15
> ## 15.65625   6.031124e-15 6.389848e-15 6.803347e-15 7.261310e-15 7.527491e-15 8.031559e-15
> ## ..........
> ## ..........
> 
> myRversion <- paste(R.version.string, "--", osVersion)
> if((mach <- Sys.info()[["machine"]]) != "x86_64")
+     myRversion <- paste0(myRversion, "_", mach)
> if(!capabilities("long.double"))
+     myRversion <- paste0(myRversion, "_no_LDbl")
> myRversion
[1] "R Under development (unstable) (2024-11-15 r87338 ucrt) -- Windows Server 2022 x64 (build 20348)_x86-64"
> 
> rngP <- function(y, M = 1e-14) { yr <- range(y); if(yr[2] < M) yr[2] <- M; yr }
> 
> boxplot(abs19(m.relE), at = shs, log="y", ylim = c(7e-17, rngP(abs(m.relE))[2]), yaxt="n")
> eaxis(2); drawEps.h(); mtext(myRversion, adj=1, cex=3/4)
> 
> matplot(shs, qrelE, type="l", log="y", yaxt="n", ylim = rngP(qrelE))
> title("|relErr( dgamma(x, sh) |   for  x / sh  in [.5, 1.25]")
> eaxis(2); drawEps.h(); mtext(myRversion, adj=1, cex=3/4)
> 
> ## take *one* of these:
> plot(abs(m.relE[, shs == 14.53125]), type="l", log="y", ylim = c(1e-16, 1.5e-14))
> drawEps.h()
> 
> sh <- 14.53125
> stopifnot(identical(sh, 465 / 32))
> 
> x14.5 <- sh2x_gamma(sh, nx = 21) # 21 points
> xM <- mpfr(x14.5, 512)
> dg1 <- stats::dgamma(x14.5, sh)
> dgM <- Rmpfr::dgamma(xM,    sh)
> cbind(x14.5, relE = asNumeric(relErrV(dgM, dg1))) # very "constant" ~=~  - 1e-14
          x14.5          relE
 [1,]  8.265625  6.785869e-17
 [2,]  8.810547  3.161981e-17
 [3,]  9.355469  1.052582e-16
 [4,]  9.900391 -6.737270e-17
 [5,] 10.445312  4.426931e-17
 [6,] 10.990234 -1.488916e-16
 [7,] 11.535156  5.594019e-17
 [8,] 12.080078 -9.633252e-20
 [9,] 12.625000 -1.152292e-16
[10,] 13.169922  1.130550e-16
[11,] 13.714844 -1.834899e-17
[12,] 14.259766 -1.045007e-16
[13,] 14.623047  1.190495e-17
[14,] 14.804688 -5.127820e-17
[15,] 15.077148 -3.645089e-17
[16,] 15.304199  7.765576e-17
[17,] 15.349609 -3.751054e-17
[18,] 15.417725 -1.539056e-16
[19,] 15.474487 -3.722260e-17
[20,] 15.502869  4.550291e-16
[21,] 15.517059  1.310339e-16
[22,] 15.524155  1.352235e-17
[23,] 15.527702  1.097760e-16
[24,] 15.531250  1.280574e-16
[25,] 15.534798  1.042788e-16
[26,] 15.538345  1.010527e-16
[27,] 15.545441 -6.785626e-17
[28,] 15.559631  7.267617e-17
[29,] 15.588013 -3.031525e-16
[30,] 15.644775 -7.823027e-17
[31,] 15.758301 -2.498776e-17
[32,] 15.894531 -2.327030e-16
[33,] 15.985352 -2.231037e-17
[34,] 16.439453  1.310205e-16
[35,] 16.984375  6.502887e-17
[36,] 17.347656 -9.150441e-17
[37,] 17.529297  4.144616e-17
[38,] 18.074219  7.840489e-16
[39,] 18.619141  1.156866e-16
[40,] 19.164062  7.993767e-17
>                                         #            in R-devel   around  1e-16 !!
> ## try easier x:
> sh <- 14.53125 ; stopifnot(identical(sh, 465 / 32))
> x0 <- 1/4 + 8:20
> xM <- mpfr(x0, 512)
> dg1 <- stats::dgamma(x0, sh)
> dgM <- Rmpfr::dgamma(xM, sh)
> relE <- asNumeric(relErrV(dgM, dg1))
> signif(cbind(x0, relE, abs(relE)), 4) # R <= 4.3.*: very "constant" ~=~  - 1e-14
         x0       relE          
 [1,]  8.25  1.276e-16 1.276e-16
 [2,]  9.25  1.294e-16 1.294e-16
 [3,] 10.25 -1.408e-16 1.408e-16
 [4,] 11.25 -2.108e-17 2.108e-17
 [5,] 12.25 -1.306e-17 1.306e-17
 [6,] 13.25  1.464e-16 1.464e-16
 [7,] 14.25 -8.908e-17 8.908e-17
 [8,] 15.25 -5.852e-18 5.852e-18
 [9,] 16.25 -3.029e-17 3.029e-17
[10,] 17.25  1.900e-16 1.900e-16
[11,] 18.25 -1.511e-17 1.511e-17
[12,] 19.25 -5.779e-17 5.779e-17
[13,] 20.25  1.848e-16 1.848e-16
> ## R-devel:                   |  no-long-double  == *same* numbers
> ##    x0       relE           |       relE
> ##  8.25  1.276e-16 1.276e-16 |  1.276e-16 1.276e-16
> ##  9.25  1.294e-16 1.294e-16 |  1.294e-16 1.294e-16
> ## 10.25 -1.408e-16 1.408e-16 | -1.408e-16 1.408e-16
> ## 11.25 -2.108e-17 2.108e-17 | -2.108e-17 2.108e-17
> ## 12.25 -1.306e-17 1.306e-17 | -1.306e-17 1.306e-17
> ## 13.25  1.464e-16 1.464e-16 |  1.464e-16 1.464e-16
> ## 14.25 -8.908e-17 8.908e-17 | -8.908e-17 8.908e-17
> ## 15.25 -5.852e-18 5.852e-18 | -5.852e-18 5.852e-18
> ## 16.25 -3.029e-17 3.029e-17 | -3.029e-17 3.029e-17
> ## 17.25  1.900e-16 1.900e-16 |  1.900e-16 1.900e-16
> ## 18.25 -1.511e-17 1.511e-17 | -1.511e-17 1.511e-17
> ## 19.25 -5.779e-17 5.779e-17 | -5.779e-17 5.779e-17
> ## 20.25  1.848e-16 1.848e-16 |  1.848e-16 1.848e-16
> 
> if(getRversion() >= "4.4.0") # *not* true for  R <= 4.3.3 :
+     stopifnot(abs(relE) < 4e-16) # seen max = 1.900e-16
> 
> cat('Time elapsed: ', proc.time(),'\n') # "stats"
Time elapsed:  15.79 0.2 15.98 NA NA 
> if(!interactive()) warnings()
Warning messages:
1: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
2: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
3: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
4: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
5: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
6: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
7: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
8: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
9: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
10: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
11: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
12: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
13: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
14: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
15: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
16: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
17: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
18: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
19: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
20: In xy.coords(x, y, xlabel, ylabel, log) :
  5 x values <= 0 omitted from logarithmic plot
> 
> proc.time()
   user  system elapsed 
  15.79    0.20   15.98